第五章 结构动力学中常用的数值解法1

第五章结构动力学中常用的数值解法§5.1概述数值分析技术为结构的动态分析提供了有力的保障，为工程结构在各种复杂的动力学环境下的模拟和仿真提供了有效工具。

工程结构的动态分析主要包括两个方面：结构的动态特性分析和结构动态响应分析标准特征值问题和广义特征值问题1 雅可比方法（Jacobi）、2.Rayleigh-Ritz3.子空间迭代法4. 行列式搜索法行列式搜索法是求解大型特征值问题的另一种方法。

它的特点是综合运用多项式加速割线迭代，移轴向量逆迭代，Sturm序列的性质以及Gram－Schmidt正交化过程，直接计算所需要的任意特征对，通常是计算最小的部分特征值及相应的特征向量。

因此，它是一种计算部分特征对的特殊求解方法。

此方法具有计算速度快，精度高，灵活等优点。

nczos法Lanczos方法目前被认为是求解大型矩阵特征值问题的最有效方法，与子空间迭代法相比，其计算量要少得多。

响应数值分析：1.中心差分法2.Wilson -θ法3.Newmark 法响应求解方法的选择取决的因素有：载荷、结构、精度要求、非线性影响程度、方法的稳定性等。

综合各方面的因素，比较、权衡，才能判定所应采取的方法；有时为了互相验证，也可以同时采取两种以上的方法来处理动响应分析对于载荷，一般分为波传导载荷与惯性载荷。

对结构过于复杂的情况，宜采用直接积分法，结构较简单的情况可采用模态迭加法。

对精度要求较低的初步设计阶段，可采用取少数模态的模态迭加法。

对精度要求较高的最后设计阶段，宜采用直接积分法§ 5.2 求解系统固有频率主振型的近似解法1.邓柯利法：是邓柯利首先通过实验方法建立起来的一个计算公式，后来才得到完整的数学证明。

[]M []δ设质量矩阵，柔度矩阵为则有{}[][]{}0x M x δ+=1894年邓柯利：提出一种近似计算多圆盘轴横向振动基频的实用方法（偏小）设系统作j 阶主振动,则有：2()2{}{}sin {}j j j j x A t x ωωω=-=-代入得特征方程：21([][][]){}0jM I x δω-=有111112*********2222222112221101n njn njn n n n nn nn jm m m m m m m m m δδδωδδδωδδδω--=-假设质量矩阵为对角阵，展开得：1111222222(1)11()nn nn n n jjm m m δδδωω--++++=根据多项式的根与系数之间的关系21jω211ω22211nωω的n 个根，之和为1111222222212111nn nnnm m m δδδωωω+++=+++由于二阶频率往往比基频高得多22221111n ωωω111122222111nnn nn ii ii i m m m m δδδδω==+++=∑22211n ωω得忽略111nii iii mωδ==∑ii ii m ω表示仅有质量单独存在时（原多自由度系统变成单自由度系统）的固有频率1ii ii ii ii iik m m ωδ==设2222111221111nnωωωω=+++如例题1m 2m 3m 3331122339169768768768l l l EI EI EI δδδ===22113331117689EI m l m ωωδ===⨯222322176816EI m l m ωδ==⨯333211921634768768768l m l m l mEI EI EIω⨯=+=134.752EImlω=134.933EImlω=精确解2.雅可比（Jacobi ）法求特征方程[]A 设为对称阵，[]{}{}A x x λ=12[][][][](,,)Tn S A S D diag d d d ==即可断定[D]的n 个对角元素就是[A]的n 个特征值，而[S]的第i 列就是[D]中第i 个对角元素所对应的特征向量，[S]为坐标变换矩阵。

