考慮邊界穩(wěn)定的結(jié)構(gòu)振動(dòng)顯式多時(shí)間步長(zhǎng)計(jì)算方法

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

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

關(guān)鍵詞：多體系統(tǒng)動(dòng)力學(xué);顯式積分;多時(shí)間步長(zhǎng);穩(wěn)定性分析;能量校驗(yàn)

中圖分類(lèi)號(hào)：0313.7;0242.2

文獻(xiàn)標(biāo)志碼：A

文章編號(hào)：1004-4523 （2019）06-1041-09

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

引言

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

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

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

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

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

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

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

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

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

1 顯式多時(shí)間步長(zhǎng)計(jì)算方法

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

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

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

3.2 能量校驗(yàn)方法

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

4 數(shù)值算例

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

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

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

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

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

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

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

為了定量地分析不同算法位移計(jì)算結(jié)果，定義最大誤差MaxEr與計(jì)算誤差和SumEr兩個(gè)計(jì)算指標(biāo)。其中最大誤差MaxEr=max‖Xm -Xref‖i，誤差和SumEr=∑ Xm - Xref‖i。其中n為時(shí)

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

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

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

采用本文3.2節(jié)能量計(jì)算方法計(jì)算模型在給定時(shí)間部?jī)?nèi)的外力功、內(nèi)能、動(dòng)能與檢驗(yàn)?zāi)芰?，檢驗(yàn)?zāi)芰糠从车氖且驗(yàn)椴糠种丿B網(wǎng)格產(chǎn)生的附加能量，這可以從側(cè)面檢驗(yàn)穩(wěn)定性結(jié)論。圖7為10倍步長(zhǎng)比條件下，實(shí)體梁模型能量分析結(jié)果。

在模型計(jì)算時(shí)間步內(nèi)，采用顯式多時(shí)間步長(zhǎng)計(jì)算方法，當(dāng)各個(gè)分區(qū)步長(zhǎng)均滿(mǎn)足CFL條件時(shí)，分區(qū)產(chǎn)生檢驗(yàn)?zāi)芰肯鄬?duì)于模型內(nèi)能、動(dòng)能而言，數(shù)量級(jí)極?。?0-4），未出現(xiàn)數(shù)值不穩(wěn)定現(xiàn)象。

為驗(yàn)證步長(zhǎng)比對(duì)計(jì)算效率的影響，統(tǒng)計(jì)實(shí)體梁在不同步長(zhǎng)比下所需計(jì)算時(shí)間，計(jì)算資源采用普通微機(jī)CPU i7-4790K。構(gòu)造統(tǒng)一的計(jì)算時(shí)間統(tǒng)計(jì)標(biāo)準(zhǔn)，此時(shí)分區(qū)1的時(shí)間步長(zhǎng)為6. 25×10-5 s，分區(qū)3為1. 25×10-6 s，分區(qū)2步長(zhǎng)變動(dòng)。分區(qū)2與3的步長(zhǎng)比為m。仿真時(shí)間為0.4 s。

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

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

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

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

5 結(jié) 論

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

（1）在子循環(huán)時(shí)間步內(nèi)，預(yù)測(cè)波形得以完整的傳遞，不涉及插值、截?cái)嗟?，提高?jì)算精度;可以進(jìn)一步擴(kuò)大分區(qū)步長(zhǎng)比使用范圍，縮短有限元空間多尺度問(wèn)題的求解時(shí)間。

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

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

由于顯式積分方法的解耦特性，該方法容易并行化求解大規(guī)模以及超大規(guī)模結(jié)構(gòu)動(dòng)力學(xué)問(wèn)題。這將在提高精度的基礎(chǔ)上，進(jìn)一步提升計(jì)算效率，這也是下一步研究的重點(diǎn)。

參考文獻(xiàn)：

[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] 陳麗華，程建鋼，姚振漢，沖擊動(dòng)力問(wèn)題的混合積分并行算法及應(yīng)用[J].工程力學(xué)，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]高 暉，李光耀，鐘志華，等，汽車(chē)碰撞計(jì)算機(jī)仿真中的子循環(huán)法分析[J].機(jī)械工程學(xué)報(bào)，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]繆建成，朱 平，陳關(guān)龍，等，多柔體系統(tǒng)響應(yīng)計(jì)算的子循環(huán)計(jì)算方法研究[J].力學(xué)學(xué)報(bào)，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.