1、空間動力計算問題的有限元程序編制    摘 要：本文利用大型有限元分析軟件ANSYS 進行前處理，編制了適用于空間動力計算問題的子空間迭代程序，隨后應用AutoCAD 內嵌的VBA 語言和Fortran 語言混合編程的方法，實現了結構振型圖的可視化，然后通過懸臂梁數值考題，與ANSYS 計算結果進行對比，驗證了程序求解結構自振特性的正確性；本文介紹了線性加速度法、Willson- 法和Newmark法這三種常用的逐步積分方法各自的特點，推導了振型分解時程分析法的相關公式，隨后編制了能同時實現這三種逐步積分方法的計算程序，并通過數值考題驗證了程序的準確性。

2、結果表明自編動力計算程序能得到較好的精度，可以有效地應用于空間動力計算問題。關鍵詞：動力計算；子空間迭代法；振型分解時程分析法本文利用ANSYS 進行前處理建模并導出模型數據文件，然后用已有程序計算并輸出結構集中質量矩陣、剛度矩陣、結點自由度指示矩陣等所需數據。在此基礎上，用Fortran 語言編制了子空間迭代法程序以及振型分解時程分析法程序，并通過AutoCAD 內嵌的VBA和Fortran 語言混合編程的方法實現了結構振型圖的可視化。1 子空間迭代法程序本程序使用子空間迭代法來計算結構的自振特性，適合于求解大型結構的前若干階自振特性。本程序計算時，首先讀入前述數據及求解自振特性的控制參數，

3、然后利用子空間迭代法進行求解，最后輸出結構的自振特性，包括自振頻率、振型、振型參與系數等數據。1.1 程序說明子空間迭代法的具體計算原理以及流程，本文不再贅述。現僅對程序中質量矩陣形式的選擇以及初始向量的選取做一點說明。結構的質量矩陣有一致質量矩陣和集中質量矩陣兩種，本程序采用的是后者。當采用一致質量矩陣時，所求得的頻率偏高，一般是真實頻率的上限值；而采用集中質量矩陣時，所求得的頻率偏低，為真實頻率的下限值。因為集中質量矩陣假定的單元位移模式與真實變形是有差距的，這等于人為地引入了約束，增加了結構的剛度，從而可能一起頻率值的偏低，但由于采用有限元離散結構，將質量集中到若干結點上，人為地去掉了一

4、些約束，等于減少了結構的剛度，兩者結合可抵消一部分對結構的影響。故若恰當地分配質量，用集中質量矩陣也可以得到較為精確的結果。此外，由于一致質量矩陣為滿陣，數值計算費時，而集中質量矩陣為對角陣，占用內存少，計算簡單且省時。因此，在實際動力分析中一般都采用集中質量矩陣。在使用時要注意：若使用高精度單元，網格剖分較粗，將會使計算結果偏差較大，故一般要控制單元的大小，單元越小，計算結果越精確。初始迭代向量組的選取，將直接影響迭代的次數以及結果的精度。通過分析比較1，本程序使用如下方法來形成初始迭代向量組Y0=MX0。首先令Y0中的第一列為M中的全部對角元mii，這就保證了全部有質量的自由度都受到激發而

5、不致遺漏振型。其他各列為單位向量，在比值mii/kii 的最大坐標上置。1.2 懸臂梁考題- 2 -由于水利工程中的結構大都體形龐大，又由于劃分網格時往往受到單元數的限制，故每個單元的尺寸都較大。為了解自編程序在單元尺寸較大時的計算情況，故考題取一3000mm*2000mm*20000mm 的懸臂梁。首先將其劃分為120 個邊長都為1000mm 的立方體單元；然后為了考察單元細分對計算結果的影響，再將其劃分為7680 個單元，每個單元邊長都為250mm。a.粗網格（120 個單元） b.細網格（7680 個單元）圖1. 懸臂梁網格劃分示意圖分別用不同的方法計算兩種劃分情況下此梁的前十階自振頻率

6、，計算結果見下表。其中，用ANSYS 中的子空間迭代法進行計算時，無論此梁采用粗網格還是細網格，在迭代100 次以后，仍無法收斂，且迭代過程中的頻率值與準確值相差甚遠。所以在用ANSYS 計算時，采用了Lancozs 法。一維懸臂梁理論自振工程頻率值的計算公式為：1 2 3 =3.51602 f , =22.0345f, = 61.6972f，其中，212f EIl A= 。代入各參數，求得1 2 3 =3.83320E04, = 2.40222E 03, = 6.72629E 03，分別對應三維自振頻率的第2,5,9 階。將表1 中的計算值與相應的理論值對照，可以看出，當采用集中質量矩陣時，

