馬志強 樓云鋒 金先龍

摘要：多時間步長方法常用于求解有限元動力學空間多尺度問題，可分為單元重疊與非重疊兩種形式。單元非重疊形式的多時間步長方法涉及邊界數(shù)據(jù)的插值過程，造成計算不穩(wěn)定、精度不高。為了提高多時間步長方法計算穩(wěn)定性與求解精度，基于節(jié)點分割與重疊邊界單元，提出一種改進的多時間步長計算方法。結合顯式積分格式數(shù)據(jù)傳遞規(guī)律與重疊單元，分區(qū)邊界數(shù)據(jù)由動力學方程顯式求出，算法實施過程簡單。采用能量方法分析了算法的穩(wěn)定性能，給出了算法一般性穩(wěn)定條件。穩(wěn)定性分析結果表明改進方法的穩(wěn)定性只由分區(qū)內(nèi)部單元特性決定，重疊單元不影響算法的穩(wěn)定性能，這擴展了改進方法的適用范圍。數(shù)值算例結果表明該算法具有較高的計算精度，能量校驗進一步證實了算法的穩(wěn)定性能。

關鍵詞：多體系統(tǒng)動力學;顯式積分;多時間步長;穩(wěn)定性分析;能量校驗

中圖分類號：0313.7;0242.2

文獻標志碼：A

文章編號：1004-4523 （2019）06-1041-09

DOI：10. 16 385/j. cnki. issn. 1004-4523. 2019. 06. 013

引言

多體系統(tǒng)動力學以及車輛碰撞有限元仿真等問題中常涉及不同時間與空間尺度的耦合計算。這類問題需要根據(jù)結構設計特點、載荷分布以及分析區(qū)域劃分不同屬性的網(wǎng)格。采用中心差分法分析時整個模型采用單一時間步長。受限于穩(wěn)定性要求，臨界時間步長受模型最小單元屬性的限制，部分尺寸較小的網(wǎng)格會使得整個仿真模型采用較小的時間步長，這會導致計算效率大幅下降。

多時間步長方法或者子循環(huán)方法是處理這類問題的常用方法。所謂多時間步長方法是指根據(jù)需要將模型劃分成若干區(qū)域，不同區(qū)域根據(jù)分區(qū)單元特性選擇時間步長。一個多時間步長計算過程包含一個系統(tǒng)時間步與多個子循環(huán)時間步。

Liu和Belytschko首先采用網(wǎng)格分割將多時間步長方法應用于瞬態(tài)動力學問題求解[1]。而后眾多學者進一步研究改進方法，并將其運用到特定工程中。根據(jù)模型單元分割特點可以分為重疊與非重疊單元形式，不同多時間步長方法的區(qū)別在于分區(qū)邊界數(shù)據(jù)的處理方法。

非重疊單元形式的多時間步長方法研究中，以Smolinski和Daniel為代表的學者，在系統(tǒng)時間步內(nèi)，以位移、速度或者加速度滿足特定插值規(guī)律來處理各子區(qū)域邊界節(jié)點間的相互作用[2-4]。Gra-vouil和Combescure以FETI（the finite elementtearing and interconnecting method）方法為基礎，約束邊界節(jié)點速度在子循環(huán)時間步內(nèi)連續(xù)，提出了顯隱混合多時間步長的GC方法[5]。Prakash和Hj elmstad改進了GC方法，通過拉格朗日乘子約束邊界節(jié)點速度在系統(tǒng)時間步連續(xù)，提高了求解效率[6]。Lindsay等使用多時間步長方法處理諸如裂紋之類的結構非連續(xù)問題[7]。Ruparel等結合域分解方法與多時間步長方法處理非協(xié)調(diào)邊界網(wǎng)格有限元模型，通過能量分析給出算法穩(wěn)定的分區(qū)邊界數(shù)據(jù)連續(xù)性條件[8]。

國內(nèi)學者中，陳麗華等使用區(qū)域分解方法，構造了沖擊動力問題的混合時間步長顯式積分并行算法[9]。高暉等比較了幾種常用多時間步長方法的穩(wěn)定性與精度，提出了一種應用于汽車碰撞過程仿真的基于常速度的阻尼子循環(huán)法[10]。繆建成等基于中心差分方法提出一種適用于柔性多體系統(tǒng)仿真的顯式子循環(huán)方法[11]。

