稀疏求解的動力學(xué)顯隱混合異步長交錯計算方法

馬志強 孔令爽 樓云鋒 金先龍

(1.上海交通大學(xué)機械與動力工程學(xué)院， 上海 200240； 2.上海交通大學(xué)機械系統(tǒng)與振動國家重點實驗室， 上海 200240)

0 引言

直接積分法是求解有限元結(jié)構(gòu)動力學(xué)問題的常用方法，可分為顯式與隱式兩大類。顯式方法適合求解沖擊載荷引起的波傳播問題，隱式方法適合求解結(jié)構(gòu)低頻振動問題[1]。復(fù)雜結(jié)構(gòu)動力學(xué)分析中常涉及不同時間與網(wǎng)格尺度，結(jié)合顯式與隱式方法優(yōu)點的顯隱混合算法是求解此類問題的經(jīng)典方法[2]。這類顯隱混合方法可以追溯到BELYTSCHKO等[3]的研究。在近40年的研究歷程中，大致可以分為2個階段。早期的研究以顯隱同步長為起點，結(jié)合顯式、隱式在求解動力學(xué)問題上的優(yōu)勢處理復(fù)雜結(jié)構(gòu)動力學(xué)問題[3-5]。國內(nèi)外學(xué)者嘗試將顯隱混合積分方法運用到分布式發(fā)電系統(tǒng)、流固耦合等領(lǐng)域[6-8]。節(jié)點分割和單元分割是兩種有限元模型的分區(qū)方法，顯隱分區(qū)邊界多采用共享節(jié)點形式[9-11]。顯式方法多采用中心差分法，隱式積分方法常采用Newmark方法。不同的顯隱混合方法多集中保證不同分區(qū)邊界數(shù)據(jù)連續(xù)性上，尤其對顯隱異步長而言，顯式分區(qū)計算時常需要對隱式分區(qū)邊界節(jié)點數(shù)據(jù)作插值處理[12]，這在一定程度上增加了算法對邊界數(shù)據(jù)的連續(xù)性要求[13]，也降低了顯隱最大步長比的選擇。

隨后，以FETI(Finite element tearing and interconnect method)方法為基礎(chǔ)，采用拉格朗日乘子耦合分區(qū)邊界的GC方法與PH方法被先后提出[14-15]。不同分區(qū)異步長插值處理由節(jié)點變量(位移、加速度等)變成拉格朗日乘子[16]，同樣存在上述邊界數(shù)據(jù)連續(xù)性問題。作為共享邊界節(jié)點的一種替換方案，重疊網(wǎng)格形式的Arlequin模型體現(xiàn)了在求解結(jié)構(gòu)動力系問題上的穩(wěn)定性優(yōu)勢[17]。

本文參考Arlequin結(jié)構(gòu)動力學(xué)模型，以重疊網(wǎng)格為基礎(chǔ)，提出一種改進的顯隱混合異步長計算方法。采用節(jié)點分割有限元模型，分區(qū)邊界節(jié)點與外部節(jié)點構(gòu)成耦合區(qū)域。

1 結(jié)構(gòu)動力學(xué)顯式與隱式計算方法

結(jié)構(gòu)經(jīng)有限元離散，含阻尼的線彈性體結(jié)構(gòu)動力學(xué)方程可以寫為[1]

(1)

式中M、C、K——結(jié)構(gòu)質(zhì)量、阻尼與剛度矩陣

fext——節(jié)點外力向量

阻尼矩陣C常用Rayleigh阻尼形式。為了統(tǒng)一顯式與隱式求解格式，采用基于預(yù)測校正格式的Newmark格式。預(yù)測校正格式求解過程分為預(yù)測步、加速度求解以及校正步3個求解過程。預(yù)測步中，第n+1時間步節(jié)點速度以及位移的預(yù)測值可以由第n時間步節(jié)點的位移、速度以及加速度表示為

(2)

β、γ——Newmark時間積分參數(shù)

Δt——積分時間步長

(3)

1.1 顯式求解格式

顯式求解過程是將時間離散后的位移與速度的預(yù)測值式(2)直接代入方程(1)中，此時節(jié)點加速度求解公式為[1]

