(1.大連理工大學(xué) 工程力學(xué)系 工業(yè)裝備結(jié)構(gòu)分析國家重點實驗室，大連 116024； 2.英特工程仿真技術(shù)(大連)有限公司，大連 116023)

1 引言

直接積分法是結(jié)構(gòu)瞬態(tài)分析中常用的方法，適應(yīng)范圍較廣[1]，適合比較復(fù)雜的非線性問題以及多物理場問題[2]。其中，無條件隱式的直接積分法由于允許采用較大的時間步長，具有比較重要的地位，如常用的Newmark法[3]、Wilson法[4]以及Houbolt法[5]等。在線性結(jié)構(gòu)動力學(xué)問題中，Bathe等[6]對其計算穩(wěn)定性進行了分析。

線性情況下無條件穩(wěn)定的積分算法在大變形結(jié)構(gòu)動力學(xué)計算中可能出現(xiàn)計算崩潰的問題[7]。為了提高算法的穩(wěn)定性，消除計算崩潰的問題，常用的解決方法是增大算法的數(shù)值阻尼。如在Newmark法中，可以通過調(diào)節(jié)參數(shù)β和γ來增大算法的數(shù)值阻尼。盡管引入數(shù)值阻尼提升了非線性計算的穩(wěn)定性，可是不合理的數(shù)值阻尼常導(dǎo)致算法精度的較大損失，如在計算結(jié)構(gòu)體的自由振動的過程中，振幅會出現(xiàn)明顯衰減。上述的Newmark法、Wilson法及Houbolt法等在用于結(jié)構(gòu)動力學(xué)分析時，為了計算穩(wěn)定性而引入的數(shù)值阻尼往往導(dǎo)致精度的損失過大。因此,如何合理地引入數(shù)值阻尼，即在保證計算穩(wěn)定性的情況下盡可能減小計算精度的損失，對時間積分算法意義重大。

Hilber等[8]提出的HHT方法有效地克服了以上方法的缺點，可以通過調(diào)整計算參數(shù)，減小數(shù)值阻尼對低階模態(tài)的影響。Chung等[9]進一步提出了Generalized-α法，廣泛應(yīng)用于計算力學(xué)中，表現(xiàn)了優(yōu)良的數(shù)值性能。彭海軍等[10,11]將HHT和Generalized-α法應(yīng)用到基于SiPESC平臺的多體系動力學(xué)仿真中。郭晛等[12]將HHT應(yīng)用到接觸約束多體系統(tǒng)動力學(xué)仿真中。另外，國外的商業(yè)軟件ABAQUS和ANSYS在結(jié)構(gòu)動力學(xué)有限元分析中均采用了HHT法。近年來，針對大變形動力學(xué)的時間積分的數(shù)值不穩(wěn)定問題，Bathe等[7]提出了一種無參數(shù)控制的時間積分方法，是 Newmark 法和Euler法的一種組合方法。劉天云等[13]也提出了一種無參數(shù)控制的Euler法。

預(yù)估-校正方法在直接積分方法中是一種比較常用的手段[14,15]。首先預(yù)估出下一個時間步上的解，然后通過非線性迭代，計算增量并校正預(yù)估解，獲得最終的解。將預(yù)估-校正方法應(yīng)用到非線性結(jié)構(gòu)動力學(xué)求解的Generalized-α法中，可以在保證Generalized-α法的優(yōu)越性能前提下，簡化非線性迭代公式，并易于編程實現(xiàn)。本文基于多物理場仿真平臺INTESIM[16]實現(xiàn)了一種基于預(yù)估-校正的Generalized-α法，并應(yīng)用到非線性結(jié)構(gòu)動力學(xué)的有限元分析中。通過殼和實體單元的結(jié)構(gòu)大變形動力學(xué)算例，與其他的時間積分方法進行了比較，說明了該方法的優(yōu)勢。

2 基于Generalized-α的預(yù)估-校正算法

2.1 基本控制方程

結(jié)構(gòu)動力學(xué)方程經(jīng)過有限元離散后具有以下形式，

