馬志強(qiáng), 樓云鋒, 李俊杰, 金先龍
(上海交通大學(xué) 機(jī)械與動(dòng)力工程學(xué)院; 機(jī)械系統(tǒng)與振動(dòng)國(guó)家重點(diǎn)實(shí)驗(yàn)室, 上海 200240)
大型復(fù)雜結(jié)構(gòu)的有限元分析常涉及局部沖擊載荷、接觸等問(wèn)題,需要針對(duì)特定分析部位劃分精細(xì)化的網(wǎng)格[1].這類(lèi)問(wèn)題常采用顯式算法求解,但受限于求解精度以及局部單元穩(wěn)定性的要求,整個(gè)模型需要采用較小的時(shí)間步長(zhǎng),這會(huì)導(dǎo)致計(jì)算時(shí)間大大加長(zhǎng)[2].因此國(guó)內(nèi)外眾多學(xué)者提出改進(jìn)方案,期望可根據(jù)結(jié)構(gòu)區(qū)域的特點(diǎn),將結(jié)構(gòu)分成不同屬性的區(qū)域,針對(duì)單一區(qū)域選擇合適的積分算法和時(shí)間步長(zhǎng)進(jìn)行動(dòng)力學(xué)分析.
改進(jìn)方法主要包括域分解法以及子循環(huán)方法兩大類(lèi).域分解法是將計(jì)算區(qū)域分解為若干小的子區(qū)域,各個(gè)子區(qū)域單獨(dú)計(jì)算.根據(jù)分區(qū)邊界是否重疊可以將域分解法分為重疊與非重疊兩種形式.重疊形式主要包括Schwarz迭代法和Arlequin方法[3],非重疊形式包括Schur補(bǔ)方法和dual Schur方法,其區(qū)別在于處理區(qū)域邊界約束的方式.張友良等[4]采用Schur補(bǔ)方法開(kāi)發(fā)了巖土大規(guī)模并行計(jì)算軟件.Farhat等[5]提出了以FETI (Finite Element Tearing and Interconnecting)為代表的dual Schur方法.劉鋮等[6]基于FETI方法研究大變形柔性多體系統(tǒng)動(dòng)力學(xué)方程高效求解算法.Kozubek等[7]將FETI方法推廣到大規(guī)模非線性動(dòng)力學(xué)并行仿真.
子循環(huán)算法是在一次計(jì)算步內(nèi)包含兩種時(shí)步即一次主步計(jì)算和多次子步計(jì)算的積分算法,主要用于解決分區(qū)步長(zhǎng)不匹配的問(wèn)題.早期研究者主要有Smolinski[8]與Belyschko[9],研究對(duì)象多針對(duì)兩分區(qū)的結(jié)構(gòu)動(dòng)力學(xué)問(wèn)題.Prakash等[10]研究了非線性結(jié)構(gòu)的異步長(zhǎng)計(jì)算,通過(guò)拉格朗日乘子約束分區(qū)邊界速度連續(xù).劉占生等[11]將子循環(huán)與預(yù)測(cè)校正方法相結(jié)合求解壁板動(dòng)力響應(yīng).高雙武等[12]采用子循環(huán)的緊耦合算法,分析固體火箭發(fā)動(dòng)機(jī)點(diǎn)火瞬態(tài)工作過(guò)程.
在上述研究中,域分解法可在一定程度上提高加速比,縮短計(jì)算時(shí)間,但不同計(jì)算區(qū)域采用相同的時(shí)間步長(zhǎng),在處理存在局部細(xì)分網(wǎng)格的動(dòng)力學(xué)問(wèn)題時(shí),會(huì)造成計(jì)算資源的浪費(fèi).在經(jīng)典的顯式子循環(huán)方法中,分區(qū)多采用單元分割,不同步長(zhǎng)邊界節(jié)點(diǎn)信息傳遞過(guò)程常涉及插值或截?cái)郲13],這在一定程度上造成計(jì)算不穩(wěn)定、降低計(jì)算精度,也會(huì)限制最大步長(zhǎng)比的選取[14].
為此,本文提出一種基于多重節(jié)點(diǎn)的顯式異步長(zhǎng)并行計(jì)算方法.該方法將多重節(jié)點(diǎn)引入異步長(zhǎng)并行計(jì)算,顯式小步長(zhǎng)分區(qū)子循環(huán)過(guò)程中預(yù)測(cè)波形得以在顯式重疊區(qū)域內(nèi)完整傳遞,沒(méi)有插值、截?cái)嗾`差,計(jì)算精度較高,不同網(wǎng)格區(qū)域可以采用較大的步長(zhǎng)比.顯式算法具有解耦特性(求解下一時(shí)刻的相關(guān)信息只需用到相鄰節(jié)點(diǎn)上一時(shí)步的相關(guān)信息),使并行計(jì)算過(guò)程不需要聯(lián)立求解界面方程.因此,該方法具有一般顯式并行算法易于編程實(shí)現(xiàn)的特點(diǎn),同時(shí)相對(duì)于經(jīng)典的子循環(huán)算法精度較高.
含阻尼的線彈性體結(jié)構(gòu)動(dòng)力學(xué)方程可以寫(xiě)為
(1)
顯式Newmark求解格式分為預(yù)測(cè)步與校正步.預(yù)測(cè)步中,第n+1時(shí)間步節(jié)點(diǎn)位移以及速度的預(yù)測(cè)值可以由第n時(shí)間步節(jié)點(diǎn)的位移、速度以及加速度表示[15]:
(2)
(3)
(4)
將求得的節(jié)點(diǎn)加速度值帶入顯式Newmark離散格式的校正步
(5)
結(jié)構(gòu)有限元網(wǎng)格采用節(jié)點(diǎn)分割方法,將節(jié)點(diǎn)劃分給不同分區(qū),分割單元被相鄰分區(qū)共享,并行計(jì)算過(guò)程需要共享分割單元信息.定義分區(qū)步長(zhǎng)比:
(6)
式中:Δti為第i個(gè)分區(qū)選用的時(shí)間步長(zhǎng);Δtr為所有分區(qū)最小時(shí)間步長(zhǎng).異步長(zhǎng)計(jì)算過(guò)程中涉及主時(shí)間步與子循環(huán)時(shí)間步,定義時(shí)間變量上標(biāo):
(7)
式中:X為某一變量;n為主時(shí)間步數(shù);j為子循環(huán)時(shí)間步數(shù).采用節(jié)點(diǎn)分割策略,定義節(jié)點(diǎn)布爾分割向量Oi,如果某節(jié)點(diǎn)在給定分區(qū)i,則該節(jié)點(diǎn)在分割向量中為1,否則為0.含多重節(jié)點(diǎn)的分區(qū)i布爾分割向量定義為:
(8)
則所有節(jié)點(diǎn)組成的分割向量O=INC.INC中所有分量都為1,NC為有限元節(jié)點(diǎn)自由度數(shù)目.
顯式異步長(zhǎng)求解過(guò)程需要重疊多層邊界節(jié)點(diǎn).圖1是以2倍步長(zhǎng)比的網(wǎng)格劃分示意圖.圓形表示內(nèi)部節(jié)點(diǎn),用I表示;三角形節(jié)點(diǎn)表示邊界節(jié)點(diǎn),用B表示;正方形節(jié)點(diǎn)表示外部節(jié)點(diǎn),用E表示.陰影部分為共享單元.由節(jié)點(diǎn)重疊和分割單元相鄰分區(qū)共享可知:分區(qū)邊界節(jié)點(diǎn)是相鄰分區(qū)的外部節(jié)點(diǎn).
圖1 異步長(zhǎng)計(jì)算多重節(jié)點(diǎn)劃分示意
參照分區(qū)節(jié)點(diǎn)分割向量的定義式,小步長(zhǎng)分區(qū)分割向量為
Oi=
(9)
式中:OI代表內(nèi)部節(jié)點(diǎn)分割向量(包含邊界節(jié)點(diǎn));OE為外部節(jié)點(diǎn)分割向量.對(duì)步長(zhǎng)比為mi的分區(qū)(為簡(jiǎn)化公式暫不包含阻尼),離散的結(jié)構(gòu)動(dòng)力學(xué)方程可以寫(xiě)為
(10)
對(duì)于小步長(zhǎng)分區(qū)的一個(gè)子循環(huán)時(shí)間步j(luò),內(nèi)部節(jié)點(diǎn)與邊界節(jié)點(diǎn)的加速度可以寫(xiě)為
(11)
(12)
對(duì)于大步長(zhǎng)分區(qū),內(nèi)部與邊界節(jié)點(diǎn)的加速度求解公式與式(10)和(11)相似,不同的是此時(shí)分區(qū)時(shí)間步長(zhǎng)為Δti=miΔtr,預(yù)測(cè)步與校正步中求解公式需要將時(shí)間步長(zhǎng)換成miΔtr.大步長(zhǎng)分區(qū)質(zhì)量矩陣與剛度矩陣需按照分區(qū)分割向量分塊.
算法的有效性在國(guó)家超級(jí)計(jì)算廣州中心的“天河二號(hào)”超級(jí)計(jì)算機(jī)上進(jìn)行測(cè)試,文中算例均采用CPU核心計(jì)算.并行算法的性能常用加速比與并行計(jì)算效率來(lái)進(jìn)行衡量.異步長(zhǎng)并行算法中,結(jié)構(gòu)分區(qū)數(shù)以及步長(zhǎng)比不同都會(huì)導(dǎo)致計(jì)算加速比不同.定義結(jié)構(gòu)分區(qū)數(shù)目k固定,步長(zhǎng)比m不同時(shí)加速比為
Ski=tkm/tki
(13)
定義m固定而k不同時(shí)加速比為
Smi=tmk/tmi
(14)
計(jì)算效率為衡量一個(gè)處理器的計(jì)算能力被有效利用的百分比,定義為
E=S/p
(15)
式中:S為加速比;p為使用計(jì)算核心數(shù)目.
圖2 彈性桿受軸向集中載荷
圖2為彈性桿在軸向集中荷載作用下結(jié)構(gòu)示意圖.桿左端固定,右端點(diǎn)A受軸向集中載荷.桿材料密度ρ為 1 000 kg/m3,彈性模量1 GPa,橫截面積 0.001 m2,桿的長(zhǎng)度l為1 m,自由端加載外力fext為 1 000 N.x軸代表?xiàng)U軸向坐標(biāo).整個(gè)結(jié)構(gòu)劃分為4個(gè)子分區(qū),每個(gè)子分區(qū)長(zhǎng)度為 0.25 m,每個(gè)子分區(qū)離散為200個(gè)兩節(jié)點(diǎn)桿單元,整個(gè)結(jié)構(gòu)有800個(gè)單元,801個(gè)節(jié)點(diǎn).仿真時(shí)間 0.05 s.
子分區(qū)1、2、3采用相同時(shí)間步長(zhǎng),步長(zhǎng)Δt為 2.5 μs.子分區(qū)4依次使用相同時(shí)間步長(zhǎng)Δtr(Δt=Δtr)、5倍時(shí)間步長(zhǎng)(Δt=5Δtr)與10倍時(shí)間步長(zhǎng)(Δt=10Δtr)共3種情況.
圖3為子分區(qū)4采用不同時(shí)間步長(zhǎng)比,桿右端點(diǎn)A軸向位移(D)曲線與理論位移曲線.圖4為圖3的局部放大圖.圖3和4表明, 隨著分區(qū)時(shí)間步長(zhǎng)比的增加,位移的計(jì)算精度沒(méi)有出現(xiàn)明顯下降,依舊保持較高精度.
圖3 右端點(diǎn)A軸向位移曲線
為比較網(wǎng)格特性對(duì)計(jì)算精度的影響,將子分區(qū)4進(jìn)行細(xì)分離散,離散2倍,此時(shí)含有400個(gè)單元,而子分區(qū)1~3維持200個(gè)單元不變.子分區(qū)1~3仍采用步長(zhǎng)Δt為 2.5 μs,子分區(qū)4采用步長(zhǎng)比為10的異步長(zhǎng)方法,計(jì)算桿右端點(diǎn)A的軸向位移.Belyschko方法是基于多重節(jié)點(diǎn)計(jì)算的經(jīng)典方法,點(diǎn)A位移與Belyschko方法[9]進(jìn)行比較.
采用2種方法計(jì)算得到的右端點(diǎn)位移曲線如圖5所示,10倍步長(zhǎng)比下位移曲線局部放大如圖6所示.
(16)
式中:n表示主時(shí)間步數(shù);Dthe為節(jié)點(diǎn)A理論位移;Dpro為采用給定方法計(jì)算的節(jié)點(diǎn)位移.統(tǒng)計(jì)的節(jié)點(diǎn)位移精度指標(biāo)見(jiàn)表1所示.
圖6 圖5的局部放大
表1 節(jié)點(diǎn)A位移計(jì)算精度指標(biāo)
圖7 多層建筑受沖擊載荷(m)
從表1中可見(jiàn):采用相同計(jì)算模型,隨著步長(zhǎng)比的增加,2種算法的最大計(jì)算誤差以及誤差和都在逐漸增加,計(jì)算精度相應(yīng)降低.在相同步長(zhǎng)比的條件下,與Belyschko方法相比,本文方法的最大誤差以及誤差和在數(shù)值上更小,計(jì)算指標(biāo)變化的幅度也較小,表明本文方法具有較高的計(jì)算精度與穩(wěn)定性.
多層鋼結(jié)構(gòu)建筑模型如圖7所示,整體建筑高 39.5 m,樓層底面長(zhǎng)度為48 m,寬為16 m,建筑共有13層,其中最低層為高 3.5 m,其他樓層為3 m.整個(gè)建筑結(jié)構(gòu)采用8節(jié)點(diǎn)六面體單元離散,離散后的該建筑模型含有 6 918 660 節(jié)點(diǎn),5 212 768 單元,20 755 980 自由度.仿真工況計(jì)算采用對(duì)結(jié)構(gòu)左側(cè)面的右三角脈沖載荷.動(dòng)力響應(yīng)分析總時(shí)間為 0.6 s,大步長(zhǎng)分區(qū)時(shí)間步長(zhǎng)Δt=200 ms,總計(jì) 30 000 個(gè)主計(jì)算時(shí)間步.
為了驗(yàn)證本文并行算法的計(jì)算效率,建筑結(jié)構(gòu)最右端(圖7中虛線框區(qū)域)接受沖擊載荷部位子分區(qū)采用較小的時(shí)間步長(zhǎng)計(jì)算,計(jì)算建筑高層點(diǎn)B的位移曲線,P為沖擊載荷,Pmax為三角沖擊載荷的最大值.將整個(gè)建筑單元分別劃分為120、240、480與960個(gè)子區(qū)域,小步長(zhǎng)分區(qū)網(wǎng)格在120子分區(qū)時(shí)占據(jù)6個(gè)分區(qū).
圖8為建筑結(jié)構(gòu)在步長(zhǎng)比為4的情況下節(jié)點(diǎn)B的位移時(shí)間曲線,該實(shí)例計(jì)算結(jié)果與LS-DYNA軟件計(jì)算結(jié)果較相符.建筑結(jié)構(gòu)節(jié)點(diǎn)位移沒(méi)有理論解,位移求解精度只做定性比較.
圖8 多層建筑節(jié)點(diǎn)B的位移曲線
異步長(zhǎng)并行算法中,分區(qū)數(shù)目與分區(qū)步長(zhǎng)比都會(huì)影響并行計(jì)算效率.表2為各個(gè)分區(qū)采用相同的時(shí)間時(shí)根據(jù)顯式異步長(zhǎng)算法在不同計(jì)算核數(shù)下得到的并行效率.從表2中可知,隨著計(jì)算核數(shù)的增加,計(jì)算總時(shí)間逐漸減少,計(jì)算加速比增加,但是從計(jì)算效率的角度來(lái)看,并行計(jì)算效率在逐漸下降.
表3為在圖8中小步長(zhǎng)分區(qū)m=4的條件下,隨著分區(qū)數(shù)目變化的顯式異步長(zhǎng)算法并行效率.表3所得并行計(jì)算效率結(jié)論與表2類(lèi)似,由兩表格數(shù)據(jù)相比較可知:異步長(zhǎng)并行條件計(jì)算下,步長(zhǎng)比的選取對(duì)于并行效率沒(méi)有實(shí)質(zhì)性提升;但是,局部區(qū)域采用較小的時(shí)間步長(zhǎng),其他區(qū)域采用較大時(shí)間步長(zhǎng),計(jì)算總時(shí)間大幅度減少.
表2 分區(qū)同步長(zhǎng)條件下計(jì)算效率(m=1)
表3 同步長(zhǎng)比條件下計(jì)算效率(m=4)
表4為在核數(shù)等于480時(shí),隨著最右端分區(qū)步長(zhǎng)比變化所得的并行計(jì)算效率.由表4可以看出,隨著步長(zhǎng)比的增加,計(jì)算總時(shí)間與計(jì)算時(shí)間比率大幅下降,計(jì)算時(shí)間比率減少幅度在逐漸下降.步長(zhǎng)比增加時(shí),大步長(zhǎng)分區(qū)計(jì)算時(shí)間降低有限,在一個(gè)異步長(zhǎng)通信周期內(nèi),分區(qū)負(fù)載不均衡現(xiàn)象顯現(xiàn),因此研究負(fù)載均衡算法可以進(jìn)一步提升顯式異步算法性能,降低計(jì)算時(shí)間比率.
表4 相同分區(qū)數(shù)不同步長(zhǎng)比下計(jì)算時(shí)間
本文提出一種基于多重節(jié)點(diǎn)的預(yù)測(cè)校正形式顯式異步結(jié)構(gòu)動(dòng)力學(xué)并行計(jì)算方法,該方法具有如下特點(diǎn):
(1) 采用顯式Newmark求解格式,對(duì)時(shí)間離散形式的結(jié)構(gòu)動(dòng)力學(xué)方程自由度解耦,邊界信息傳遞只與相鄰分區(qū)有關(guān),提高了并行計(jì)算過(guò)程信息傳遞效率.
(2) 采用節(jié)點(diǎn)分割將計(jì)算模型分割成若個(gè)子區(qū)域,各個(gè)子區(qū)域選擇各自的時(shí)間步長(zhǎng),解決了因局部精細(xì)網(wǎng)格造成的臨界時(shí)間步長(zhǎng)增加這一難題.
(3) 采用多重節(jié)點(diǎn)的分區(qū)耦合策略,相鄰分區(qū)的部分區(qū)域設(shè)定為耦合分區(qū),子循環(huán)過(guò)程中,預(yù)測(cè)波形得以完整地在整個(gè)耦合區(qū)域傳播,而不是傳統(tǒng)方法中經(jīng)過(guò)截?cái)嗟牟ㄐ?這在一定程度上提高了計(jì)算精度.
研究結(jié)果也進(jìn)一步表明,隨著步長(zhǎng)比增加,分區(qū)負(fù)載不平衡在一定程度上阻礙了計(jì)算效率的提升.因此,下一步研究的重點(diǎn)是提高異步并行負(fù)載的均衡性能.