由上可知，以非重疊邊界為基礎的多時間步長方法，處理不同步長分區(qū)邊界數(shù)據(jù)傳遞，涉及邊界變量插值，這會惡化多時間步長方法的穩(wěn)定性能，限制分區(qū)步長比的選擇[12]?，F(xiàn)有的隱式或者顯隱子循環(huán)算法都無法保證無條件穩(wěn)定，邊界數(shù)據(jù)插值會降低算法的穩(wěn)定性，需要嚴格限制分區(qū)邊界數(shù)據(jù)的連續(xù)性要求，對于顯式子循環(huán)方法而言，還要限制各自分區(qū)內(nèi)部臨界時間步長。

同時，這類方法存在著變量不一致性誤差，進一步降低計算精度[13]。需要指出的是采用重疊單元形式的多時間步長方法的研究中，Arlequin方法是典型的方法[14]。Ghanem等將Arlequin方法用于結構瞬態(tài)動力學顯隱混合多尺度分析[15]，而后Fernier等研究顯式積分格式的Arlequin模型[16]。Arlequin方法使用拉格朗日乘子約束邊界連續(xù)性，拉格朗日乘子可以理解為不同分區(qū)的邊界耦合力。處理時間多尺度問題時Arlequin方法引入特定插值形式的拉格朗日乘子，這同非重疊形式的方法一樣需要嚴格限制邊界數(shù)據(jù)連續(xù)性條件。Rama等將剛度、質(zhì)量與阻尼矩陣拆分為塊對角矩陣與耦合矩陣，耦合的分區(qū)內(nèi)力通過重疊單元獲得[17]，但該方法只適用于分區(qū)同時間步長的計算。

為了改善現(xiàn)有多時間步長方法對邊界數(shù)據(jù)連續(xù)性的要求，提高不同分區(qū)步長比的選擇范圍，本文提出一種基于邊界單元重疊形式的顯式多時間步長計算方法。有限元仿真模型采用節(jié)點分割辦法被劃分為若干分區(qū)，分區(qū)邊界重疊多層單元。各分區(qū)采用預測校正形式的顯式Newmark積分格式。重疊邊界單元方法符合顯式動力學求解過程數(shù)據(jù)傳遞規(guī)律，重疊單元其實質(zhì)為相鄰分區(qū)預測波形的傳遞區(qū)域，因而處理分區(qū)多時間步長時不涉及邊界數(shù)據(jù)插值過程，其實質(zhì)是在子循環(huán)過程中，以重疊邊界的變化代替邊界數(shù)據(jù)的插值過程。

由于不同分區(qū)采用不同的時間步長，一般形式的多時間步長方法的穩(wěn)定性研究依舊是個熱點課題。本文基于能量法分析算法的穩(wěn)定性，給出了一般形式的穩(wěn)定性條件。最后以典型數(shù)值算例驗證方法的有效性。

1 顯式多時間步長計算方法

條件穩(wěn)定的顯式求解格式，其臨界時間步長受CFL條件（Courant-Friedrichs-Lewy condition）限制。為了保證顯式計算的數(shù)值穩(wěn)定性，采用的時間步長要小于臨時時間步長。所有分區(qū)最小時間步長為△tr，由分區(qū)單元屬性決定

特別地，當γE一γL一1/2時，滿足式（18）的穩(wěn)定性條件即含有阻尼形式的中心差分法的穩(wěn)定性條件。中心差分格式也適用于本文的多重疊網(wǎng)格形式的多時間步長計算方法。

采用能量增量形式計算的重疊網(wǎng)格多時間步方法。重疊分區(qū)矩陣BEB同樣滿足式（17）的分析結論，除一般顯式算法臨界步長限制外，無需額外約束重疊分區(qū)內(nèi)節(jié)點變量的連續(xù)性要求，就可以滿足穩(wěn)定性條件。

3.2 能量校驗方法

多時間步長方法中，數(shù)值計算過程中的不穩(wěn)定現(xiàn)象可以直觀的由能量校核方法體現(xiàn)。對線性彈性結構有限元分析模型中，模型的內(nèi)能W int、動能Wkin以及外力做的功wext可以表示為式中 K為剛度矩陣;M為質(zhì)量矩陣;u，v為節(jié)點位移與速度矢量;fext為等效節(jié)點外力矢量。若計算過程數(shù)值穩(wěn)定，則在系統(tǒng)時間步n內(nèi)，能量校驗滿足式中 ε為一很小的誤差容許極限（10-2）。通過檢測任意時間步有限元模型能量的變化，可進一步驗證穩(wěn)定性分析的結論。

