結(jié)構(gòu)動(dòng)力學(xué)方程常用數(shù)值解法

1、結(jié)構(gòu)動(dòng)力學(xué)方程常用數(shù)值解法對(duì)于一個(gè)實(shí)際結(jié)構(gòu)，由有限元法離散化處理后，動(dòng)力學(xué)方程可寫為:M x C x Kx = F （t）從數(shù)學(xué)角度看，這是一個(gè)常系數(shù)的二階線性常微分方程組，計(jì)算數(shù)學(xué)領(lǐng)域，常微分?jǐn)?shù) 值算法常用的有兩大類：-、針對(duì)一階微分方程數(shù)值積分法發(fā)展的歐拉法，中點(diǎn)法， Rugge-kutta （龍格一庫(kù)塔）方法。二、直接基于二階動(dòng)力學(xué)方程發(fā)展的方法。對(duì)結(jié)構(gòu)動(dòng)力學(xué)問(wèn)題的數(shù)值求解，常用的有兩大類：一是坐標(biāo)變換法，它是對(duì)結(jié)構(gòu)動(dòng)力 方程式，在求解之前，進(jìn)行模態(tài)坐標(biāo)變換，實(shí)際上就是一種Rize變換，即把原物理空間的動(dòng)力方程變換到模態(tài)空間中去求解?，F(xiàn)在，普遍使用的方法是模態(tài)（振型）迭加法。 二是直接積

2、分法，它是對(duì)結(jié)構(gòu)動(dòng)力方程式在求解之前不進(jìn)行坐標(biāo)變換，直接進(jìn)行數(shù)值積 分計(jì)算。這種方法的特點(diǎn)是對(duì)時(shí)域進(jìn)行離散，然后將該時(shí)刻的加速度和速度用相鄰時(shí)刻 的各位移線性組合而成。通常又稱為逐步積分法。模態(tài)迭加方法，比較常用，但如下情況通常使用直接積分方法（即求解之前不進(jìn)行模 態(tài)分析）一、非比例阻尼，非線性情況。二、有沖擊作用，激起高頻模態(tài)，力作用持續(xù) 時(shí)間較短，模態(tài)迭加計(jì)算量太大。振型迭加法與 Duhamel積分?jǐn)?shù)值解按照有限單元法的一般規(guī)則，經(jīng)過(guò)邊界條件的約束處理，結(jié)構(gòu)在強(qiáng)迫振動(dòng)時(shí)多自由度 體系的運(yùn)動(dòng)平衡方程可以表示為：MU CU KU =R其中，M是體系的質(zhì)量矩陣，C是體系的阻尼矩陣，而K則是剛度矩

3、陣.R為外荷載向 量.U、U和U則分別是體系單元節(jié)點(diǎn)的位移、速度和加速度向量 .上述動(dòng)力平衡方程 實(shí)質(zhì)上是與加速度有關(guān)的慣性力MU和與速度有關(guān)的阻尼力CU及與位移有關(guān)的彈性力KU在時(shí)刻t與荷載的靜力平衡。振型疊加法是把多自由度體系的結(jié)構(gòu)的整體振動(dòng)分解為與振型次數(shù)相對(duì)應(yīng)的單自由度 體系，求得各個(gè)單自由度體系的動(dòng)力響應(yīng)后，再進(jìn)行疊加得出結(jié)構(gòu)整體響應(yīng)振型疊加 法原理是利用結(jié)構(gòu)無(wú)阻尼自由振動(dòng)的振型矩陣作為變換矩陣，將結(jié)構(gòu)動(dòng)力方程式（1）式 變換成一組非耦合的微分方程逐個(gè)地求解這些方程后，將解疊加即可得到動(dòng)力方程的 解。將體系單元節(jié)點(diǎn)的位移向量表示為如下的變換形式：U (t) _X(t)式中的變換矩陣是

