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

1、第五章结构动力学中的常用数值方法5.1.结构动力响应的数值算法.Mx-cx-kx=F(t)x(0)=a。Ix(0)=V。当c为比例阻尼、线性问题T模态叠加最常用。但当C无法解耦，有非线性存在，有冲击作用(激起高阶模态，此时模态叠加法中的高阶模态不可以忽略)。此时就要借助数值积分方法，在结构动力学问题中，有一类方法称为直接积分方法最为常用。所识直接是为模态叠加法相对照来说，模态叠加法在求解之前，需要对原方程进行解耦处理，而本节的方法不用作解耦的处理，直接求解。(由以力学，工程中的力学问题为主要研究对象的学者发展出来的)中心差分法的解题步骤1 .初始值计算(1) 形成刚度矩阵K,质量矩阵M和阻尼矩

2、阵Co(2) 定初始值x0,x。，x。(3) 选择时间步长&,使它满足At<Atcr,并计算a11 .,a22a。2 .t计算x_t.：=x。(5)形成等效质量阵M=a°MaC(6)对M阵进行三角分解TM=LDL2.对每一时间步长(1)计算时刻t的等效载荷Qt=Q(K2aMtx(0a-M1a)p.:tx(2)求解t+&时刻的位移(LDLT)xt=Qt(3)如需要计算时刻t的速度和加速度值，则xta1(xt-.t-x02al2a0t.v)若系统的质量矩阵和阻尼矩阵为对角阵时，则计算可进一步简化。纽马克法的解题步骤1 .初始值计算(1)形成系统刚度矩阵K,质量矩阵M

3、和阻尼矩阵C.(2)定初始值x0,X0,X0。(3)选择时间步长”，参数¥、仃。并计算积分常数_0.50.25(0.5-)2一1一一1a0-_2>a1?a2_、；t2、：t、；t1 ,At,一a31,a41,a5(-2)2、:2-a6=；：t(1一&),a7=.'：t(4)形成等效刚度矩阵KK=Ka0Ma1c(5)K7矩阵进行三角分解K=LDtL2.对第一时间步长(1)计算t十&时刻的等效载荷Qt.t=QtM(a6xta2xta3«t)Cax1tax4tax5t)(2)求解t+&时刻的位移(LDLT)xt=Qt.t.：(3)计算t+&a

4、mp;时刻的加速度和速度.xt：A=a0(xt-t-xt)-a2xt-a3<t.xt:由=xta6xt-a7xt：n威尔逊-9法的解题步骤1.初始值计算(1)形成系统刚度矩阵K,质量矩阵M和阻尼矩阵C.(2)定初始值x0,x0,x0。(3)选择时间步长&,并计算积分常数1-1.46一3c->ai=5a2=2aiu；：t2(-.:t)t'ta0-a2a4=，a5=003a6=1一,ayet."：t2=，a8=26(4)形成等效刚度KK=Ka0Ma1c(5)将等效刚度K进行三角分解一TK=LDL2.对每一个时间步长(1)计算t十次时刻的等效载荷Rt<1、

5、=QtRQt.t-Qt)M(a0xta2xt2xt).-C(a1xt-2xt-a3xt)(2)求解t+&时刻的位移(LDLT)xt.：”Rt不(3)计算在t十At时刻的加速度、速度和位移.xtt=a4(xt.二t.：-xt)a5xta6Kt5.2结构动力响应数值算法性能分析对公式(5.1)描述的线性系统结构动力学问题，已经有证明对整个多自由度的积分，等价于将模态分解后对单自由度的积分的结果进行模态叠加，因此可以通过对单自由度问题的分析，来说明算法的特性，其中阻尼均假设为比例阻尼，这样，模态分解后的单自由度结构动力学方程为：-2x+2-C0xx=f(t)(5-29)以下算法的性能分析，均