4 數(shù)值算例

采用經(jīng)典的實體梁承受集中載荷模型為數(shù)值算例，梁兩端簡支。實體梁結構模型長L為20 m，高2m，厚1m。采用六面體網(wǎng)格離散，為了分析不同時間尺度對計算結果的影響，模型劃分了3種不同尺寸的有限元網(wǎng)格，離散后的實體梁模型有13720個單元，17370個節(jié)點，52110個自由度。實體梁模型有限元網(wǎng)格截面如圖2所示。

實體梁賦予各向同性材料，密度p=2200 kg/m3，彈性模量E=3×1010 Pa，泊松比v—0.22。梁沿著厚度方向劃分5層單元，點A所在梁截面上所有節(jié)點施加等效集中階躍載荷fext，載荷總大小為4×l06 N。仿真時間為0.4 s，采用顯式多時間步長方法計算，不同網(wǎng)格尺寸施加不同的時間積分步長，計算梁中間位置點B的橫向位移。模型的臨界時間步長受分區(qū)3單元尺寸的約束，臨界步長為1. 37×10-5 s。

不同分區(qū)采用的時間步長如表1所示。其中m為分區(qū)1與3步長比

采用步長比條件2即m=5計算梁中間點橫向位移，參考方法采用經(jīng)典的多時間步長方法[2]，采用模態(tài)疊加法計算梁中點響應的理論解。參考方法將顯式與隱式積分格式統(tǒng)一，在子循環(huán)過程中，小步長所需邊界節(jié)點位移實質(zhì)為相鄰大步長節(jié)點位移的線性插值。參考方法中求解格式采用預測校正New-mark格式，對大步長時間步邊界節(jié)點位移校正值做線性插值傳遞到小步長分區(qū)邊界節(jié)點。

圖3與4為本文所述方法與參考方法比較的結果，圖4為計算結果的局部放大圖。圖3與4可以看出經(jīng)典的多時間步長方法與本文方法在一定條件下都可以計算出梁橫向位移。從圖4中可以直觀地看出，在相同步長比的條件下，多重疊形式的網(wǎng)格允許動力學預測波形在整個分區(qū)內(nèi)傳遞，不涉及插值與波形截斷，參考方法中子循環(huán)過程邊界節(jié)點需要大步長分區(qū)節(jié)點的線性插值。相比于非重疊形式的參考解法，本文方法具有更高的計算精度。

分別以2倍步長比、5倍步長比與1 0倍步長比計算梁中間節(jié)點的位移，梁中間位置橫向位移計算結果局部放大圖如圖5所示。選取時間段0. 295-0.3 s，3條放大橫向位移曲線表明，隨著步長比的增加，梁中間位置橫向位移曲線基本保持一致，計算結果表明本文方法具有較高的精度一致性。

仿真時間為0. 297 s時，梁橫向位移計算云圖如圖6所示。

為了定量地分析不同算法位移計算結果，定義最大誤差MaxEr與計算誤差和SumEr兩個計算指標。其中最大誤差MaxEr=max‖Xm -Xref‖i，誤差和SumEr=∑ Xm - Xref‖i。其中n為時

間步數(shù)，Xref為某一變量參考值（理論值）。

由于分區(qū)采用的時間步長不同，在相同仿真時間下，時間步長越小，計算時間步越多，這樣就不具有相同的度量標準，因此全部采用步長比為m =2的條件下，大步長分區(qū)所需的時間步數(shù)作為標準n。統(tǒng)計梁中間點B橫向位移的不同步長比下的方法精度如表2所示。

從表2的計算結果可以看出，隨著步長比的增加，本文方法最大誤差與誤差和均逐步增加。也就說明給定網(wǎng)格劃分，采用步長比越大，計算精度相對降低。但是與經(jīng)典的非重疊多時間步長相比，本文所述的方法在相同計算步長的前提下具有較高的計算精度，同時步長比越大，本文方法所述計算精度保持更高。

采用本文3.2節(jié)能量計算方法計算模型在給定時間部內(nèi)的外力功、內(nèi)能、動能與檢驗能量，檢驗能量反映的是因為部分重疊網(wǎng)格產(chǎn)生的附加能量，這可以從側面檢驗穩(wěn)定性結論。圖7為10倍步長比條件下，實體梁模型能量分析結果。