4、由動(dòng)力方程對(duì)應(yīng)的無(wú)阻尼自由振動(dòng)方程解出的前m階振型矩陣.即：J = 1 2. m； X（t）是與時(shí)間有關(guān)的m階向量,X的各分量稱為廣義位移。將式（2）代入動(dòng)力方程（1）并左乘以叮T，則可得廣義位移為未知數(shù)的方程: Mx（t）CX （t） Kx（t）二 R（t）（3）式中M =otm，C=oTc，K=otk，R=otr（4）現(xiàn)在進(jìn)一步考察式（4）.考慮到特征向量的正交性，可得TM=1tK6=A（5）于是對(duì)應(yīng)于振型的廣義位移的平衡方程（3）可改寫為X（t）亠處TC：X（t）上X（t） = ： JR（t）（6）其中，上為特征值（7）將式（2）稍加運(yùn)算可得廣義位移用有限元位移表示的形式x - ：tmu

5、（8）在（6）式中，當(dāng)忽略了阻尼的影響，平衡方程為互不耦合的，可以對(duì)每個(gè)方程逐個(gè)地進(jìn) 行時(shí)間積分出于相同的考慮，在對(duì)有阻尼的體系進(jìn)行分析時(shí)仍然希望采用相同的計(jì)算 過(guò)程去求解互不耦合的平衡方程式問(wèn)題是式（6）中的阻尼陣C通常不能象體系的質(zhì)量 陣和剛度陣那樣由單元的剛度陣和質(zhì)量陣裝配而成但當(dāng)假定阻尼與固有頻率成比例，即假定'Tc = 2jj(9) 式中，是振型阻尼參數(shù)；v是Kronecker符號(hào)（當(dāng)i = j時(shí)，'ij = 1.當(dāng)i 一 j時(shí)，j = 0）這時(shí)式(6)可簡(jiǎn)化為如下形式的若干個(gè)方程式X(t) 2 打 jX(t)2川)訂(10)其中x(t)的初始條件為下式X t= &#

6、176;i MU 0 Xi t=l MU 0(11)2 ”式(10)表示了一個(gè)具有單位質(zhì)量，剛度為 '的自由度體系當(dāng)阻尼比為i時(shí)的運(yùn)動(dòng)平衡 控制方程。這個(gè)平衡方程的求解可通過(guò)計(jì)算 Duhame積分求得。X (t)=丄 fori (t)ec°)si(t)di +e_°fi sin旨£ + Pi cost)(式中(13)同時(shí)在計(jì)算上也避免計(jì)算阻尼陣而只需計(jì)當(dāng)利用式(9)來(lái)考慮阻尼的影響時(shí)意味著假設(shè)結(jié)構(gòu)的總阻尼是每個(gè)振型的阻尼之和, 而每個(gè)振型上的阻尼是能夠量測(cè)的，況且在大多數(shù)情況下結(jié)構(gòu)的阻尼比更易于量測(cè)。因 而便于用來(lái)近似地反映結(jié)構(gòu)體系的阻尼特性。 算剛度陣和