6、将算法用于这个方程。5.2.1算法用于结构动力学方程的有限差分表示将数值计算方法应用于(5-29),即分别在相邻的不同时刻应用算法可得如下一般形式yk+=Ayk+Lk(5-30)A为放大矩阵或称逼近算子，Lk为载荷逼近算子。ykk,Xk工,xk_m1I,yk：；1例如将Newmak方法应用于方程(5-29)有:1A二AAd(5-33)At二矩阵A的特征多项式为其中aih2(1-2-).22-_h(1一为0h1-h(12、；)，1-2h(1,)-.2一det(A-I)=,一2AlA2Ai,A为该矩阵的两个特征向量，分别为矩阵的迹的一半和矩阵的行列式11A1=traceA=(A11，A22)22A

7、2=detA=61A22_A12A21'1,2=A1二t,A1-A2Newmak方法有:v121-(2-1)'J(J.)'J24D191-(2一2)-(c._-)'?A2二2其中h为时间步长，C=0h,D=1+2馅C+2。Newmak方法放大矩阵的规模是因此特征值也只有两个，可以根据它们进行分析。有的算法放大矩阵是三维的，例如Wilson-0方法，在无阻尼情况下放大矩阵为:1A=一D：tf-1)Q2+68pfk.2j-6O女(6-2'J3-'J2)”2(12-3)-6)6h函h2(60+e3Q2-2-02Q2)Ih(60+03Q2-3-302Q2

8、/2)60+。冶3-302Q2-6)(5-34)(5-35)(5-36)(5-37)(5-38)二维的，(5-39)D=(，.FT6)1放大矢I阵A的特征多项式为:一det(A-',-I)=，3-2A1，2-A2f.-A3=0(540)其中Ai,A2,A3为该矩阵的三个特征向量，分别为矩阵的迹的一半、各阶主子式的和以及矩阵的行列式，对Wilson-日方法有A118c63I-i2-3,F2一3.，、A2A3222c(？6)4.12Qi"3一6ii212181,12一引.匕2222232-66二_，一3,1v3,1v(5-41)c(i12u26)此外，在几个不同时刻应用数值算法，

9、然后将方程中的速度和加速度项消去，可得数值算法关于位移的差分方程，例如Newmak方法，有.2.12(12一、口再1_21(2_1)-J(;.-).Jxn2412出1+(2Y-2)£Q+(5-y+2)Q2Xn=0(5-42)很显然，其特征方程与其放大矩阵A的特征方程是相同的，使用关于位移的线性多步方式和放大矩阵来说明算法性能是一样的，只不过各有方便之处。5.2.2算法的稳定性分析设九，i=1,2m为放大矩阵A的特征值，则P=max兀定义为A的谱半径，若P<1。如果算法的稳定特征值互异，则PE1的算法是稳定的，但若有重特征根，则要求性要求对步长的选取有限制，称算法是有条件稳定的，

10、反之为无条件稳定的。放大矩阵的谱半径小于等于1成立的充分条件是1-2Al-A2_0彳1+2A1+A2之0(5-43)|1-A2_0对3M3的放大矩阵1-2A1+2A2-A3>03-2Al-A2+3A3>013+2AlA23A3>0(5-44)1+2A1+A2+A3之01-A2-A3(2A1-A3),0上两式是关于算法自由参数Uc的不等式，由它可以判断算法是否无条件稳定，若不是，将给出稳定条件。例51分析Newmak方法、Wilson-6方法的稳定性解：将(5-38)代入(5-43)有2二.f2(-1)_022.(-(1-2),J_1<02显然，当(5-45)算法无条件稳

11、定。1()-J.-'2(-)2222_Jc=2-223一(5-46)算法稳定，但为条件稳定，其中Qc为临界采样频率。由于（5-43）式仅仅是充分条件，所以可进一步按照稳定性的定义得到5.1.2节叙述的无条件稳定条件。对Wilson-6方法,将(5-41)代入(5-44)得6c2_06,12(21),0<12+Q2(_1+602-60)>0_3_2_2(40+1-60)Q+24c-12>0建4(2c2+1-3c)占0(5-47)容易看出，其中第一，二，五不等式恒成立，对第三，四不等式若希望对任意的则有：Q均成立,求解上述不等式得实际使用中通常选取1=1.45.2.3算法

12、的相容性和收敛性直接积分算法的相容性、2_2_-1+6g-66>04日36日2+1>01.31一1.372(5-48)收敛性分析同样要使用其位移型的差分方程,或对应的单步多值形式。在算法（5-30）式中，用精确解代替近似解，即可得到局部截断误差表达式，用符号e（tk）表不y(tki)=Ay3),Lk.he3)(5-49)(5-50)局部截断误差表达式用放大矩阵的特征量以最常用的线性三步法为例可表示为2e(tk)=以超h)2Alx(tk)Azx(tkh)人3乂询2h)/hx(tkh),x(tk-2h)其中A1,A2,A3分别为应的3X3的放大矩阵的三个特征向量，然后将在tk点进行泰勒