(4)

質(zhì)量矩陣采用集中質(zhì)量矩陣形式，加速度求解過程無需進行矩陣求逆。

1.2 隱式求解格式

與顯式求解格式相兼容的隱式方法是將速度與位移的校正式(3)代入方程(1)中，隱式方法中的加速度求解寫成

(5)

其中

Keff=M+γΔtC+βΔt2K

式中Keff——求解節(jié)點加速度的等效剛度矩陣

由于剛度矩陣為稀疏非對角矩陣，節(jié)點加速度求解過程涉及等效剛度矩陣的求逆過程。式(4)與式(5)為格式兼容的顯式與隱式預(yù)測校正Newmark加速度求解格式，此格式用于顯隱混合求解的優(yōu)勢在于顯式與隱式求解可以融合在統(tǒng)一的程序中。

2 交替格式的混合異步長耦合計算方法

交替格式的顯隱式適用于多個分區(qū)的求解過程，為了方便描述，本節(jié)以兩個分區(qū)為例。有限元網(wǎng)格采用節(jié)點分割分成顯式與隱式兩個分區(qū)。與經(jīng)典的顯隱混合計算方法不同，交替格式的方法采用邊界重疊多重網(wǎng)格的方法實現(xiàn)異步長顯隱邊界數(shù)據(jù)的交換。

2.1 節(jié)點劃分的重疊網(wǎng)格

隱式分區(qū)采用大時間步長Δt1，顯式分區(qū)采用較小的子循環(huán)步長Δt2。假定隱式分區(qū)步長是顯式分區(qū)步長的整數(shù)m倍，即Δt1=mΔt2。圖1是以2倍步長比為例的重疊網(wǎng)格分區(qū)示意圖。隱式與顯式分區(qū)均包含內(nèi)部節(jié)點、邊界節(jié)點與外部節(jié)點，其中邊界節(jié)點與外部節(jié)點構(gòu)成重合區(qū)域。一個完整的系統(tǒng)時間步包含一個隱式時間步和多個顯式子循環(huán)時間步。內(nèi)部節(jié)點、邊界節(jié)點與外部節(jié)點用符號I、B、E表示。隱式分區(qū)的邊界節(jié)點即為顯式分區(qū)的外部節(jié)點，顯式分區(qū)的邊界節(jié)點為隱式分區(qū)的外部節(jié)點。

圖1 多重邊界網(wǎng)格顯隱分區(qū)示意圖Fig.1 Explicit-implicit partitioned schematics with multiple boundary mesh

2.2 顯式小步長分區(qū)計算

顯式計算過程中，求解某節(jié)點下一時刻數(shù)據(jù)只與當前時刻該節(jié)點的相鄰節(jié)點信息有關(guān)。在顯式子循環(huán)過程中，可正確求得的外部節(jié)點數(shù)據(jù)逐層遞減。子循環(huán)結(jié)束后顯式分區(qū)內(nèi)部與邊界節(jié)點數(shù)據(jù)可以正確求出。邊界節(jié)點加速度求解公式可以寫成

(6)

式中MB——與邊界節(jié)點相對應(yīng)的質(zhì)量矩陣

m——隱式分區(qū)與顯式分區(qū)的步長比

2.3 隱式大步長分區(qū)計算

隱式分區(qū)的節(jié)點加速度求解公式可以寫成

(7)

式中KE——隱式分區(qū)等效剛度矩陣Keff按照外部節(jié)點的分塊矩陣

KE(B+I)、K(B+I)E——隱式分區(qū)等效剛度矩陣Keff按照邊界節(jié)點的分塊矩陣

KB+I——等效剛度矩陣Keff與內(nèi)部節(jié)點對應(yīng)的分塊矩陣

對于隱式分區(qū)而言，隱式內(nèi)部與邊界節(jié)點加速度求解是域分解方法中回代求解過程[18]

(8)