Ma+Cv+Q=F

(1)

式中M為質(zhì)量矩陣，C為阻尼矩陣，F(xiàn)為外力向量，Q為內(nèi)力列向量，a為加速度列向量，v為速度列向量。對于線性問題，內(nèi)力可以寫為Q=Kd。但是對于非線性問題，Q和位移d是非線性關(guān)系，假設(shè)

Q=Q(d)

(2)

對Q在d0處進行泰勒展開，并截取一次項，代入式(1)，得到

Ma+Cv+T(d0)Δd=F-Q(d0)

(3)

式中Δd=d-d0，T(d0)為d0處的切線剛度矩陣，Q(d0)為d0對應(yīng)的內(nèi)力。

2.2 時間離散格式

(1) 在第n到n+1時刻內(nèi)，Newmark積分方法對位移和速度作如下假設(shè)，

(4)

vn + 1=vn+[(1-γ)an+γan + 1]Δt

(5)

式中Δt為時間步長，β和γ是時間積分參數(shù)。γ=0.5，β=0.25對應(yīng)常平均加速度法。Newmark法的穩(wěn)定性條件是γ≥0.5，β≥0.25(0.5+γ)2。γ和β越大，Newmark法的數(shù)值耗散效應(yīng)越強。

根據(jù)非線性動力學(xué)方程(1)，Newmark法在第n+1時刻的當(dāng)前迭代步i建立離散變量的動力學(xué)方程：

(6)

相應(yīng)地，根據(jù)非線性動力學(xué)方程(3)，建立非線性迭代格式：

(7)

(2) Chung等[9]提出的Generalized-α方法也屬于Newmark家族類的方法。不同于傳統(tǒng)的 New-mark 法，其在半離散點上建立動力學(xué)方程，得到

(8)

(9)

(10)

(11)

(12)

tn + 1 - αf=(1-αf)tn + 1+αftn

(13)

αm和αf為新引入的參數(shù)，積分參數(shù)一般取為

(14)

(15，16)

2.3 基于預(yù)估-校正的Generalized-α法

在第n+1時間步時，已知第n個時間步的計算結(jié)果an，vn和dn，現(xiàn)在需要計算第n+1步的an + 1，vn + 1和dn + 1。

(1) 在非線性迭代求解之前，即i=0時(i表示非線性迭代的步數(shù))，預(yù)估第n+1步的加速度初始解為

(17)

將式(17)代入式(4，5)，得到位移和速度的初始預(yù)估解為

(18)

(19)

(20)

(21)

(22)

最后，結(jié)合式(9～12，20～22)，代入離散方程(8)，得到

(23)

根據(jù)式(9～11)，可以進一步將式(23)寫成