7、有限元計算所得的頻率值比理論值偏低，為真實頻率的下限值，且當單元數增加時，計算值可以較好的收斂于理論值，這也與前文的分析相一致。表1. 懸臂梁自振工程頻率比較單元數120 個單元 7680 個單元求解 ANSYS Lancozs 法 ANSYS Lancozs 法方法 SOLID45 單元SOLID95 單元自編子空間迭代法 SOLID45 單元SOLID95 單元自編子空間迭代法1 2.54320E-04 2.56240E-04 2.68902E-04 2.54280E-04 2.54390E-04 2.55237E-042 3.77910E-04 3.79160E-04 3.87706E-

8、04 3.77760E-04 3.77830E-04 3.78403E-043 1.52730E-03 1.53830E-03 1.61007E-03 1.52850E-03 1.52920E-03 1.53404E-034 2.02150E-03 2.02180E-03 2.02276E-03 2.16960E-03 2.17000E-03 2.17306E-035 2.16940E-03 2.17650E-03 2.22042E-03 2.17460E-03 2.17460E-03 2.17492E-036 3.96000E-03 3.95980E-03 3.96089E-03 3.958

9、90E-03 3.95880E-03 3.95903E-037 4.03360E-03 4.05660E-03 4.23666E-03 4.03710E-03 4.03840E-03 4.05088E-038 5.45090E-03 5.46680E-03 5.56679E-03 5.45180E-03 5.45280E-03 5.45983E-039 6.05600E-03 6.06870E-03 6.06171E-03 6.52540E-03 6.52640E-03 6.52666E-03各階自振工程頻率(Hz) 10 7.35210E-03 7.38330E-03 7.69092E-03

10、 7.35760E-03 7.35960E-03 7.38150E-03由表1 可以看出，當此梁采用粗網格時，自編程序的計算結果與ANSYS 的偏差較大；當采用細網格時，兩者的計算結果差別不大。此外，當使用ANSYS 中Lancozs 法進行計算- 3 -時，網格的粗細以及單元形式的選取，對結果影響不大；而使用自編子空間迭代程序進行計算時，網格加密后可以使自振頻率值更加準確，此例中，兩者誤差為2%7%。網格加密可以使自振頻率計算值更加準確，但是其計算耗時將成級數倍增長。此例中計算前10 階自振頻率，粗網格耗時半分鐘不到，而細網格則需將近半小時。故有必要找到一個合適的網格尺寸，使計算精度滿足條件

11、，且計算耗時也在可以接受的范圍以內。可見，用自編子空間迭代法程序可以計算出較精確的結構的自振頻率。另外，此例中用ANSYS 的子空間迭代法計算時，結果無法收斂，而自編程序則可以較好的求解，這說明在計算體形較大結構的自振頻率時，自編程序比ANSYS 的子空間迭代法具有一定的優勢。為驗證振型結果的準確性，通過AutoCAD 內嵌的VBA 和Fortran 語言混合編程的方法實現了結構振型圖的可視化。經過對比，自編程序計算得到的振型與ANSYS 的結果完全一致。限于篇幅，下面只給出自編程序算得并繪出的此梁細網格的前三階振型圖。a.第一階振型圖 b. 第二階振型圖 c. 第三階振型圖圖2. 懸臂梁細網

12、格前三階振型圖2 振型分解時程積分法2.1 逐步積分法簡介程序采用彈性時程分析法，就是將結構作為彈性振動系統，輸入地面地震加速度記錄，由初始狀態開始逐步對運動方程積分，直至地震終止，從而求出結構在整個地震過程中的地震反應時程曲線。逐步積分法的特點是，給定初始時刻的位移、速度和加速度后，可求得t1 時刻的位移、速度和加速度，而后逐步求得t2，t3，tn 時刻的解。故只需推導出由t 時刻的解，來求解下一時刻的位移、速度和加速度的計算公式，就建立了逐步積分的算法。一般常見的逐步積分方法有線性加速度法，Willson- 法和Newmark 法。線性加速度法假設質點的加速度在任一微小時段t 內，呈線性變

13、化。此算法是有條件穩定的，要求時間步長t 不能太大，當t 過大時，系統的反應會出現振蕩現象2。同時，為了保證計算結果的精度，一般要求t 小于系統最小自振周期的十分之一。然而，結構的動力響應主要取決于長周期（低頻）分量的影響，步長過短必然會大大增加計算工作量3，因此引入無條件穩定的Willson- 法。Willson- 法是線性加速度法的推廣，它假設加速度在t，t+t內仍是線性的。當=1時，它退化為線性加速度法；當1.37 時，此算法無條件穩定；當=2 時，就是著名的雙步長法3。有文獻表明， 的最佳值為1.420815，但通常取為1.44。此算法會使振動周期延長、振幅縮減，這種效應是由于此算法本