采用稀疏直接求解器計算出隱式分區(qū)內(nèi)部節(jié)點與邊界節(jié)點加速度，賦值給顯式分區(qū)外部節(jié)點，完成一個系統(tǒng)時間步的顯式與隱式分區(qū)節(jié)點數(shù)據(jù)的求解。交替格式的混合異步長耦合計算方法求解順序為顯式子循環(huán)時間步、隱式分區(qū)回代求解以及顯式分區(qū)外部節(jié)點賦值3個串行交替步驟。顯隱異步長串行交替求解流程如圖2所示。虛線框內(nèi)是顯式分區(qū)節(jié)點與隱式分區(qū)節(jié)點加速度數(shù)據(jù)賦值過程。

圖2 顯隱異步長串行交替求解流程Fig.2 Sequential format process for explicit-implicit mixed multi-time step

2.4 隱式分區(qū)的稀疏行壓縮存儲格式的實現(xiàn)

顯式分區(qū)求解中質(zhì)量矩陣為對角矩陣，節(jié)點加速度直接根據(jù)節(jié)點自由度編號形成方程。方程(4)右端載荷項也是在單元層級計算按照節(jié)點自由度累加形成。采用坐標存儲對應(yīng)節(jié)點自由度質(zhì)量項與載荷項即可。

而對于隱式分區(qū)而言等效剛度矩陣為稀疏對稱矩陣，求解內(nèi)部與邊界節(jié)點加速度過程涉及矩陣求逆。采用稀疏矩陣的行壓縮CRS存儲格式有較好的存儲和求解效率[19]。隱式分區(qū)節(jié)點按照先外部節(jié)點后邊界與內(nèi)部節(jié)點的順序編號，將分區(qū)等效剛度矩陣與右端載荷項分塊。

3 精度驗證

采用彈簧質(zhì)量系統(tǒng)驗證算法的計算精度與收斂性。5節(jié)點彈簧質(zhì)量系統(tǒng)如圖3所示。采用的彈簧質(zhì)量參數(shù)為ki=1 N/m(i=1,2,…,6),mj=50 kg(j=1,2,…,5)。初始位移條件為u0=(0,0,1,0,0),各個節(jié)點初始速度均為0。顯式與隱式分區(qū)參數(shù)均為β=0.5，γ=0.25。分區(qū)同步長計算時節(jié)點3、4與彈簧k4為重疊區(qū)域。為了比較顯隱分區(qū)步長比對精度的影響，異步長計算時選擇步長比m=3。此時節(jié)點2～5，彈簧3～5為重疊區(qū)域，節(jié)點m1為顯式分區(qū)內(nèi)部節(jié)點。集中參數(shù)的彈簧質(zhì)量系統(tǒng)的解析解作為計算結(jié)果的參考對象。圖中ui表示第i個節(jié)點的位移。

圖3 5節(jié)點彈簧質(zhì)量系統(tǒng)Fig.3 Five-node spring mass system

顯式分區(qū)時間步長為0.001 s，分別采用顯隱同步長、3倍步長比以及GC 3倍步長比方法[14]計算節(jié)點3的位移與理論計算結(jié)果相比較。圖4為使用不同方法的位移計算結(jié)果。其中圖4a為節(jié)點3位移計算結(jié)果，因為節(jié)點5彈簧質(zhì)量系統(tǒng)初始條件為節(jié)點3的位移，節(jié)點2與節(jié)點4的理論位移應(yīng)該重疊，u2-u4可以直觀地反映不同方法下的位移計算誤差，如圖4b所示。

圖4 不同方法彈簧質(zhì)量系統(tǒng)位移節(jié)點結(jié)果Fig.4 Nodal displacement results of spring mass system with different methods

從圖4a可以看出，3種方法均可以計算出彈簧質(zhì)量系統(tǒng)的節(jié)點位移。由圖4b可知，顯式分區(qū)采用相同步長，隱式分區(qū)由同步長到3倍步長條件下，位移誤差有所增加，但是誤差在一個有限的區(qū)間內(nèi)波動。GC方法在計算節(jié)點位移時，涉及拉格朗日乘子處理，同時隨著仿真時間步的進行，位移誤差逐步增大。相同步長比條件下，本文所述方法具有較高的求解精度。

位移收斂特性采用相對位移計算結(jié)果[20]

(9)

式中ui——i時刻5個節(jié)點位移向量

