结构动力学中的常用数值方法精

1、第五章 结构动力学中的常用数值方法5. 1.结构动力响应的数值算法.M x ex kx = F (t)彳x(0) = aoI.x(0) =Vo当e为比例阻尼、线性问题> 模态叠加最常用。但当 C无法解耦，有非线性存在，有冲击作用(激起高阶模态，此时模态叠加法中的高阶模态不可以忽略)。此时就要借助数值积分方法，在结构动力学问题中， 有一类方法称为直接积分方法最为常用。所识直接是为模态叠加法相对照来说， 模态叠加法在求解之前， 需要对原方程进行解耦处理，而本节的方法不用作解耦的处理，直接求解。(由以力学，工程中的力学问题为主要研究对象的学者发展出来的)中心差分法的解题步骤12(氏)1，a1

2、：2：ta2=2ao1.初始值计算(1)形成刚度矩阵K，质量矩阵M和阻尼矩阵Co(2)定初始值x0 ,Xo , Xo o(3)选择时间步长At ,使它满足.,：£，并计算(4) 计算 X=Xo x xo2a2a°(5) 形成等效质量阵 M -aoM - a1C(6) 对M阵进行三角分解M =LDLt2.对每一时间步长(1) 计算时刻t的等效载荷Qt=Q-( K- 2a M t x ( o a -M 1 a)Q：tX(2) 求解t 氏时刻的位移(LDLT)Xt -. t -Qt(3) 如需要计算时刻t的速度和加速度值，则xt 二a°(Xt t -2冶 x-t) 若系

3、统的质量矩阵和阻尼矩阵为对角阵时，则计算可进一步简化。纽马克法的解题步骤1初始值计算（1） 形成系统刚度矩阵K，质量矩阵M和阻尼矩阵C（2） 定初始值x0， Xo， Xo。（3） 选择时间步长过，参数、二。并计算积分常数y>0. 5 6 启0.25（0.5+Y）2小1畀c1a 2 ， a ， a2 0;.-:t21;.-:t 21YAt 丫a31，a4 一 -1，a5（- -2）3 2：， 4， 5 2 I氏=：t（1 - ），= t（4） 形成等效刚度矩阵KK = K a0MaQ（5）K矩阵进行三角分解K 二 L DtL2.对第一时间步长（1） 计算时刻的等效载荷Qt t =QtM （

4、a6Xta2 xt a3 xt） 。为 a4 xta5Xt）（2）求解t t时刻的位移（LDLXt 吧=Qt：费 t（3）计算t =t时刻的加速度和速度xt . t - a0（xt . t _ xt） _ a2 xt _ a3 xtxt . xt a6 xt a7 xt餉威尔逊法的解题步骤1.初始值计算（1）形成系统刚度矩阵K，质量矩阵M和阻尼矩阵C（2） 定初始值x0， X0，X0。（3） 选择时间步长氏，并计算积分常数v -1.4a° -,a? = 2ai"t6Ata厂万，a4aoa5一 a2a6 =1 一3 ,a7o :ta8：t(4)形成等效刚度KK = K a0

5、M aQ(5)将等效刚度 K进行三角分解K 二 l dtl2对每一个时间步长(1) 计算t 氏时刻的等效载荷Rt 站t =Qtv(Qt . t Qt) M (a0xt a2 xt 2xt)C(a1xt 2xt a3 xt)(2) 求解t 过时刻的位移(LDLXtY t 二 Rt P(3) 计算在t 过时刻的加速度、速度和位移xt -.x = a4(x： - Xt) a 5Xt a 6xt5. 2结构动力响应数值算法性能分析对公式(5.1)描述的线性系统结构动力学问题，已经有证明对整个多自由度的积分，等价于将模态分解后对单自由度的积分的结果进行模态叠加，因此可以通过对单自由度问题的分析，来说明算

6、法的特性，其中阻尼均假设为比例阻尼，这样，模态分解后的单自由度结构动力学方程为：x 2 X 2X = f (t)(5-29)以下算法的性能分析，均将算法用于这个方程。5. 2. 1算法用于结构动力学方程的有限差分表示将数值计算方法应用于(5-29),即分别在相邻的不同时刻应用算法可得如下一般形式yk 1 二 Ayk Lk( 5-30)A为放大矩阵或称逼近算子，Lk为载荷逼近算子。XkJ,*xk 1，X k ,xk _m例如将Newmak方法应用于方程5-29)有:(5-33)2、2 2h2、.，1 2h 'Adh221(1 一2眈22-h(Y)oh1 h(1 _2、;厂矩阵A的特征多项

7、式为2-det(A- I)仝2 -2A A2(5-34)其中A1, A为该矩阵的两个特征向量，分别为矩阵的迹的一半和矩阵的行列式1 1A - traced (Ah A22)(5-35)对Newmak方法有:1 (2 一1)VA?二 det A A11A22 - A12'1,2 = A1 二A1 'A21 (2 一2)一 -)2 A2 _2(5-36)(5-37)(5-38)其中h为时间步长，寫壬:h,D =1 2'S 。 Newmak方法放大矩阵的规模是二维的,因此特征值也只有两个，可以根据它们进行分析。有的算法放大矩阵是三维的，例如Wilson-二方法，在无阻尼情况下