13、展开，然后利用运动平衡方程化简即可。若局部截断误差表达式为步长的s(s>0)阶小量，则称算法是s阶相容的。对2M2的放大矩阵，可仿照上述步骤，来验证算法的相容性。在经典的数值算法收敛性分析理论中，一个重要的结论就是相容加稳定等于收敛，其相容的阶数就是算法的精度阶。收敛性的含义也是当时间步长趋于零，算法的数值解趋于精确解。对直接积分算法该定理同样可以证明是成立的。例52分析Newmak法的相容性和精度解：其局部误差仿照(550)式得:e(tk)x(tkh)2Ax(tk)-A23h)h2(5-51)也可以由(5-42)是直接求得，即将x(tk+h),x(tk_h)在tk点泰勒展开，并注意到在

14、tk时刻的运动方程有:121.1(3)e(tk)=(¥)夏+(/6+¥)£Q_(¥)x()h262(5-52)12-1.112(4)23-(_、：,一)-1.(.：.-)1'.1x()h-o(h)261212241显然，当物理阻尼为零时，选择¥=1算法是二阶的，即截断误差是步长的二阶小量。物理21阻尼的存在，使算法精度降了一阶，但若同时选择a=1,算法精度仍然是二阶的，一般称6为Newmak线加速度法。显然Newmak方法中有两个参数待定，每种特定的选取都是一个特定的算法，最常用的几个算法见表5-1表5-1:常用的Newmak族直接积分

15、算法16方法名称稳定条件无阻尼问题精度阶有阻尼问题精度阶1/21/4平均加速度方法(梯形法)无条件211/21/6线性加速度方法'c=23=3.46221/201ifc=221显式如果在一个时间步内需要求解一个隐式的方程组，则称算法是隐式的，反之不需要求解方程，直接计算即可得到下一时刻的值，则称算法是显示的。从5.1节的Newmak方法的计算步骤可以看出，这类方法是隐式的，但对于中心差分方法，若质量矩阵和阻尼矩阵都是对角矩阵就可以显示地计算。显然显示方法计算量要小得多。读者可自行分析Wilson-日方法的精度，不难分析，无论是无阻尼还是有阻尼其精度都是2阶的，它也是隐式方法。5.2.4

16、算法耗散和弥散特性算法的精度，在小步长的情况下可以通过局部截断误差分析来说明比较，但是，在实际计算过程中，步长的选取可能不是很小，此时如何来度量算法的计算精度，当然可以针对有解析解的问题进行大量的数值计算，将数值解与解析解进行比较来分析算法的计算精度。理论上还可以通过数值耗散(disspation)和弥散(dispersion)来辅助度量与分析，为引出这两个概念的含义，我们仍然以单自由度有阻尼自由振动问题为例，该问题的解析解为：x(t);e一,t(c1cos-dtc2sin，dt)(5-53)当直接积分算法用于这样的问题，前小节已经讲述过它可以写成形如(5-42)的关于位移的有限差分形式，就是

