基于預(yù)估-校正的Generalized-α法及其在非線性結(jié)構(gòu)動(dòng)力學(xué)中的應(yīng)用

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

1 引言

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

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

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

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

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

2.1 基本控制方程

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

Ma+Cv+Q=F

(1)

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

Q=Q(d)

(2)

對(duì)Q在d0處進(jìn)行泰勒展開，并截取一次項(xiàng)，代入式(1)，得到

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

(3)

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

2.2 時(shí)間離散格式

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

(4)

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

(5)

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

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

(6)

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

(7)

(2) Chung等[9]提出的Generalized-α方法也屬于Newmark家族類的方法。不同于傳統(tǒng)的 New-mark 法，其在半離散點(diǎn)上建立動(dòng)力學(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時(shí)間步時(shí)，已知第n個(gè)時(shí)間步的計(jì)算結(jié)果an，vn和dn，現(xiàn)在需要計(jì)算第n+1步的an + 1，vn + 1和dn + 1。

(1) 在非線性迭代求解之前，即i=0時(shí)(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)，可以進(jìn)一步將式(23)寫成

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

(24)

式(24)可以寫為

(25)

其中殘差為

(26)

3 數(shù)值算例

算例1薄板受迫振動(dòng)問題

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

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

圖2 薄板受到的沖擊載荷的時(shí)間履歷

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

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

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

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

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

Fig.3 Comparison of displacement response curves of four cases

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

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

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

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

表1 不同增量步長(zhǎng)三種方法穩(wěn)定性對(duì)比(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)定計(jì)算崩潰計(jì)算崩潰0.0075穩(wěn)定計(jì)算崩潰穩(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)計(jì)算崩潰的現(xiàn)象，Bathe法在計(jì)算時(shí)間為0.0065 s時(shí)出現(xiàn)計(jì)算崩潰現(xiàn)象(圖中三角形標(biāo)記)，當(dāng)參數(shù)取為γ=0.5，β=0.25時(shí)，Newmark法在計(jì)算時(shí)間為0.019 s時(shí)出現(xiàn)計(jì)算崩潰現(xiàn)象(圖中方塊形標(biāo)記)。另外，當(dāng)改變Newmark法的參數(shù)為γ=0.6，β=0.3025，雖然增加了Newmark法的穩(wěn)定性，但是數(shù)值耗散效應(yīng)明顯，精度較差。

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

Fig.6 Comparison of displacement response curves of four cases

圖7 增加計(jì)算時(shí)間后三種方法位移響應(yīng)對(duì)比(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 不同時(shí)間步長(zhǎng)下三種方法的穩(wěn)定性對(duì)比(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)定計(jì)算崩潰計(jì)算崩潰0.0004穩(wěn)定計(jì)算崩潰數(shù)值耗散明顯0.0003穩(wěn)定穩(wěn)定穩(wěn)定

4 結(jié) 論

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