[M(1-αm)+C(1-αf)γΔt+

(24)

式(24)可以寫為

(25)

其中殘差為

(26)

3 數(shù)值算例

算例1薄板受迫振動問題

圖1 薄板模型尺寸及有限元模型

Fig.1 Module size of thin plate and finite element model

圖2 薄板受到的沖擊載荷的時間履歷

Fig.2 Time history of the transient load on the sheet

通過對比四種方法計算結(jié)果的位移響應(yīng)曲線可以看出，四種情況的計算結(jié)果具有一致性，但是當(dāng)參數(shù)取為γ=0.5，β=0.25時，Newmark方法在0.52 s時出現(xiàn)計算崩潰(圖中方塊形標(biāo)記)，Bathe方法在0.23 s時出現(xiàn)了計算崩潰(圖中三角形標(biāo)記)，而Generalized-α方法順利地完成計算，沒有出現(xiàn)計算崩潰的問題。

可以看出，當(dāng)參數(shù)取為γ=0.6，β=0.3025時，Newmark方法也順利地完成了計算，沒有出現(xiàn)計算崩潰的問題，但是由于這種方法引入較大的數(shù)值阻尼，致使板自由振動的振幅發(fā)生較明顯的衰減，產(chǎn)生較大的誤差。

改變工況，將圖2的沖擊載荷的峰值降低到F=50fN/m2(其他保持不變)，并將計算時間延長到5 s，采用Newmark方法(γ=0.5，β=0.25)、Bathe方法和Generalized-α方法分別對該問題重新計算，得到檢測點處垂直位移響應(yīng)曲線如圖4所示?？梢娫谳d荷減小之后，Newmark方法(γ=0.5，β=0.25)和Bathe法穩(wěn)定計算的時間延長了，但是Newmark法在4.7 s時還是出現(xiàn)了計算崩潰。相比之下，Generalized-α法和Bathe法計算保持穩(wěn)定，沒有出現(xiàn)計算崩潰的問題。

圖3 四種情況位移響應(yīng)曲線對比(F=500fN/m2，Δt=0.002，t=2 s)

Fig.3 Comparison of displacement response curves of four cases

圖4 增加計算時間位移響應(yīng)曲線對比(F=50fN/m2，Δt=0.002，t=5 s)

Fig.4 Comparison of displacement response curve of increasing the calculation time

進一步考察不同時間步長Δt對三種方法的影響。由表1可知，當(dāng)時間步長較大時，Newmark法和Bathe法均會出現(xiàn)計算崩潰，而Generalized-α方法相對穩(wěn)定，沒有出現(xiàn)計算崩潰的問題。

算例2曲柄結(jié)構(gòu)受迫振動問題

表1 不同增量步長三種方法穩(wěn)定性對比(50fN/m2，2 s)

Tab.1 Stability comparison of different time incremental step in three different methods

Δt/s預(yù)估-矯正的Generalized-α(ρ=0.6)Newmark(γ=0.5,β=0.25)Bathe0.01穩(wěn)定計算崩潰計算崩潰0.0075穩(wěn)定計算崩潰穩(wěn)定0.005穩(wěn)定穩(wěn)定穩(wěn)定

圖5 曲柄結(jié)構(gòu)模型尺寸及邊界條件

Fig.5 Crank model size and boundary conditions

可以看出，Generalized-α法沒有出現(xiàn)計算崩潰的現(xiàn)象，Bathe法在計算時間為0.0065 s時出現(xiàn)計算崩潰現(xiàn)象(圖中三角形標(biāo)記)，當(dāng)參數(shù)取為γ=0.5，β=0.25時，Newmark法在計算時間為0.019 s時出現(xiàn)計算崩潰現(xiàn)象(圖中方塊形標(biāo)記)。另外，當(dāng)改變Newmark法的參數(shù)為γ=0.6，β=0.3025，雖然增加了Newmark法的穩(wěn)定性，但是數(shù)值耗散效應(yīng)明顯，精度較差。

圖6 四種情況位移響應(yīng)曲線對比(F=10×106fN/m2，Δt=0.0003，t=0.03 s)

Fig.6 Comparison of displacement response curves of four cases

圖7 增加計算時間后三種方法位移響應(yīng)對比(F=5×106fN/m2，Δt=0.0004，t=0.1 s)

Fig.7 Comparison of displacement response of thre methods when increasing computing time

表2 不同時間步長下三種方法的穩(wěn)定性對比(F=5×106fPa，t=0.1 s)

Tab.2 Stability comparison of different time incremental step in three different methods

Δt/s預(yù)估-矯正的Generalized-α(ρ=0.6)Newmark(γ=0.5,β=0.25)Bathe0.0005穩(wěn)定計算崩潰計算崩潰0.0004穩(wěn)定計算崩潰數(shù)值耗散明顯0.0003穩(wěn)定穩(wěn)定穩(wěn)定

4 結(jié) 論

本文提出一種基于預(yù)估-校正的Generalized-α法，基于INTESIM平臺實現(xiàn)該算法，并將其應(yīng)用于非線性結(jié)構(gòu)動力學(xué)問題。通過數(shù)值算例，與 Newmark 法和Bathe法進行了比較，說明了本文方法具有較好的穩(wěn)定性和精度。