pq a [(,,)]S p q θ在[A]中非对角线元素中选取一个绝对值最大的元素，设为，利用平面旋转矩阵A 对进行正交变换：1()()()[][(,,)][][(,,)]i i i i A S p q A S p q θθ+=其中cos cos sin sin 0,,pp qq pq qp ij s s s s s i j p qθθθθ===-==≠22pqpp qqa tg a a θ=-(1)(2)()[][][][]n S S S S =()1,ii s i p q =≠用雅可比法求n 阶对称矩阵[A]的特征值和特征向量的步骤pq a n S I =①设为单位矩阵②在A 中选取非对角线元素中绝对值最大的元素pqx a =-1()2pp qq y a a =-22xx yω=+sin 2θω=2sin 2(11)ωθω=+-2cos 1sin θθ=-④Vb 实现pq a ε<③tan 2xyθ=[]11cos sin 11sin cos 11p q p S q θθθθ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第列第列第行第行3.瑞利（Rayleigh ）法已知系统的刚度[K]，质量[M]，并设定系统的j 阶主振型为(){}j A 对于作简谐运动的多自由度系统，其动能T 与势能V 1{}[]{}2T T x M x =1{}[]{}2T V x K x =系统作j 阶主振动时()(){}{}sin()j j j x A t ωϕ=+()(){}{}cos()j j j j x A t ωωϕ=+()()2{}{}sin()j j j j x A t ωωϕ=-+速度及加速度2()()max 1{}[]{}2T j j j T A M A ω=()()max 1{}[]{}2T j j V A K A =()()2()(){}[]{}{}[]{}T T j j j j j A K A A M A ω=()(1,2)j j r ψ=里兹法：是瑞利法的改进将瑞利法使用的单个假设模态改进为若干个独立的假设模态的线性组合()1{}[]{}r j j j A a a ψψ===∑12{}{}Tr a a a a ={}[]{}{}[]{}(){}[]{}{}[]{}T T T T A K A a K a R a A M A a M a ψ==[][][][]T K K ψψ=[][][][]T M M ψψ=瑞利商在真实模态处取驻值0(1,2)j R j r a ∂==∂（参见刘延柱振动力学107页）得2([][]){}0K M a ω-=问题又归结为矩阵的本征值问题，但与原系统的本征值比，矩阵的阶数r 小于原系统的阶数n.()2()[]{}[]{}j j j K A M A ω=4.矩阵迭代法()2()[]{}[]{}j j j K M φωφ=标准化()1()21{}[][]{}j j j K M φφω-=整理得：()()21{}[][]{}[][]j j j M M D φδφδω==对于正定系统：可写为11{}[]{}i i iA D φφ++=迭代210[]121011K k -⎡⎤⎢⎥=--⎢⎥⎢⎥-⎣⎦00[]0000m M m m ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦假设一个初值11111[][]122123K k δ-⎡⎤⎢⎥==⎢⎥⎢⎥⎣⎦111[][][]122123m D M k δ⎡⎤⎢⎥==⎢⎥⎢⎥⎣⎦解：取初值：01{}11φ⎧⎫⎪⎪=⎨⎬⎪⎪⎩⎭110111113{}[]{}1221 1.6612312m m A D k k φφ⎡⎤⎧⎫⎧⎫⎪⎪⎪⎪⎢⎥===⎨⎬⎨⎬⎢⎥⎪⎪⎪⎪⎢⎥⎣⎦⎩⎭⎩⎭221114.66{}[]{}[]1.66 1.782 2.21m A D D k φφ⎧⎫⎧⎫⎪⎪⎪⎪===⎨⎬⎨⎬⎪⎪⎪⎪⎩⎭⎩⎭33214.99{}[]{} 1.7992.24m A D k φφ⎧⎫⎪⎪==⎨⎬⎪⎪⎩⎭44315.039{}[]{} 1.82.24m A D k φφ⎧⎫⎪⎪==⎨⎬⎪⎪⎩⎭振型值趋于稳定。

211 5.039mk ω=10.445k mω=迭代终止：{}(1)11.82.24φ⎧⎫⎪⎪=⎨⎬⎪⎪⎩⎭几种数值算法的比较：根据瑞利商的性质，原则上可用瑞利商计算任意阶固有频率，但由于高阶主振型很难合理假设，所以瑞利商一般用于求基频，瑞利商求的基频是真实值的上限，这是因为假设的一阶主振型与真实振型的偏差，相当于对系统附加了某些约束，从而提高了系统的刚度，使基频有所提高。

矩阵迭代法只能求前几阶频率，求较高阶的频率和模态需要把前几阶的模态剔除才能收敛。

邓柯莱公式建立的是质量矩阵为对角阵的情况，而且第二阶及第二阶以上的固有频率通常远大于基频，邓柯莱公式计算出的基频显然是精确值的下限。