7、質(zhì)量陣。積分遞推公式對(duì)以上方程式(10)，考慮某一模態(tài)的振動(dòng),并略去下標(biāo)i可寫為x (t) 2 x(t) 2X訂(t)(14)在初始條件t4X0X心0 =X0(15)下的定解為x(t) = eff ttsin (t-t0)+co3 (t-匚)人(t-t0)*x)+丄(t-5«cot0(16)式中J- 2 ,將上式對(duì)時(shí)間求導(dǎo)，得I%21)sin t - to X。 e-；t0x(t) = b z- % .cos t - to - =sin t - to x°丄r( 亠° t0- sint - cos_ t- d(17)cox(t 2 ) =-22z Isin 2也

8、x(t)將上式中的to , t分別代以t , t V (其中2厶為時(shí)間步長(zhǎng))，并按拋物線法則計(jì)算式中的 卷積，有下式：_2蘆八人COx(t 2 ) = e=sin 2cos2 p x(t)co丄 e，sin 2 _* x(t)e_2sin 2 _ r34._e_sin : r(t t)(18)3e cos2- =sin2： x(t)©e(cos2_ - sin 2 p r(t)3-sin：)r(t ) r(t 2 )必3(19)以上兩式以矩陣表示為:21ai2a22ai3a23ai4a24a15a25x(t)x(t)r(t)r(t+A)j(t+2®(20)記*e =e q=

9、q1 2 co 二盲c, =cosAq =si n 矗 s,? ? ?二 2 § qC2 二 2c； -1（21）則上述A矩陣元素為a1 二 e (c?C0S2)耳2 = e2S2 / ；-*13=2*25a14 = 4e1 §a25 / a21 二-a12*22 二勺(c2 - 4 勺)*23 = *22*25 *24 = 4ei (c1 - $ § )*25a25 =厶 / 3當(dāng)a求出后，以應(yīng)用上述遞推公式，以前一時(shí)刻來(lái)求后一時(shí)刻的結(jié)果。計(jì)算不重復(fù) 后在時(shí)域中的步進(jìn)求解只是一些簡(jiǎn)單的數(shù)組相乘。計(jì)算速度很快。Newmark 類方法F面討論直接積分法1 1&quo

10、t;IQ « «IO « « «AX(tn 1) = X(tn)t X(tn)t X(J )t X(tn)0( t )(1)2 6X(tn 1)= X(tn)t X(tnX(J)0( t3)(2)X (tn 1)=X(tn )X(tn(3)將（3）代入（1）, （2）得:(tn l) OC t4)也 t 用變權(quán)來(lái)調(diào)節(jié) X(tn 1) = X(tn)t X(tn) X二 t .3三 X(tn 1)0 ( t3)二 tX（tn 1 ） = X（tn ）乙 X（tn ）-1性U (廠卄x,"2然后假設(shè)在J 1時(shí)刻近似滿足運(yùn)動(dòng)方程M X n

11、1 CX n 1 KX通過(guò)變換將速度和加速度用位移表示， 得代入運(yùn)動(dòng)方程，只剩n+1時(shí)刻位移一個(gè)未知數(shù),參數(shù)不同選取包含著三個(gè)經(jīng)典算法(1)Newmark平均加速度法2Newmark線加速度法(3)中心差分法Newmark法的一般步驟:1初始值計(jì)算(1)形成系統(tǒng)矩陣K, M和C定初始值X。，X。和X。(3)選擇時(shí)間步長(zhǎng):t，參數(shù)、并計(jì)算積分常數(shù)：a0 = a1a廠廠一1yAt 丫a4 二 -1a52(丁一2)比=小(1一？)a7=t形成等效剛度矩陣KK a o M a 1 C(5)K矩陣進(jìn)行三角分解二 LDL T 2對(duì)每一時(shí)間步(1) 計(jì)算t 厶t時(shí)刻的等效載荷Q t：t= Q t M (a6

12、 Xta2xt a3x t)C (axta4xta5 x t)(2) 求解：t時(shí)刻的位移Xt CtXt a7（Xt tXt）(3) 計(jì)算t rt時(shí)刻的加速度和速度xt：t= a°(Xt t 一 x一 a? xt - a3 xtx t t 二 x t a 6 x t a 7 x tWilson- r法的一般步驟:1初始值計(jì)算(1)形成系統(tǒng)矩陣K, M和C定初始值x。，x oxo 0(3)選擇時(shí)間步長(zhǎng),并計(jì)算積分常數(shù)v -1.4ao(切)2aiF:ta2a3 - 2a4a。a23a6 =1a7 二 ta8 =：t2形成等效剛度(5)將等效剛度進(jìn)行三角分解二 LDL2對(duì)每一個(gè)時(shí)間步長(zhǎng)(1)

13、計(jì)算t *時(shí)刻的等效載荷Rt 珂=Q (Q t -Q) M(a0x a2xt 2xt) C(x 2xt a3xt)(2)求解t * nVt時(shí)刻的位移(LDL T)Xtr：t(3)計(jì)算在t負(fù)時(shí)刻的加速度、速度和位移xt t 二 a4(xt r t 一 xt) a5 xt a6 xt為 t 二 Xttxt a8（Xt tXt）三結(jié)構(gòu)動(dòng)力響應(yīng)數(shù)值算法性能分析算法數(shù)值計(jì)算結(jié)果如何評(píng)價(jià)，針對(duì)不同的結(jié)構(gòu)動(dòng)力響應(yīng)計(jì)算問(wèn)題應(yīng)該如何選擇更合適 的算法等是非常重要的問(wèn)題。這就需要深入研究算法的數(shù)值計(jì)算性能，也就是算法的計(jì) 算精度、穩(wěn)定性等。對(duì)線性結(jié)構(gòu)動(dòng)力學(xué)問(wèn)題，已經(jīng)有證明對(duì)整個(gè)多自由度的積分，等價(jià)于將模態(tài)分解后對(duì)

14、 單自由度的積分的結(jié)果進(jìn)行模態(tài)疊加，因此可以通過(guò)對(duì)單自由度問(wèn)題的分析，來(lái)說(shuō)明算 法的特性，其中阻尼均假設(shè)為比例阻尼。算法用于結(jié)構(gòu)動(dòng)力學(xué)方程的有限差分表示為：x 2x2x= f(t)以下算法的性能分析，均將算法用于這個(gè)方程。分別在相鄰的不同時(shí)刻應(yīng)用算法可得如 下一般形式：y k 1 二 Ay k L k其中A為放大矩陣或稱逼近算子,Lk為載荷逼近算子。2yk - X k , X k -1 ,yk1 一 X k 1 , X k ,對(duì)自由振動(dòng)情況有y Any°顯然計(jì)算的第n步的值與A直接有關(guān)。例如，Newma方法:-1A 二 AtAd：1 +h2 甩 22h203At =21 +2h

15、9;o _矩陣A的特征多項(xiàng)式為h2|1_（1_2®co2h1_h（1_2P）31 2h（1 -）'-蛍h（1 丫）皎-det（ A - I H ' 22 A< A2 = 0A1traceA2（A11A22 ）A? = det A = A11A22 - A12 A21 對(duì)Newma方法有：1 (2 一1)(A、1 21(2-2)門()門A2其中為時(shí)間步長(zhǎng),Wilson-二方法,放大矩陣為:-31'-1 22丁(門 2f 23)6)2'-j/2)6" 11-3二 2,；16)6h 171h2h(6，D 二('1放大矩陣A的特征多項(xiàng)式

16、為:-det(A - 1) = 3 - 2A)2 A? 1 - A3 - 0其中A1, A2, A3為該矩陣的三個(gè)特征向量，分別為矩陣的跡的一半、各階主子式的和以 及矩陣的行列式，分別表達(dá)如下：A318 二-63 -2 一 3, , 2二 2 _ 31 r2日(0 2日 2 + 6)4 門 23門 2二 3 一 62r 218- 12丁(門 2r6)_66，2_ 32寸 2.門 2丁 33： 12二珥門2寸26)此外，在幾個(gè)不同時(shí)刻應(yīng)用數(shù)值算法，然后將方程中的速度和加速度項(xiàng)消去，可得數(shù)值 算法關(guān)于位移的差分方程，例如 Newma方法，有y 11(2-2)門('1 2f + 2)。 Xn

17、_1 =0(12'-門 2)Xn .1 - 21(2 一 1)(八 2Xn3.1算法的穩(wěn)定性分析穩(wěn)定性定義：設(shè)漏iTZ'm為放大矩陣A的特征值，則卩= maxh定義為A的譜半徑,若特征值互異，則'叨的算法是穩(wěn)定的，但若有重特征根，則要求 門。如果算法的穩(wěn) 定性要求對(duì)步長(zhǎng)的選取有限制，稱算法是有條件穩(wěn)定的，反之為無(wú)條件穩(wěn)定的。判斷方法：放大矩陣的譜半徑小于等于 1成立的充分條件是1-2A +A >0*1+2A+A2 蘭01 A 2 X 0對(duì)3 3的放大矩陣1-2A 1+2A2 _A 3-03-2A 1A23A3-032A 1A2-3A3-012A 1+A2A3-01

18、A2 -A 3(2 A 1A 3)-0上兩式是關(guān)于算法自由參數(shù)的不等式，由它可以判斷算法是否無(wú)條件穩(wěn)定，若不是,將給出穩(wěn)定條件。算法數(shù)值穩(wěn)定性的物理解釋：物理上，對(duì)一個(gè)無(wú)阻尼或者有阻尼自由振動(dòng)系統(tǒng)，系統(tǒng)的能量隨著時(shí)間不應(yīng)該增加，有 阻尼情況還應(yīng)該減小。因此，一個(gè)數(shù)值方法的計(jì)算結(jié)果也不應(yīng)該放大初始能量，如果經(jīng)過(guò) 若干步的數(shù)值計(jì)算以后，計(jì)算結(jié)果遠(yuǎn)比初始條件大，那就是數(shù)值算法本身計(jì)算是不穩(wěn)定的。分析Newma方法、Wilson-二方法的穩(wěn)定性將Newma方法放大矩陣特征量代入穩(wěn)定性分析表達(dá)式門2(-)-022(2)(_ 2 )-0顯然，當(dāng)> 1 ,6 > 722，算法無(wú)條件穩(wěn)定。對(duì) Wi

19、lson-2602 H0602(2日-1P012 +02(-1 +6日2 -6日)工0(4日'+1 -6日2)02 +24日-12 啟004(282 +1-30)30容易看出，其中第一，二，五不等式恒成立，對(duì)第三，四不等式若希望對(duì)任意的均成立，則有：工1 6乎-6： _04二3 -6二21 _0求解上述不等式得1 + V3丄 3 1.372實(shí)際使用中通常選取，=1.4。3.2算法的精度(相容性和收斂性)直接積分算法的相容性、收斂性分析同樣要使用其位移型的差分方程，用精確解代替 近似解，得到局部截?cái)嗾`差表達(dá)式，用符號(hào)e(tk)表示y(tkJ =Ay(tQ Lk he(tQ以最常用的線性三

20、步法為例，局部截?cái)嗾`差用放大矩陣的特征量可表示為2e(tQ 珂x(th) 2Ax(tk) A2X(tk -h) 一 Ax(tk -2h)/ h然后將x(tk -h)，x(tk -2h)在tk點(diǎn)進(jìn)行泰勒展開(kāi)，利用運(yùn)動(dòng)平衡方程化簡(jiǎn)即可。算法精度定義：若局部截?cái)嗾`差表達(dá)式為步長(zhǎng)的s階小量，則稱算法是s階相容的。相容加穩(wěn)定等于收 斂，其相容的階數(shù)就是算法的精度階。收斂性的含義是當(dāng)時(shí)間步長(zhǎng)趨于零，算法的數(shù)值 解趨于精確解。分析Newma法的相容性和精度其局部誤差表達(dá)式得：“、x(tk+h)-2Ax(tJ+A2(tkh)e(tQ 二2h即將x(tk h),x(tk -h)在tk時(shí)刻點(diǎn)泰勒展開(kāi)，并注意到在該時(shí)刻的運(yùn)動(dòng)方程有：e(tk)之 )i(-2、一 -丄)門-(-1)x(3)h (!) 2 1 丄(!rx h2 o(h3)26226121224顯然，當(dāng)物理阻尼為零時(shí)，選擇 2算法是二階的。物理阻尼的存在，使算法精度降了 一階。Newmak方法中有兩個(gè)參數(shù)待定，每種特定的選取都是一個(gè)特定的算法，最常用的 幾個(gè)算法見(jiàn)表常用的Newma族直接積分算法16方法名稱穩(wěn)定條件1/21/4平均加速度方法