統(tǒng)計顯式分區(qū)步長從0.001、0.002、0.006、0.01 s的位移誤差曲線，采用雙對數(shù)坐標系，統(tǒng)計時間為9 s，n取900，結(jié)果如圖5所示。

圖5 位移誤差收斂曲線Fig.5 Displacement error convergence rate curves for different methods

計算曲線結(jié)果表明，在顯隱分區(qū)積分參數(shù)均為β=0.5，γ=0.25的條件下，顯隱混合異步長計算方法具有2階的位移收斂精度。同顯式分區(qū)步長條件下，步長比m越大所得的計算誤差越大。

4 數(shù)值算例

為了進一步驗證算法在精度和效率方面的特性，本節(jié)采用U型管承受局部沖擊載荷fext作為數(shù)值算例。圖6為結(jié)構(gòu)有限元模型與外載荷時間歷程曲線。

圖6 U型管有限元模型及載荷時間曲線Fig.6 Finite element model and load time curve for U-type tube

模型左右對稱，計算時采用整體模型的一半。采用各向同性材料，材料彈性模型E為210 GPa，材料密度ρ為7 800 kg/m3,泊松比υ為0.28。U型管半徑為1 m，管壁厚度為8 mm。管中間部位承受彎曲載荷fext，載荷最大為1.5×105N。結(jié)構(gòu)采用三角形殼單元離散，沖擊部位與管彎曲部分網(wǎng)格尺度較小。顯式計算部分如圖6中所畫網(wǎng)格部位，重疊部分單元為管壁若干圈單元。模型含有13 574個節(jié)點，26 992個單元。

采用顯式積分方法時受限于穩(wěn)定性條件，采用積分時間步長Δt=2×10-6s。仿真時間0.1 s。隱式分區(qū)分別采用時間積分步長為mΔt。圖7為在步長比m為3、12這兩種情況下沖擊點豎向位移曲線。參考曲線為商業(yè)軟件LS-DYNA顯式計算結(jié)果。從圖7可以看出，隨著步長比的增加，算法計算結(jié)果符合位移計算規(guī)律，精度穩(wěn)定性較好。

圖7 不同方法計算節(jié)點豎向位移結(jié)果Fig.7 Vertical displacement results for impact node by different methods

為了研究算法計算效率，將本文方法與經(jīng)典顯式中心差分方法用于U型管的沖擊計算。模型在共享內(nèi)存模式計算機上計算，CPU主頻4.2 GHz，內(nèi)存16 GB。統(tǒng)計的有限元模型計算時間見表1。表中m表示隱式分區(qū)積分時間步長是顯式分區(qū)的倍數(shù)。

從表1可以看出，隨著步長比的增加，模型求解時間得以降低。由于采用串行計算格式，一個系統(tǒng)時間步內(nèi)顯式分區(qū)先計算，隱式分區(qū)后計算，提高顯隱分區(qū)步長比意味著在顯式分區(qū)步長不變的情況下，隱式分區(qū)采用了更大的時間步長。考慮重疊單元的邊界處理，在不同步長比下，隱式分區(qū)單元保持一致，顯式分區(qū)單元規(guī)模需增加部分分區(qū)邊界單元。計算時間的減少主要來自隱式分區(qū)計算步數(shù)的降低。也需要看到隨著顯隱步長比增加，時間比率減少幅度減少，此時隱式分區(qū)計算時間減少有限，計算時間占比中顯式計算所占比重逐漸增大。

表1 顯隱分區(qū)不同步長比計算時間Tab.1 Computational times for proposed explicit-implicit method with different time step ratios

5 結(jié)論

(1)分區(qū)邊界數(shù)據(jù)傳遞不涉及插值過程，這一改進提高了計算精度，積分參數(shù)β=0.5，γ=0.25，方法具有二階收斂精度。

(2)顯隱分區(qū)根據(jù)單元屬性選擇時間積分步長，降低計算時間。一定程度上，步長比越大，計算所需時間越小。

(3)顯隱分區(qū)采用兼容的Newmark格式，基于稀疏存儲CRS格式，顯隱計算程序可以具有統(tǒng)一的計算格式，這為顯隱式混合積分擴展至并行化提供了便利。