5.子空间迭代法子空间迭代法实质就是对一组试验向量反复地使用里兹法和矩阵迭代法，它将矩阵迭代法每次迭代一个假设模态，发展为同时迭代系统的前r 阶假设模态，因而提高了计算效率，迭代过程中各阶假设模态的正交性由里兹法保证。

(){}(1,2)j j r ψ=(){}(1,2)j j r ψ=n r ⨯设系统的前r 阶模态构成全部n 阶模态所张成的线性空间的一个子空间，任选r 个独立的向量作为子空间的假设模态,组成阶矩阵(1)(2)()[]{}{}{}r ψψψψ⎡⎤=⎣⎦各个假设模态(){}(1,2)j j r ψ=总能表示成真实模态的线性组合。

()()()1{}{}(1,2)n j j i i i a j r ψφ===∑②式()()[]{}{}j j j D φλφ=[]D 代入②式两边乘以阵()()()()()11111[]{}[{}{}](1,2)r n j j i j i i i i i i i r D a a j r λλψλφφλλ==+⎛⎫⎛⎫=+= ⎪ ⎪⎝⎭⎝⎭∑∑如此迭代k 次后()()()1[]{}(){}r k j j ki j i i D a ψλφ==∑第二项包含的高于r 阶模态的成分比第一个求和式更快趋近于零。

{}[]{}()()j j j D λφφ=子空间迭代的步骤0[][]n r A ψ⨯=假设初始模态矩阵(1)(0)[][][]D A ψ=将上式左乘矩阵[D](1)[]ψ(1)[]A 利用各列的线性组合表示子空间基的一次近似(1)(1)[][][]A a ψ=(1)(2)()[]{}{}{}Tr a a a a ⎡⎤=⎣⎦[]a (1)[]A 以为假设模态，进行里兹法计算，以确定系数矩阵(1)(1)[][][][]T K K ψψ=(1)(1)[][][][]T M M ψψ=先作出以下r 阶方阵瑞利商取驻值[][][]([][])[][][]T T a K a R a a M a ψ=[]0[]R a ∂=∂2([][])[]0K M a ω-*=r r {}a 解此本征值，得到个本征值和个本征向量(1,2,)i a i r =(1)(1)[][][]A a ψ=(1)A 代入式中则子空间的一次近似完全确定，构成的各个模态满足正交性条件，至此完成第一次迭代，(1)[]A (0)[]A 将代替进行第二次迭代(2)(1)[][][]D A ψ=(2)(2)[][][]A a ψ=得：100[]012002M m ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦210[]132022K k -⎡⎤⎢⎥=--⎢⎥⎢⎥-⎣⎦1112[][][][][]124125m D M K M k δ-⎡⎤⎢⎥===⎢⎥⎢⎥⎣⎦例取前二阶假设模态，归一化后作为子空间基的零次近似00.52[]0.9111A -⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦01120.52 3.41[][]1240.91 6.30125117.31m m D A k k --⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥=-=⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦(1)0.4651[]0.863011ψ-⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦代入式得：0.41230.2054[]0.20544K k ⎡⎤=⎢⎥⎣⎦ 2.9617 1.5342[] 1.53423M m ⎡⎤=⎢⎥⎣⎦设2k m γω=代入得：20.4123 2.96170.2054 1.53421.607(7.7495 4.0643)00.2054 1.534243γγγγγγ--=-+=--127.18370.5658γγ==解出本征值：120.3731 1.3295k k m mωω==438.70.519611a -⎡⎤=⎢⎥⎣⎦(1)0.4625 2.58540.86100.933411A -⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦得到前两阶固有频率的一次近似值及对应的系数矩阵(1)1120.4625 2.58541240.86100.9334125113.3235 1.51886.18450.45227.18450.5478m DA k m k -⎡⎤⎡⎤⎢⎥⎢⎥=-⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦-⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦(2)0.4626 2.77250.86080.825511ψ-⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦0.41130.000956[]0.00095618.142K k ⎡⎤=⎢⎥⎣⎦ 2.9550.00685[]0.0068510.268M m ⎡⎤=⎢⎥⎣⎦代入得：本征方程为：20.4113 2.95500.0009560.006850.0009560.0068518.14210.3687.462(7.756 4.106)0γγγγγγ----=-+=解出本征值127.18450.5715γγ==前二阶固有频率的二次近似值，更接近真实值120.3730 1.3228k k m mωω==作业：5-1（1）-（5）。