17、可以得到一个关于位移的有限差分方程，对于一个收敛的且有一定精度的算法，这个差分方程通常有一对共扼复根，1,2=e(<d)h(5-54)其中孰际1金,该两根称为主根，其它根称为寄生根(spuriousroots)。解的一般形式可写为(5-55)-一tn,.一,.=e(c1cos-dtn-c2sin.dt式中的£称为算法阻尼比，当有物理阻尼存在时，它还包括了物理阻尼的影响，5称为算(5-53)式形式是相同的,法频率，对应的T=2冗/缶称为算法周期。可以看到上式前两项与这给了我们与精确解进行比较的可能，如果1)寄生根的影响较小，即i=3'm：'1.22)解表达式中得常

18、数Ci，C2与Ci，C2差别不太大。这样，我们就可以通过比较不同算法的算法阻尼比和相对周期误差，来比较算法的计算精度。此时，对无阻尼问题，可以很明显看到算法阻尼比会使得数值解曲线的幅值与解析解相比要降低而产生振幅衰减，这就是所谓的算法的数值耗散。同时，不同的算法的算法周期与精确的周期会有一定的误差，这个误差一般用相对周期误差来表示(T-T)/T,它会使得数值解曲线上产生周期的延长或缩短，即所谓的数值弥散。实际分析时，可以首先通过求解放大矩阵的特征方程得到特征方程的主根和寄生根，若主根可以表示为：1,2(5-56)并注意到式(5-54)有d-arctan(b/a)(5-57),-arctan(b

19、/a)/.1.(2(5-58)(5-59)-199_二Ln(a2b2)(5-60)2-对结构动力学问题，一般总希望算法在低频段有较小的耗散和弥散，而且在低频段，当Ct。时，可以很方便地获得它们的近似解析表达式。在高频段算法的耗散特性，用谱半径来说明更适合，一般用'=m9(A)来度量算法对高频分量的耗散特性，特别地，当p瓷=0时，称算法具有高频渐进消去特性，即当算法计算一步以后高频极限完全地被耗散掉，而其它高频分量由高到低渐进地被耗散。由前面的叙述可以看到，算法放大矩阵的特征向量A1,A2或A1,A2,A3决定了算法对应特征方程的根，也就决定了算法的稳定性，同时确定了谱半径、以及算法耗散

20、和弥散特性。这些特性有时也称算法的谱特性。放大矩阵相同的不同算法，称为互相相同的，放大矩阵不同但特征值相同，称算法互相相似，或称算法频谱等价的。对于算法的耗散特性，应该说明的是高频耗散特性对实际的结构动响应求解是有益的，因为实际结构进行有限元离散以后计算出的高频行为并不真正代表系统的物理行为，它是结构系统在空间进行有限元离散的结果，是虚假的行为，而不具备高频耗散特性的算法，是将所有频率上的响应全部进行了积分，尽管步长取得相对较大时，高频的积分不准确，这样的计算结果显然与系统实际的反应不一致。但是同时要注意，它在低频也不同程度地引入了数因为长时间以后应该精确计算的低频响因此，有耗散特性的直接积分

21、方法只适算法才可能有二阶精度，我们仅讨论这一值耗散，这样这些算法就不适合进行长时间的计算，应，由于耗散特性的存在，已经被耗散得面目全非,合计算瞬态的、短时间内的低频动力响应。例5-3分析Newmak族算法频谱特性1解：对Newmak族算法来说，当¥=时,2种情况。此时算法放大矩阵的两个特征量为:Ai1丫：二；2(5-61)1,2aiib则谱半径P=_a2,b2=A2(5-62)不考虑物理阻尼时，对任意的有P=1,=0,(5-63)(5-64)T-T-111111T-1212arctanc、1(.：.i/(1-(；.一*i2)就是无阻尼时，Newmak族算法不存在数值耗散，但有一定的相

22、对周期误差。对有阻尼问题，P8=1,£=七十0(C2)5.2.4算法的超调特性谱半径这个指标对算法性能的影响还需要进一步说明的是，它只决定算法的长期特性,即pW1可以保证随着算法计算步数的增加计算过程是数值稳定的。但对无条件稳定的算法，由于步长大小选择没有限制，一般在满足指定精度的条件下，尽可能取较大的时间步长，对于非零初始条件问题，在计算开始的几步可能会出现初始数据及其误差(如初始位移，速度的测量误差，初始加速度的计算误差)被放大的现象，这称为超调(overshoot)。这种现象是放大矩阵A病态，有较大的条件数而产生的。实际应用时，由于当Ct0算法是收敛的，不会出现超调。一般为简单

23、起见，只分析当Ct必时，在计算的第一步是否会出现超调。例54分析Newmak平均加速度法的超调特性为分析简便起见，将Newmak平均加速度法用于无阻尼自由振动问题，此时其放大矩阵为：A二一4+。24EC+C3e/24=2vi4-2X0V042在Ct笛时，可得近似等式x1=o(1)xO,Vi=o(Q)x0-O(1)V0(5-65)其中o(1),o(C)分别表示关于的零次和一次关系式。由此可知算法在位移上无超调，但由于初位移的影响，在速度上有关于。线性超调现象。前面提过Wilson-日方法有很强的超调现象，对无阻尼问题，从放大矩阵的各元素的表达式中，很容易得到xi(二-u2/14.(2V一6)：】

24、212v1122u.i2(31一2丁)-I.1312，122-2'-?2xoh(V2-1)'J2-6)x0'77-v06”J2(U2-3)-J2-6iV067小2(5-66)(5-67)在Ct妙时，可得近似等式Xi=o(LI2)X0-o(h)v0Vi=o(i.1)V0-o(1)X0其中o（c2）,0（夏）分别表示关于的二次和一次关系式。由此可知Wilson-9方法在位移上关于初位移有二次超调，同时关于初始速度有一次超调。在速度上有关于初位移一次超调。另外，也可以用数值计算的方法对指定的初始条件，计算出近似解x1,v1,然后与精确解比较，或计算系统能量范数:En22Mvn

25、-Kxn(5-68)然后将E0VE1比较。显然，由于直接积分方法适合于短时间的瞬态问题计算，因此超调现象也是必须加以注意的。综上所述，对一个数值积分算法理论上要分析其相容性，稳定性，数值的耗散与弥散特性，对无条件稳定的算法还要分析其超调特性。这样才可能对算法的本质有深入的了解，进而指导数值计算结果的解释与分析。此外，由于直接积分方法对结构运动平衡方程进行数值积分的目的在于估计结构真实的动力响应。为了精确地预计结构的动力响应，要求模态分解以后所有的单自由度平衡方程都必须被精确地积分，但在直接积分法中，对所有的方程积分都相当于采用相同的步长h,所以时间步长的选择必须针对系统的最小周期。如果Tn是系

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

27、频响应就可以了。也就是步长可选择为Tp/10,比Tn/10大Tp/Tn倍。由于实际情况中Tp/Tn可能会非常大，这样条件稳定算法的稳定条件就可能无法满足，而无条件稳定算法对步长的选取就没有稳定性的限制，因此对于实际的结构动力响应计算，多数都使用无条件稳定算法。5.3矩阵特征值问题及解法5.3.1 .问题分类（K-2M）=0K='M'这原本是广义特征值问题，但可以化为标准特征值问题，前乘M,，得12.（M*-，I）=i令A=M-K,则有（A一，I）=这是代数里的标准的矩阵特征值问题，对结构动力学问题而言M对角tM,K=A对称M对称非对角tM,K=A不一定对称非对称矩阵的结构特征值

28、问题计算量很大，此时作如下处理：1对M进彳tCholesky分解（因为M对称正定）M=LUL（下三角阵）U=LT（上三角阵）代入方程（K-2M）=（K-2LU）前乘L工得（LlK-.>:：2U）二0提出U得（L工KU/；.?2）U=0令A=L°KU",S=LT4则有（A-2I产=0可以征明A与a有相同的特征值，因为它们是相似的矩阵，但A是对称的T因为A=（L,KL）T二l_kTL二二LKl二二A所以A对称。只不过特征向量有所变化，求出=L至于Cholesky分解，可直接由1I11I1211*22。L=n11n21nn的各非零元素j_L、,一，2、1/21jj=jj乙1

29、jr)r.Xj-1ij=(mij-%1ir1jr)/1jjrA上式依次取j=1,2,n,艮%111v12v220L二Jn1Vn2Vnn一则其下三角之i=1,2,.,n行元素,由以后(或K)的元素计算得到i=j+1,j+2,n求得L的下三角部分各列元素。又若记依次为（4）Am阵的特征值是,特征向量仍为x,m为整数；（5）若A非奇异，则A的逆矩阵存在，记为A-,是A的特征值为1/九，特征向量为x;（6）若矩阵B与A相似，即有可逆阵P存在，使_1一B=P-AP则B的特征值也是九，特征向量是P-x4.特征值的积与积若矩阵A的特征值为兀、九2、.，鼠，则有n1''2.：''

30、;''n二、aiii4.i"2.Rn=det（A）它们可作为校核、估计甚至计算特征值的手段，读者可以自行验证。5.3.3特征值问题的基本计算方法由于高于四次的一元代权方程无法求精确解，所以必然要采用迭代方法求计算求近似解。目前的迭代方法，只不过迭代的方法和技巧不同。已有很多成熟方法。方法选择时主要取决于1）所要求的特征对数目2）K、M的阶，3）K、M的带宽以及是否带状。根据用到的基本关系大致可分四类.向量迭代法（哥迭代法Poweriterationmethod）基本关系Ki-Mi根据迭代格式不同又分：正向迭代（哥迭代）逆迭代（反哥迭代）主要用于求特征向量。.变换法基本

31、关系：1TK:：J-.:M中-I中质量归化为的模态矩阵A=diag（',与，3,，'n）雅可比迭代（Jacobi）（小模型实对称阵标准特征值问题的全部特征对方法）广义雅可比迭代（求全部特征对，小模型的广义特征值问题，或者有大量非对角线零元素和少量对角线零元素问题）豪斯霍尔德（Householder）.多项式迭代法（不单独使用）基本关系：；()=det(K-M)显式多项式迭代隐式多项式迭代.斯图姆（Sturm）序列法利用P（九）=det（KUM）的Sturm序列性质来求解。实际工程结构的动力学问题绝大多数使用有限元法方法离散求解。随着计算机存储和速度能力的成倍提高，所建模型越来越

32、精确，几十万、上百万个自由度的计算成为可能。为此发展了一些针对大型特征值问题的方法，它们综合上述典型方法和技巧，常用的有：多项式迭代行列式探索法（带宽较窄的低阶特征值问题。子空间迭代法（Subspaceiteration）（与行列式法相似，但计算量小）兰索斯（Lanczos）法向量迭代HQRI（HouseholderQR迭代方法）（针对大型满阵的大多数或全部的特征值和向量）逆迭代方法任意给定初始的第一阶振型向量和第一阶特征值，一般常用的选取为x1=1,1,.,1=1利用关系式Kx=.Mx按下列方式实施迭代1）选Xi=1,1,1T,%二12）令K=1,2,.i）KXk1=MxkK正定利用K的乔莱斯基分解，进行线性代数代权方程组的求解，解出£中ii）对xk中实施对M的正定化，得新的下一次迭代向量Xk1xk1'Zr2（xk1MXk1）工1=,/x：1KXk1然后返回i）使ktk+1进行下轮迭代，直至相邻的两次迭代值的差6%中4<6=10-6normxk+-xk<w=10k1k雅可比（Jacbi）法雅可比（Ja