馬志強(qiáng)， 樓云鋒， 李俊杰， 金先龍

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

大型復(fù)雜結(jié)構(gòu)的有限元分析常涉及局部沖擊載荷、接觸等問題，需要針對(duì)特定分析部位劃分精細(xì)化的網(wǎng)格[1].這類問題常采用顯式算法求解，但受限于求解精度以及局部單元穩(wěn)定性的要求，整個(gè)模型需要采用較小的時(shí)間步長，這會(huì)導(dǎo)致計(jì)算時(shí)間大大加長[2].因此國內(nèi)外眾多學(xué)者提出改進(jìn)方案，期望可根據(jù)結(jié)構(gòu)區(qū)域的特點(diǎn)，將結(jié)構(gòu)分成不同屬性的區(qū)域，針對(duì)單一區(qū)域選擇合適的積分算法和時(shí)間步長進(jìn)行動(dòng)力學(xué)分析.

改進(jìn)方法主要包括域分解法以及子循環(huán)方法兩大類.域分解法是將計(jì)算區(qū)域分解為若干小的子區(qū)域，各個(gè)子區(qū)域單獨(dú)計(jì)算.根據(jù)分區(qū)邊界是否重疊可以將域分解法分為重疊與非重疊兩種形式.重疊形式主要包括Schwarz迭代法和Arlequin方法[3]，非重疊形式包括Schur補(bǔ)方法和dual Schur方法，其區(qū)別在于處理區(qū)域邊界約束的方式.張友良等[4]采用Schur補(bǔ)方法開發(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ū)步長不匹配的問題.早期研究者主要有Smolinski[8]與Belyschko[9]，研究對(duì)象多針對(duì)兩分區(qū)的結(jié)構(gòu)動(dòng)力學(xué)問題.Prakash等[10]研究了非線性結(jié)構(gòu)的異步長計(jì)算，通過拉格朗日乘子約束分區(qū)邊界速度連續(xù).劉占生等[11]將子循環(huán)與預(yù)測(cè)校正方法相結(jié)合求解壁板動(dòng)力響應(yīng).高雙武等[12]采用子循環(huán)的緊耦合算法，分析固體火箭發(fā)動(dòng)機(jī)點(diǎn)火瞬態(tài)工作過程.

在上述研究中，域分解法可在一定程度上提高加速比，縮短計(jì)算時(shí)間，但不同計(jì)算區(qū)域采用相同的時(shí)間步長，在處理存在局部細(xì)分網(wǎng)格的動(dòng)力學(xué)問題時(shí)，會(huì)造成計(jì)算資源的浪費(fèi).在經(jīng)典的顯式子循環(huán)方法中，分區(qū)多采用單元分割，不同步長邊界節(jié)點(diǎn)信息傳遞過程常涉及插值或截?cái)郲13]，這在一定程度上造成計(jì)算不穩(wěn)定、降低計(jì)算精度，也會(huì)限制最大步長比的選取[14].