14、身的誤差所導致的，可以看成是一種附加到實際阻尼中的“算法阻尼”，它可以有效的消除不正確的高階模態響應，這也正是Willson- 法的優- 4 -點5。Newmark 法對線性加速度法的假定作了修正，引入了兩個參數r、。當r=0.5 時，只保留，則為Newmark- 法， 值一般取00.25；當r=0.5，=1/6 時，則為線性加速度法；當r=0.5，=0.25 時，相當于加速度在t 內為常量，其值為t 兩端加速度的平均值，即為平均加速度法，為無條件穩定算法，且可以給出滿意的精度6。2.2 振型分解時程積分法由于在計算水工結構的動力問題時，規模一般都十分龐大，如直接采用逐步積分法，每一時間步都需

15、要求解一個大型線性方程集，這將耗費巨大的機時，故本程序采用振型分解時程分析法。它用振型分解的思想，先將運動方程解耦為n 個相互獨立的單自由度方程，然后分別對每一個方程進行逐步積分并求出其響應，最后將各方程的響應進行線性疊加，以求出        結構總的地震響應。振型分解的計算公式為：* ,* 2 *,* 2 *, * ( ) T T T Tj j j j j j j j j j j j j j j j m= M k= K =mc= C = mP= pt以上公式依次計算廣義質量、廣義剛度、廣義阻尼和廣義力。為了減少

16、誤差，程序中計算廣義剛度及廣義阻尼時不采用簡便公式，而分別采用如下公式：* ,* * * Tj j j j j j k= K c=am+bk。其中，阻尼采用瑞利比例阻尼假設1，a、b 是其計算參數，( ) ( ) 1 2 1 2 2 1 2 2 1 12 2 2 22 1 2 12 2a ,b = = 。這樣只需假設結構前兩階的阻尼系數，可以避免假設高階振型的阻尼系數而將帶來的計算誤差。鑒于線性加速度法、Willson- 法和Newmark 法這三種逐步積分方法在計算時只是計算參數取值不同，故本文將這三種方法集成在一個程序中編寫。下面的公式中出現的三個參數、r 和，在選取不同的積分方法時，有不

17、同的取值。具體取值見程序框圖(圖3)。從t 到t+t 時段內的速度和位移增量為：y(t)=y(t)t+y(t)t i ii ii (1)( ) ( ) 1 ( ) ( )2 ( ) ( )22y t=y tt+ y tt +y t t i ii ii (2)式中，在“”上面加一道橫線表示與延長時間步長相對應的增量。從上兩式中解出用y(t)表示 y(t) i和 y(t) ii表達式( ) ( ) ( ) ( 1) ( )2 y t y t y t ty tt = i i ii (3)2( ) 1 ( ) 1 ( ) 1 ( )( ) 2 y t y t y t y t t t = ii i ii

18、 (4)再代入用增量形式表示的動力平衡方程：My(t)+Cy(t)+Ky(t)=P(t) ii i (5)可以得到- 5 -K式中21( ) t tK K M C = + + - 6 -圖3. 振型分解時程分析法程序框圖讀入數據計算動力荷載計算外荷載增量，然后計算擬靜力荷載增量計算幾何坐標系下的位移、速度和加速度，并輸出所需信息結束線性加速度法 =1.0, =0.5,=1/ 6Willson-法=1.42,=0.5,=1/ 6Newmark-法 =1.0, =1.0,=1.0計算時程分析控制參數調用振型分解時程分析法按METHOD 的值選擇1 2 3求出廣義質量矩陣、廣義剛度矩陣、廣義瑞利阻尼

19、矩陣及廣義干擾力對使用的振型階數循環求出初始振型位移、速度和加速度，并計算擬靜力剛度增量對時間步循環計算延長時間間隔的加速度、速度和位移增量計算下一時刻的振型位移、速度和加速度- 7 -2 2 Pt Pt ( yt yt) ( yt tyt)tM C = + + + +i ii i ii由式(6)求出延長時間步長上的位移增量y(t)后，可代入式(4)求出 y(t) ii。然后，可以用內插法求出對應于時間步長t 的加速度增量： y(t) 1 y(t) = ii ii。這樣，就可以求得時間步長t 上的速度和位移增量，進而得到下一時段的初始條件，并可以求得下一時段的位移、速度和加速度向量。重復以上步驟，便可求得各個時段的位移、速度和加速度增量。2.3 程序計算主要步驟本程序的具體計算步驟如下:(1). 讀入計算所需要的數據，主要包括地震波信息、剛度矩陣和集中