馬志強(qiáng) 孔令爽 樓云鋒 金先龍

(1.上海交通大學(xué)機(jī)械與動(dòng)力工程學(xué)院， 上海 200240； 2.上海交通大學(xué)機(jī)械系統(tǒng)與振動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室， 上海 200240)

0 引言

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

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

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

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

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

(1)

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

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

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

(2)

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

Δt——積分時(shí)間步長(zhǎng)

(3)

1.1 顯式求解格式

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

(4)

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

1.2 隱式求解格式

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

(5)

其中

Keff=M+γΔtC+βΔt2K

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

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

2 交替格式的混合異步長(zhǎng)耦合計(jì)算方法

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

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

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

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

2.2 顯式小步長(zhǎng)分區(qū)計(jì)算

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

(6)

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

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

2.3 隱式大步長(zhǎng)分區(qū)計(jì)算

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

(7)

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

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

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

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

(8)

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

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

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

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

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

3 精度驗(yàn)證

采用彈簧質(zhì)量系統(tǒng)驗(yàn)證算法的計(jì)算精度與收斂性。5節(jié)點(diǎn)彈簧質(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),各個(gè)節(jié)點(diǎn)初始速度均為0。顯式與隱式分區(qū)參數(shù)均為β=0.5，γ=0.25。分區(qū)同步長(zhǎng)計(jì)算時(shí)節(jié)點(diǎn)3、4與彈簧k4為重疊區(qū)域。為了比較顯隱分區(qū)步長(zhǎng)比對(duì)精度的影響，異步長(zhǎng)計(jì)算時(shí)選擇步長(zhǎng)比m=3。此時(shí)節(jié)點(diǎn)2～5，彈簧3～5為重疊區(qū)域，節(jié)點(diǎn)m1為顯式分區(qū)內(nèi)部節(jié)點(diǎn)。集中參數(shù)的彈簧質(zhì)量系統(tǒng)的解析解作為計(jì)算結(jié)果的參考對(duì)象。圖中ui表示第i個(gè)節(jié)點(diǎn)的位移。

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

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

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

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

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

(9)

式中ui——i時(shí)刻5個(gè)節(jié)點(diǎn)位移向量

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

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

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

4 數(shù)值算例

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

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

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

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

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

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

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

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

5 結(jié)論

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

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

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