為此，本文提出一種基于多重節(jié)點(diǎn)的顯式異步長并行計(jì)算方法.該方法將多重節(jié)點(diǎn)引入異步長并行計(jì)算，顯式小步長分區(qū)子循環(huán)過程中預(yù)測(cè)波形得以在顯式重疊區(qū)域內(nèi)完整傳遞，沒有插值、截?cái)嗾`差，計(jì)算精度較高，不同網(wǎng)格區(qū)域可以采用較大的步長比.顯式算法具有解耦特性(求解下一時(shí)刻的相關(guān)信息只需用到相鄰節(jié)點(diǎn)上一時(shí)步的相關(guān)信息)，使并行計(jì)算過程不需要聯(lián)立求解界面方程.因此，該方法具有一般顯式并行算法易于編程實(shí)現(xiàn)的特點(diǎn)，同時(shí)相對(duì)于經(jīng)典的子循環(huán)算法精度較高.

1 預(yù)測(cè)校正形式的顯式Newmark求解格式

含阻尼的線彈性體結(jié)構(gòu)動(dòng)力學(xué)方程可以寫為

(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)

2 多重節(jié)點(diǎn)的顯式異步長并行求解方法

2.1 異步長計(jì)算變量定義

結(jié)構(gòu)有限元網(wǎng)格采用節(jié)點(diǎn)分割方法，將節(jié)點(diǎn)劃分給不同分區(qū)，分割單元被相鄰分區(qū)共享，并行計(jì)算過程需要共享分割單元信息.定義分區(qū)步長比：

(6)

式中：Δti為第i個(gè)分區(qū)選用的時(shí)間步長；Δtr為所有分區(qū)最小時(shí)間步長.異步長計(jì)算過程中涉及主時(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ù)目.

2.2 分區(qū)異步長顯式并行求解

顯式異步長求解過程需要重疊多層邊界節(jié)點(diǎn).圖1是以2倍步長比的網(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 異步長計(jì)算多重節(jié)點(diǎn)劃分示意

參照分區(qū)節(jié)點(diǎn)分割向量的定義式，小步長分區(qū)分割向量為

Oi=

(9)

式中：OI代表內(nèi)部節(jié)點(diǎn)分割向量(包含邊界節(jié)點(diǎn))；OE為外部節(jié)點(diǎn)分割向量.對(duì)步長比為mi的分區(qū)(為簡化公式暫不包含阻尼)，離散的結(jié)構(gòu)動(dòng)力學(xué)方程可以寫為

(10)

對(duì)于小步長分區(qū)的一個(gè)子循環(huán)時(shí)間步j(luò)，內(nèi)部節(jié)點(diǎn)與邊界節(jié)點(diǎn)的加速度可以寫為

(11)

(12)

對(duì)于大步長分區(qū)，內(nèi)部與邊界節(jié)點(diǎn)的加速度求解公式與式(10)和(11)相似，不同的是此時(shí)分區(qū)時(shí)間步長為Δti=miΔtr，預(yù)測(cè)步與校正步中求解公式需要將時(shí)間步長換成miΔtr.大步長分區(qū)質(zhì)量矩陣與剛度矩陣需按照分區(qū)分割向量分塊.

3 數(shù)值算例分析

3.1 并行計(jì)算性能參數(shù)

算法的有效性在國家超級(jí)計(jì)算廣州中心的“天河二號(hào)”超級(jí)計(jì)算機(jī)上進(jìn)行測(cè)試，文中算例均采用CPU核心計(jì)算.并行算法的性能常用加速比與并行計(jì)算效率來進(jìn)行衡量.異步長并行算法中，結(jié)構(gòu)分區(qū)數(shù)以及步長比不同都會(huì)導(dǎo)致計(jì)算加速比不同.定義結(jié)構(gòu)分區(qū)數(shù)目k固定，步長比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 彈性桿受軸向集中載荷

3.2 彈性桿承受集中軸向荷載

圖2為彈性桿在軸向集中荷載作用下結(jié)構(gòu)示意圖.桿左端固定，右端點(diǎn)A受軸向集中載荷.桿材料密度ρ為 1 000 kg/m3，彈性模量1 GPa，橫截面積 0.001 m2，桿的長度l為1 m，自由端加載外力fext為 1 000 N.x軸代表?xiàng)U軸向坐標(biāo).整個(gè)結(jié)構(gòu)劃分為4個(gè)子分區(qū)，每個(gè)子分區(qū)長度為 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í)間步長，步長Δt為 2.5 μs.子分區(qū)4依次使用相同時(shí)間步長Δtr(Δt=Δtr)、5倍時(shí)間步長(Δt=5Δtr)與10倍時(shí)間步長(Δt=10Δtr)共3種情況.

圖3為子分區(qū)4采用不同時(shí)間步長比，桿右端點(diǎn)A軸向位移(D)曲線與理論位移曲線.圖4為圖3的局部放大圖.圖3和4表明， 隨著分區(qū)時(shí)間步長比的增加，位移的計(jì)算精度沒有出現(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仍采用步長Δt為 2.5 μs，子分區(qū)4采用步長比為10的異步長方法，計(jì)算桿右端點(diǎn)A的軸向位移.Belyschko方法是基于多重節(jié)點(diǎn)計(jì)算的經(jīng)典方法，點(diǎn)A位移與Belyschko方法[9]進(jìn)行比較.

采用2種方法計(jì)算得到的右端點(diǎn)位移曲線如圖5所示，10倍步長比下位移曲線局部放大如圖6所示.

(16)

式中：n表示主時(shí)間步數(shù)；Dthe為節(jié)點(diǎn)A理論位移；Dpro為采用給定方法計(jì)算的節(jié)點(diǎn)位移.統(tǒng)計(jì)的節(jié)點(diǎn)位移精度指標(biāo)見表1所示.

圖6 圖5的局部放大

表1 節(jié)點(diǎn)A位移計(jì)算精度指標(biāo)

圖7 多層建筑受沖擊載荷(m)

從表1中可見：采用相同計(jì)算模型，隨著步長比的增加，2種算法的最大計(jì)算誤差以及誤差和都在逐漸增加，計(jì)算精度相應(yīng)降低.在相同步長比的條件下，與Belyschko方法相比，本文方法的最大誤差以及誤差和在數(shù)值上更小，計(jì)算指標(biāo)變化的幅度也較小，表明本文方法具有較高的計(jì)算精度與穩(wěn)定性.

3.3 多層建筑受沖擊載荷

多層鋼結(jié)構(gòu)建筑模型如圖7所示，整體建筑高 39.5 m，樓層底面長度為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，大步長分區(qū)時(shí)間步長Δt=200 ms，總計(jì) 30 000 個(gè)主計(jì)算時(shí)間步.

為了驗(yàn)證本文并行算法的計(jì)算效率，建筑結(jié)構(gòu)最右端(圖7中虛線框區(qū)域)接受沖擊載荷部位子分區(qū)采用較小的時(shí)間步長計(jì)算，計(jì)算建筑高層點(diǎn)B的位移曲線，P為沖擊載荷，Pmax為三角沖擊載荷的最大值.將整個(gè)建筑單元分別劃分為120、240、480與960個(gè)子區(qū)域，小步長分區(qū)網(wǎng)格在120子分區(qū)時(shí)占據(jù)6個(gè)分區(qū).

圖8為建筑結(jié)構(gòu)在步長比為4的情況下節(jié)點(diǎn)B的位移時(shí)間曲線，該實(shí)例計(jì)算結(jié)果與LS-DYNA軟件計(jì)算結(jié)果較相符.建筑結(jié)構(gòu)節(jié)點(diǎn)位移沒有理論解，位移求解精度只做定性比較.

圖8 多層建筑節(jié)點(diǎn)B的位移曲線

異步長并行算法中，分區(qū)數(shù)目與分區(qū)步長比都會(huì)影響并行計(jì)算效率.表2為各個(gè)分區(qū)采用相同的時(shí)間時(shí)根據(jù)顯式異步長算法在不同計(jì)算核數(shù)下得到的并行效率.從表2中可知，隨著計(jì)算核數(shù)的增加，計(jì)算總時(shí)間逐漸減少，計(jì)算加速比增加，但是從計(jì)算效率的角度來看，并行計(jì)算效率在逐漸下降.

表3為在圖8中小步長分區(qū)m=4的條件下，隨著分區(qū)數(shù)目變化的顯式異步長算法并行效率.表3所得并行計(jì)算效率結(jié)論與表2類似，由兩表格數(shù)據(jù)相比較可知：異步長并行條件計(jì)算下，步長比的選取對(duì)于并行效率沒有實(shí)質(zhì)性提升；但是，局部區(qū)域采用較小的時(shí)間步長，其他區(qū)域采用較大時(shí)間步長，計(jì)算總時(shí)間大幅度減少.

表2 分區(qū)同步長條件下計(jì)算效率(m=1)

表3 同步長比條件下計(jì)算效率(m=4)

表4為在核數(shù)等于480時(shí)，隨著最右端分區(qū)步長比變化所得的并行計(jì)算效率.由表4可以看出，隨著步長比的增加，計(jì)算總時(shí)間與計(jì)算時(shí)間比率大幅下降，計(jì)算時(shí)間比率減少幅度在逐漸下降.步長比增加時(shí)，大步長分區(qū)計(jì)算時(shí)間降低有限，在一個(gè)異步長通信周期內(nèi)，分區(qū)負(fù)載不均衡現(xiàn)象顯現(xiàn)，因此研究負(fù)載均衡算法可以進(jìn)一步提升顯式異步算法性能，降低計(jì)算時(shí)間比率.

表4 相同分區(qū)數(shù)不同步長比下計(jì)算時(shí)間

4 結(jié)語

本文提出一種基于多重節(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ì)算過程信息傳遞效率.

(2) 采用節(jié)點(diǎn)分割將計(jì)算模型分割成若個(gè)子區(qū)域，各個(gè)子區(qū)域選擇各自的時(shí)間步長，解決了因局部精細(xì)網(wǎng)格造成的臨界時(shí)間步長增加這一難題.

(3) 采用多重節(jié)點(diǎn)的分區(qū)耦合策略，相鄰分區(qū)的部分區(qū)域設(shè)定為耦合分區(qū)，子循環(huán)過程中，預(yù)測(cè)波形得以完整地在整個(gè)耦合區(qū)域傳播，而不是傳統(tǒng)方法中經(jīng)過截?cái)嗟牟ㄐ?這在一定程度上提高了計(jì)算精度.

研究結(jié)果也進(jìn)一步表明，隨著步長比增加，分區(qū)負(fù)載不平衡在一定程度上阻礙了計(jì)算效率的提升.因此，下一步研究的重點(diǎn)是提高異步并行負(fù)載的均衡性能.