8、放大矩阵为:诃-1)069.-6豹 2力(6 V2"32)h2(6 K -2-斧2)e(02 (日2-3)+6)h(6日+e302 -3-3日202 / 2)( 5-39)6h內6日+。2日3 _3日2°2D =(门勺26门放大矩阵A的特征多项式为:(5 40)32idet(A-，l)=2A1.；” A2 - A3 - 0其中A1, A2, A3为该矩阵的三个特征向量，分别为矩阵的迹的一半、各阶主子式的和以及 矩阵的行列式，对 Wilson-二方法有A18c6十3°2C2 如3。2 日Al =222cg 0 +6)22322A402 +阻_6& 日 2刊8

9、8_12< A2 = ( 5-41 )+6)_6 +6H _02 d。2 日 2 +°2 日 3 +3° 2日A3 2 2ieg e +6)此外，在几个不同时刻应用数值算法，然后将方程中的速度和加速度项消去，可得数值算法关于位移的差分方程，例如Newmak方法，有1 2(1 2 小、)Xn .1 一21 (2 一1)(、._? _?Xn1 21 (2 -2)门()2Xn=：O( 5-42)2很显然，其特征方程与其放大矩阵A的特征方程是相同的，使用关于位移的线性多步方式和放大矩阵来说明算法性能是一样的，只不过各有方便之处。5. 2 . 2算法的稳定性分析设人，i =1,

10、2m为放大矩阵a的特征值，则P = maxi定义为A的谱半径，若特征值互异，则 乞1的算法是稳定的，但若有重特征根，则要求 r ： 1。如果算法的稳定性要求对步长的选取有限制，称算法是有条件稳定的，反之为无条件稳定的。放大矩阵的谱半径小于等于1成立的充分条件是1 _2厲 A2 _01 2A1 A2 _0(5-43)1 -A2 _0对3 3的放大矩阵彳-24 +2A2 - A3 工032舛A2 +3A3 303*2A1 - A2 -3民 30( 5-44)1+2A, +A2 +A3 X0A2 +A3(2A _AJ 王0上两式是关于算法自由参数/1的不等式，由它可以判断算法是否无条件稳定，若不是，

11、将给出稳定条件。例5 1分析Newmak方法、Wilson-二方法的稳定性解： 将(5-38)代入(5-43 )有2 八"2(_) _0门2（2一（1一2）门一1 乞°显然，当1-2- -2(5-45)算法无条件稳定。当小辽且(一1) 2( -1)22 “C二22262(5-46)算法稳定，但为条件稳定，其中'Jc为临界采样频率。由于（5-43）式仅仅是充分条件，所以可进一步按照稳定性的定义得到5.1.2节叙述的无条件稳定条件。对Wilson-二方法，将(5-41)代入(5-44 )得602 狂0602(2日-1) >0«12 +02(1+6日2 6

12、日)H0(4日3 +1 6日2)02 +24C-12 3042Q (2c +1 3c) K0(5-47)容易看出，其中第一，二，五不等式恒成立，对第三，四不等式若希望对任意的 则有：11均成立,求解上述不等式得实际使用中通常选取-=1.45.2.3算法的相容性和收敛性直接积分算法的相容性、一 1托日2 一6日0324日一6日+1启0存十7(5-48)收敛性分析同样要使用其位移型的差分方程,或对应的单步(5-49)(5-50)多值形式。在算法（5-30）式中，用精确解代替近似解，即可得到局部截断误差表达式，用符号e（tk）表示y(tk i)=Ay(tQ Lk ' he(tQ局部截断误差表

13、达式用放大矩阵的特征量以最常用的线性三步法为例可表示为e(tk)鬥x(tkh) -2Ax(tk)A2X(tk -h) -A3X(tk -2h)/h2其中Al, A2, As分别为对应的3 3的放大矩阵的三个特征向量，然后将x(tk h), x(tk 2h)在tk点进行泰勒展开，然后利用运动平衡方程化简即可。若局部截断误差表达式为步长的s(s>0)阶小量，则称算法是 s阶相容的。对2 2的放大矩阵，可仿照上述步骤，来验证算法的相容性。在经典的数值算法收敛性分析理论中，一个重要的结论就是相容加稳定等于收敛，其相容的阶数就是算法的精度阶。收敛性的含义也是当时间步长趋于零，算法的数值解趋于精确解

14、。对直接积分算法该定理同样可以证明是成立的。例5 2分析Newmak法的相容性和精度解：其局部误差仿照(5 50)式得:e(tk)匚x(tkh) -2駆怨)A2(tk -h)(5-51)也可以由(5-42)是直接求得，即将x(tk h),x(tk -h)在tk点泰勒展开，并注意到在tk时刻的运动方程有：e(tk) =( -2 (-2、-：)-( -2)x(3)h2 6 2(5-52)【(_' 厂占)2 "召、7l2x(4)h2 o(h3)261212241显然，当物理阻尼为零时，选择=2算法是二阶的，即截断误差是步长的二阶小量。物理1阻尼的存在，使算法精度降了一阶，但若同时选

15、择:.二丄，算法精度仍然是二阶的，一般称6为Newmak线加速度法。显然 Newmak方法中有两个参数待定，每种特定的选取都是一个 特定的算法，最常用的几个算法见表5-1表5-1:常用的Newmak族直接积分算法16方法名称稳定条件无阻尼问 题精度阶有阻尼问题精度阶类型1/21/4平均加速度方法(梯形法)无条件21隐式1/21/6线性加速度方法门 _叫=2. 3 =34622隐式1/20中心差分方法C221显式如果在一个时间步内需要求解一个隐式的方程组， 则称算法是隐式的，反之不需要求解方程， 直接计算即可得到下一时刻的值， 则称算法是显示的。 从5.1节的Newmak方法的计算步骤 可以看出

16、，这类方法是隐式的， 但对于中心差分方法， 若质量矩阵和阻尼矩阵都是对角矩阵就可以显示地计算。 显然显示方法计算量要小得多。读者可自行分析 Wilson-二方法的精度,不难分析，无论是无阻尼还是有阻尼其精度都是2阶的，它也是隐式方法。524算法耗散和弥散特性算法的精度，在小步长的情况下可以通过局部截断误差分析来说明比较，但是，在实际计算过程中，步长的选取可能不是很小，此时如何来度量算法的计算精度，当然可以针对有解析解的问题进行大量的数值计算，将数值解与解析解进行比较来分析算法的计算精度。理论上还可以通过数值耗散(disspation)和弥散(dispersion)来辅助度量与分析，为引出这两个

17、概念的含义，我们仍然以单自由度有阻尼自由振动问题为例，该问题的解析解为：x(t) =e-(C| cosco dt + c2 sinodt)(5-53)当直接积分算法用于这样的问题，前小节已经讲述过它可以写成形如(5-42)的关于位移的有限差分形式，就是可以得到一个关于位移的有限差分方程，对于一个收敛的且有一定精度的算法，这个差分方程通常有一对共扼复根，其中-=.2，该两根称为主根，其它根称为寄生根(spurious roots)。解的一般形式可(5-54)入,2 =a ±bi(5-56)写为mXn =eYn(C1 COS石dtn +C2Sin 空dtn)+瓦(5-55)i=3式中的称

18、为算法阻尼比，当有物理阻尼存在时，它还包括了物理阻尼的影响，称为算(5-53)式形式是相同的,法频率，对应的T =2二/称为算法周期。可以看到上式前两项与这给了我们与精确解进行比较的可能，如果1)寄生根的影响较小，即 人i=3*m « '-1,2。2)解表达式中得常数C1,C2与C1,C2差别不太大。这样，我们就可以通过比较不同算法的算法阻尼比和相对周期误差，来比较算法的计算精度。(T -T)/T，它会使得数值解此时，对无阻尼问题，可以很明显看到算法阻尼比会使得数值解曲线的幅值与解析解相比要 降低而产生振幅衰减， 这就是所谓的算法的数值耗散。同时，不同的算法的算法周期与精确

19、的周期会有一定的误差，这个误差一般用相对周期误差来表示 曲线上产生周期的延长或缩短，即所谓的数值弥散。实际分析时，可以首先通过求解放大矩阵的特征方程得到特征方程的主根和寄生根，若主根可以表示为：并注意到式(5-54)有(5-57)(5-58)(5-59)d 二 arctan(b/a)-=arctanQ/ a) /. 1 -2禺=0(门2)血 o(h)v°(5-66)v1 =o(l】)v° - o(1)X0(5-67)2 2(5-60)Ln(a2 b2)对结构动力学问题， 一般总希望算法在低频段有较小的耗散和弥散，而且在低频段，当；0时，可以很方便地获得它们的近似解析表达式。

20、在高频段算法的耗散特性，用谱半径来说明更适合，一般用二一 = lim二(A)来度量算法对高频分量的耗散特性，特别地，当二-=0时，称算法具有高频渐进消去特性，即当算法计算一步以后高频极限完全地被耗散掉，而其它高频分量由高到低渐进地被耗散。由前面的叙述可以看到，算法放大矩阵的特征向量AI，A2或A1, A2, A3决定了算法对应特征方程的根，也就决定了算法的稳定性，同时确定了谱半径、以及算法耗散和弥散特性。 这些特性有时也称算法的谱特性。 放大矩阵相同的不同算 法，称为互相相同的，放大矩阵不同但特征值相同， 称算法互相相似，或称算法频谱等价的。对于算法的耗散特性，应该说明的是高频耗散特性对实际的

21、结构动响应求解是有益的，因为实际结构进行有限元离散以后计算出的高频行为并不真正代表系统的物理行为，它是结构系统在空间进行有限元离散的结果，是虚假的行为，而不具备高频耗散特性的算法，是将所有频率上的响应全部进行了积分，尽管步长取得相对较大时，高频的积分不准确，这样的计算结果显然与系统实际的反应不一致。但是同时要注意，它在低频也不同程度地引入了数值耗散，这样这些算法就不适合进行长时间的计算，因为长时间以后应该精确计算的低频响应，由于耗散特性的存在，已经被耗散得面目全非，因此，有耗散特性的直接积分方法只适 合计算瞬态的、短时间内的低频动力响应。例5-3分析Newmak族算法频谱特性解：对Newmak

22、族算法来说，当咐时，算法才可能有二阶精度，我们仅讨论这一种情况。此时算法放大矩阵的两个特征量为：fl 1 (-护.T；2A2A22(5-61)1 W；21m2九2 =A<| 土iA2 _A2 =a±ib则谱半径:=. a2 b2 二A2不考虑物理阻尼时，对任意的Q有就是无阻尼时,对有阻尼问题,T -TT11 1=I =Q 纹 1arctan0 ；1：)02 /(1 +(6 -$02)Newmak族算法不存在数值耗散，但有一定的相对周期误差。(5-63)(5-64)耳=1，上 J： "(2)5.2.4算法的超调特性谱半径这个指标对算法性能的影响还需要进一步说明的是，它只

23、决定算法的长期特性,即乞1可以保证随着算法计算步数的增加计算过程是数值稳定的。但对无条件稳定的算法，由于步长大小选择没有限制，一般在满足指定精度的条件下，尽可能取较大的时间步长，对于非零初始条件问题，在计算开始的几步可能会出现初始数据及其误差(如初始位移，速度的测量误差，初始加速度的计算误差)被放大的现象，这称为超调(overshoot)。这种现象是放大矩阵A病态，有较大的条件数而产生的。实际应用时，由于当二0算法是收敛的，不会出现超调。一般为简单起见，只分析当门：时，在计算的第一步是否会出现超调。例5- 4分析Newmak平均加速度法的超调特性为分析简便起见，将 Newmak平均加速度法用于

24、无阻尼自由振动问题，此时其放大矩 阵为：4-Q24+02 -4 0+Q3/ 2A=142424hX12X02V4 24.门2一4曲+03切/2丄 4_2。2V12x0厂 v04在一:时，可得近似等式X1 =o(1)x°, v1 =o()x0 o(1)v0(5-65)其中o(1), o(Q分别表示关于Q的零次和一次关系式。由此可知算法在位移上无超调，但由于初位移的影响，在速度上有关于11线性超调现象。前面提过 Wils on-二方法有很强的超调现象，对无阻尼问题，从放大矩阵的各元素的表 达式中，很容易得到(甘_02)°4 +(2日2 _6)°2 +12卜 h(日2

25、_1)02 +6)12+2B2。2X06+日2。2(3V -2d2)"3 -121 -(二2 -3)门2 6X06+护02X1V1V0V0在门 > :时，可得近似等式其中oO2）, o（）分别表示关于Q的二次和一次关系式。由此可知 Wilson-r方法在位移上关于初位移有二次超调，同时关于初始速度有一次超调。在速度上有关于初位移一次超调。另外，也可以用数值计算的方法对指定的初始条件，计算出近似解x1,v1，然后与精确解比较，或计算系统能量范数:En2 2Mvn Kxn2(5-68)然后将Eo与Ei比较。显然，由于直接积分方法适合于短时间的瞬态问题计算，因此超调现象也是必须加以注

26、意的。综上所述，对一个数值积分算法理论上要分析其相容性，稳定性，数值的耗散与弥散 特性，对无条件稳定的算法还要分析其超调特性。这样才可能对算法的本质有深入的了解， 进而指导数值计算结果的解释与分析。此外，由于直接积分方法对结构运动平衡方程进行数值积分的目的在于估计结构真实的动力响应。为了精确地预计结构的动力响应，要求模态分解以后所有的单自由度平衡方程都必须被精确地积分，但在直接积分法中，对所有的方程积分都相当于采用相同的步长h,所以时间步长的选择必须针对系统的最小周期。如果Tn是系统的最小周期的话，则 h选为Tn /n，其中一般 n= 10。对条件稳定的算法，当然还要同时考虑这个选取是否满足算

27、法稳定性的要求。对有条件稳定的算法，要求n ，若n取10,多数的条件稳定算法的稳定-C条件都满足，不难验证表 5 - 1中的条件稳定算法全都满足。但需要注意的是，对于大型、复杂的实际结构，经过有限元离散以后通常都有上万，甚至几十万个自由度，其最大固有频率通常都很大，也就是系统的最小周期非常小，此时，按Tn/10来选取步长就非常小，这会大大增加计算量。而实际工程上只关心较低阶的固有频率，同时结构的响应也主要由若干较低阶的响应构成，因此在计算时高频可以不用精确积分，就积分出那些主要的， 感兴趣的低频响应就可以了。也就是步长可选择为 Tp/10，比Tn/10大Tp/Tn倍。由于实际情况中Tp/Tn可

28、能会非常大，这样条件稳定算法的稳定条件就可能无法满足，而无条件稳定算法对步长的选取就没有稳定性的限制， 因此对于实际的结构动力响应计算，多数都使用无条件稳定算法。5.3矩阵特征值问题及解法5. 3. 1 问题分类(K 一 2M ) = 0K = 2M这原本是广义特征值问题，但可以化为标准特征值问题，前乘M，得(MK 一 3)=1令A二M K ,则有(A- 2I)=这是代数里的标准的矩阵特征值问题，对结构动力学问题而言1M对角 r M K = A对称M对称非对角 一；M K = A不一定对称非对称矩阵的结构特征值问题计算量很大，此时作如下处理：1对M进行Cholesky分解(因为 M对称正定)M

29、二LUL (下三角阵)U二LT (上三角阵)代入方程(K - 2M)，(K -2 LU)前乘L得(L°K - 2U) =0提出U得(L4KU 4 - 2I)U =0令 A =LKU，=L则有(A- 2| ) = 0可以征明A与A有相同的特征值，因为它们是相似的矩阵，但A是对称的_T_因为 A =(L4KL)L4KTL -L4KL - A所以A对称。只不过特征向量有所变化，求出'以后=LJ至于Cholesky分解，可直接由 M(或K)的元素计算得到11 112 1l n2l nn的各非零元素j 4一 2、1/2ljj - (mjj _ '1 jr )i=j+1,j+2.

30、 ,nlij = (mij - ' lir l jr ) / l jjj=1,2,.,n,即可求得L的下三角部分各列元素。又若记r =4上式依次取V11V?20Vn2Vnn则其下三角之i=1,2，.,n行元素，依次为Vi =1/liii 1Vij = -(lir Vrj / lii )j=1,2, ,i-1r =1532特征值、特征向量的一些特性1对角阵、三角阵、块对角阵、块三角阵对角阵和三角阵的特征值就是这些矩阵对角元素的数值。块对角阵和块三角阵的特征值就是这些矩阵的对角线上各个子块的特征值。2. 几种特殊矩阵(1) 实对称矩阵的特征值必为实数，其特征向量也可选为实向量。(2) 反对

31、称矩阵的特征值或者为纯虚数，或者为零。(3) 正定对称矩阵的特征值全部大于零。(4) 对称矩阵对应于不同特征值的特征向量彼此正交。3. 与原矩阵A相关联的几种矩阵设矩阵A的特征值是，其对应的特征向量是 X,则(1) A u( I )阵的特征值是'：!-，特征向量是x;(2) aA阵的特征值是a丸(a式0 )，特征向量x；(3) At是A的转置矩阵，则 At的特征值就是；（4） Am阵的特征值是 m，特征向量仍为x,m为整数；（5） 若A非奇异，则A的逆矩阵存在，记为 A，是A的特征值为1/ ,特征向量 为x;（6）若矩阵B与A相似，即有可逆阵 P存在，使B = PA P则B的特征值也是

32、，特征向量是P Jx4特征值的积与积若矩阵A的特征值为、，2、.，'n，则有n. i2 .、 a”i=1r 2“ 二 det（A）它们可作为校核、估计甚至计算特征值的手段，读者可以自行验证。5. 3. 3特征值问题的基本计算方法由于高于四次的一元代权方程无法求精确解，所以必然要采用迭代方法求计算求近似解。 目前的迭代方法，只不过迭代的方法和技巧不同。已有很多成熟方法。方法选择时主要取决于1）所要求的特征对数目2） K、M的阶，3） K、M的带宽以及是否带状。根据用到的基本关系大致可分四类1. 向量迭代法（幕迭代法Power iteratio n method）基本关系K j二 M根据迭

33、代格式不同又分：正向迭代（幕迭代）逆迭代（反幕迭代）主要用于求特征向量。2. 变换法基本关系：gtk-uI门质量归化为的模态矩阵-I =diag（ 1, 2, 3,., 'n）雅可比迭代（Jacobi）（小模型实对称阵标准特征值问题的全部特征对方法） 广义雅可比迭代（求全部特征对，小模型的广义特征值问题，或者有大量非对角线零 元素和少量对角线零元素问题）豪斯霍尔德（Householder）3.多项式迭代法（不单独使用）基本关系：（ -） =det（K -，M ）显式多项式迭代隐式多项式迭代4斯图姆（Sturm）序列法利用；?（ J = det（K -，M ）的Sturm序列性质来求解。实际工程结构的动力学问题绝大多数使用有限元法方法离散求解。随着计算机存储和速度能力的成倍提高，所建模型越来越精确，几十万、上百万个自由度的计算成为可能。为此发展了一些针对大型特征值问题的方法，它们综合上述典型方法和技巧，常用的有：多项式迭代行列式探索法（带宽较窄的低阶特征值问题。子空间迭代法（Subspace iteration）（与行列式法相似，但计算量小）兰索斯（Lanczos） 法向量迭代 HQRI（Householder QR迭代方法）（针对大型满阵的大