在模型計算時間步內(nèi)，采用顯式多時間步長計算方法，當各個分區(qū)步長均滿足CFL條件時，分區(qū)產(chǎn)生檢驗能量相對于模型內(nèi)能、動能而言，數(shù)量級極小（10-4），未出現(xiàn)數(shù)值不穩(wěn)定現(xiàn)象。

為驗證步長比對計算效率的影響，統(tǒng)計實體梁在不同步長比下所需計算時間，計算資源采用普通微機CPU i7-4790K。構造統(tǒng)一的計算時間統(tǒng)計標準，此時分區(qū)1的時間步長為6. 25×10-5 s，分區(qū)3為1. 25×10-6 s，分區(qū)2步長變動。分區(qū)2與3的步長比為m。仿真時間為0.4 s。

由表3統(tǒng)計的模型計算時間可知：隨著步長比m的增加，分區(qū)2的時間步長逐漸由1. 25×10-6 s增加到2.5×10-5 s，分區(qū)2的計算時間步數(shù)相應減少，計算時間逐步減少，也應看到步長比增加到5以后，計算時間沒有明顯繼續(xù)下降，出現(xiàn)了負載不均衡現(xiàn)象。因而涉及多時間步長計算，負載均衡是要重點考慮的問題。

為了驗證時間步長效應對計算精度的影響，將圖2中分區(qū)2與分區(qū)3統(tǒng)一劃分較大尺寸網(wǎng)格，如圖8所示。梁的幾何尺寸與單元所采用材料屬性均保持不變，梁沿著厚度方向劃分3層網(wǎng)格。懸臂梁自由端所有節(jié)點施加正弦等效集中載荷F，載荷總大小為4×l06 N。載荷周期為0. 05 s，仿真時間為0.2 s。

整個模型含有1020個單元，1588個節(jié)點，顯式積分臨界時間步長受分區(qū)2網(wǎng)格尺寸限制，為4. 54×10-5 s，選取分區(qū)2的時間步長為4×10-5s，受限于分區(qū)2的臨界時間步長，此時步長比m最大可選為2。梁右端點上部節(jié)點C位移分別采用模態(tài)疊加法、本文所述方法以及顯式MGMT方法計算[8]。同步長與2倍步長條件下位移計算結果如圖9所示，參考的顯式MGMT方法采用2倍步長比，小分區(qū)步長4×10-5 s。局部放大如圖1 0所示。

由圖9中可知懸臂梁自由端受正弦載荷時，采用不同方法得到的自由端節(jié)點C位移得到大致相同的計算規(guī)律。2倍步長比MGMT求解不同分區(qū)邊界通過Mortar方法耦合節(jié)點外力，同F(xiàn)ETI方法相似，涉及插值過程。由圖1 0可以看出，本文所述重疊單元的方法相對與固定邊界的多時間步長方法，求解精度有所提高，同時提高步長比求解精度會降低。

5 結 論

多時間步長方法是處理局部網(wǎng)格多尺度問題的經(jīng)典方法，采用分區(qū)網(wǎng)格邊界重疊的方法，本文提出一種改進的顯式多時間步長求解方法。通過能量法分析了算法的穩(wěn)定性條件，給出了該方法一般性的穩(wěn)定性條件。該方法的優(yōu)點如下：

（1）在子循環(huán)時間步內(nèi)，預測波形得以完整的傳遞，不涉及插值、截斷等，提高計算精度;可以進一步擴大分區(qū)步長比使用范圍，縮短有限元空間多尺度問題的求解時間。

（2）重疊多層網(wǎng)格不影響分區(qū)穩(wěn)定性條件。這也意味著顯式子循環(huán)過程只要各自分區(qū)滿足臨界步長穩(wěn)定性要求，無需為重疊部分單元額外施加更加嚴格的穩(wěn)定性條件。這提高了多時間步長方法的適用范圍，同時方法易于編程實現(xiàn)。

（3）從邊界數(shù)據(jù)的傳遞規(guī)律來看，預測校正方法作為顯式積分方法的特例，子循環(huán)過程是以重疊邊界的變化代替邊界數(shù)據(jù)的插值過程。本文方法可以推廣到任意顯式時間積分格式（如中心差分法）。

由于顯式積分方法的解耦特性，該方法容易并行化求解大規(guī)模以及超大規(guī)模結構動力學問題。這將在提高精度的基礎上，進一步提升計算效率，這也是下一步研究的重點。

