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

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

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

0 引言

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

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

本文參考Arlequin結構動力學模型，以重疊網格為基礎，提出一種改進的顯隱混合異步長計算方法。采用節點分割有限元模型，分區邊界節點與外部節點構成耦合區域。

1 結構動力學顯式與隱式計算方法

結構經有限元離散，含阻尼的線彈性體結構動力學方程可以寫為[1]

(1)

式中M、C、K——結構質量、阻尼與剛度矩陣

fext——節點外力向量

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

(2)

β、γ——Newmark時間積分參數

Δt——積分時間步長

(3)

1.1 顯式求解格式

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

(4)

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

1.2 隱式求解格式

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

(5)

其中

Keff=M+γΔtC+βΔt2K

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

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

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

交替格式的顯隱式適用于多個分區的求解過程，為了方便描述，本節以兩個分區為例。有限元網格采用節點分割分成顯式與隱式兩個分區。與經典的顯隱混合計算方法不同，交替格式的方法采用邊界重疊多重網格的方法實現異步長顯隱邊界數據的交換。

2.1 節點劃分的重疊網格

隱式分區采用大時間步長Δt1，顯式分區采用較小的子循環步長Δt2。假定隱式分區步長是顯式分區步長的整數m倍，即Δt1=mΔt2。圖1是以2倍步長比為例的重疊網格分區示意圖。隱式與顯式分區均包含內部節點、邊界節點與外部節點，其中邊界節點與外部節點構成重合區域。一個完整的系統時間步包含一個隱式時間步和多個顯式子循環時間步。內部節點、邊界節點與外部節點用符號I、B、E表示。隱式分區的邊界節點即為顯式分區的外部節點，顯式分區的邊界節點為隱式分區的外部節點。

圖1 多重邊界網格顯隱分區示意圖Fig.1 Explicit-implicit partitioned schematics with multiple boundary mesh

2.2 顯式小步長分區計算

顯式計算過程中，求解某節點下一時刻數據只與當前時刻該節點的相鄰節點信息有關。在顯式子循環過程中，可正確求得的外部節點數據逐層遞減。子循環結束后顯式分區內部與邊界節點數據可以正確求出。邊界節點加速度求解公式可以寫成

(6)

式中MB——與邊界節點相對應的質量矩陣

m——隱式分區與顯式分區的步長比

2.3 隱式大步長分區計算

隱式分區的節點加速度求解公式可以寫成

(7)

式中KE——隱式分區等效剛度矩陣Keff按照外部節點的分塊矩陣

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

KB+I——等效剛度矩陣Keff與內部節點對應的分塊矩陣

對于隱式分區而言，隱式內部與邊界節點加速度求解是域分解方法中回代求解過程[18]

(8)

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

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

2.4 隱式分區的稀疏行壓縮存儲格式的實現

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

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

3 精度驗證

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

圖3 5節點彈簧質量系統Fig.3 Five-node spring mass system

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

圖4 不同方法彈簧質量系統位移節點結果Fig.4 Nodal displacement results of spring mass system with different methods

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

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

(9)

式中ui——i時刻5個節點位移向量

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

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

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

4 數值算例

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

圖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。結構采用三角形殼單元離散，沖擊部位與管彎曲部分網格尺度較小。顯式計算部分如圖6中所畫網格部位，重疊部分單元為管壁若干圈單元。模型含有13 574個節點，26 992個單元。

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

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

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

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

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

5 結論

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

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

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