參考文獻：

[1] Liu W K， Belytschko T. Mixed-time implicit-explicitfinite elements for transient analysis[J]. Computers&Structures， 1982，15 （4）：445-450.

[2] Smolinski P，Belytschko T，Neal M. Multi-time-stepintegration using nodal partitioning[J]. InternationalJournal for Numerical Methods in Engineering， 1988，26（2）：349-359.

[3] Daniel W J T. Analysis and implementation of a newconstant acceleration subcycling algorithm[J]. Inter-national Journal for Numerical Methods in Engineer-ing， 1997，40（15）：2841-2855.

[4]Daniel W J T. A partial velocity approach to subcy-cling structural dynamics[J]. Computer Methods inApplied Mechanics and Engineering， 2003， 192 （3）：375-394.

[5] Gravouil A， Combescure A. Multi-time-step and two-scale domain decomposition method for nonlinearstructural dynamics[J]. International Journal for Nu-merical Methods in Engineering， 2003，58 （10）：1545-1569.

[6] Prakash A， Hjelmstad K D. A FETI-based multi-time-step coupling method for Newmark schemes instructural dynamics[J]. International Journal for Nu-merical Methods in Engineering， 2004， 61 （13）： 2183-2204.

[7] Lindsay P，Parks M L，Prakash A. Enabling fast，stable and accurate peridynamic computations usingmulti-time-step integration[J]. Computer Methods inApplied Mechanics and Engineering， 2016， 306： 382-405.

[8] Ruparel T， Eskandarian A， Lee J D. Concurrentmulti-domain simulations in structural dynamics usingmultiple grid and multiple time-scale （MGMT） meth-od[J]. International Journal of Computational Meth-ods， 2018，15（4）：1850021.

[9] 陳麗華，程建鋼，姚振漢，沖擊動力問題的混合積分并行算法及應用[J].工程力學，2003， 20（2）：15-20.

Chen Lihua， Cheng Jiangang， Yao Zhenhan. Mixedtime integration parallel algorithm and its applicaitionto dynamic impact problem[J]. Engineering Mechan-ics， 2003，20（2）： 15-20.

[10]高 暉，李光耀，鐘志華，等，汽車碰撞計算機仿真中的子循環(huán)法分析[J].機械工程學報，2005， 41（11）：9 8-101.

Gao Hui， Li Guangyao，Zhong Zhihua， et al. Analysisof subcycling algorithms for computer simulation ofcrashworthiness[J]. Chinese Journal of MechanicalEngineering，2005， 41（11）：98-101.

[11]繆建成，朱 平，陳關龍，等，多柔體系統(tǒng)響應計算的子循環(huán)計算方法研究[J].力學學報，2008， 40（4）：511- 519.

Miao Jiancheng， Zhu Ping， Chen Guanlong， et al.Study on sub-cycling algorithm for flexible multi-bodysystem[J]. Chinese Journal of Theoretical and AppliedMechanics， 2008， 40（4）：511-519.

[12] Klisinski M. Inconsistency errors of constant velocitymulti-time step integration algorithms[J]. ComputerAssisted Mechanics and Engineering Sciences， 2001，8（1）：121-139.

[13] Karimi S，Nakshatrala K B.On multi-time-step mon-olithic coupling algorithms for elastodynamics[J].Journal of Computational Physics， 2014， 273： 671-705.

[14] Dhia H B， Rateau G.The Arlequin method as a flexi-ble engineering design tool[J]. International Journalfor Numerical Methods in Engineering， 2005， 62（11）：1442-1462.

[15] Ghanem A， Torkhani M， Mahjoubi N， et al.Arlequinframework for multi-model， multi-time scale and het-erogeneous time integrators for structural transient dy-namics[J]. Computer Methods in Applied Mechaniesand Engineering， 2013， 254：292-308.

[16] Fernier A， Faucher V， Jamond 0. Multi-model Arle-uin method for transient structural dynamies with ex-plicit time integration[J]. International Journal for Numerical Methods in Engineering， 2017， 112 （9）：119 4-1215.

[17] Rao A R M， Rao T V S R A， Dattaguru B. A newparallel overlapped domain decomposition method fornonlinear dynamic finite element analysis[J]. Comput-ers& Structures， 2003，81（26-27）：2441-2454.

[18] Hughes T J R. The Finite Element Method： LinearStatic and Dynamic Finite Element Analysis [M].Courier Corporation， 2012.