基于重疊邊界的有限元與離散元耦合異步長并行計算方法

李 佟 金先龍 王小衛 潘浚銘 楊培中

(1.上海交通大學機械與動力工程學院, 上海 200240;2.上海交通大學機械系統與振動國家重點實驗室, 上海 200240;3.上海航天精密機械研究所, 上海 201600)

0 引言

計算機硬件技術的發展使大規模復雜力學問題的數值計算更加高效[1]。有限元法(FEM)[2]和離散元法(DEM)[3]分別是能夠獲得結構宏觀響應和細觀演化兩種不同尺度的重要數值方法。面對復雜的力學問題,結合兩種算法的優勢實現多尺度耦合計算對提升計算效率意義重大[4-6]。充分發揮計算機硬件資源優勢,結合高效的數值算法,將計算模型進行區域分解[7]開展異步長計算,針對不同區域采用各異的計算方法[8]解決大規模力學問題也是進一步提升數值模擬計算效率的一個突破點。

很多學者研究了有限元與離散元耦合多尺度計算,取得了重要進展[9-14]。在上述研究中,采用區域分解的方式在一定程度上縮短了計算時間,然而這些工作往往由于算法穩定性限制[15-16]而使整個模型須采用單一的最小積分時間步長。異步長計算方法可有效改善單一積分時間步長造成的效率低下問題,早期的研究者主要有SMOLINSKI[17]和BELYSCHKO等18],研究對象多針對兩分區結構動力學求解。MAHJOUBI等[19]引入拉格朗日乘子約束位移連續研究結構動力學問題的多分區多時間步長算法。LINDSAY等[20]采用域分解和多時間步長方法,在斷裂區域采用細化網格研究了斷裂力學問題。陳麗華等[21]研究了混合時間步長顯式積分并行算法解決沖擊動力學問題。馬志強等[8,22-23]基于節點劃分和重疊邊界研究了結構動力學的顯式多時間步長并行算法和顯隱混合異步長交錯計算方法。

目前異步長方法的研究多應用在求解有限元問題,對離散元及其耦合格式相關的異步長計算方法報道較少。MARSHALL等[24-25]在對氣溶膠的力學行為研究中,提出了離散元多時間步長計算方法,其在流體計算、顆粒粘接計算和顆粒碰撞計算中分別采用3種不同的積分時間步長。CHAUDRY等[26]研究了顆粒材料的破碎問題,通過在全部-局部框架中采用不同積分時間步長解決了由于粒子運動不穩定導致的臨界時間步長減小的問題。文獻[24,26]均利用異步長算法縮短了整體計算時間,但是對于不同步長比下方法效率并沒有做深入的研究。張銳等[27-28]基于元/網格動量傳遞的動力學過渡算法建立了有限元與離散元耦合的時間與空間多尺度計算模型,該方法在邊界耦合過渡區含有多層離散單元,耦合區自由度較多且離散單元尺寸與有限元網格尺寸有關。而且上述研究中采用的都是串行的求解格式,事實上,多核并行計算對于提升多尺度問題計算效率同樣具有重要作用[29]。

本文基于重疊邊界提出一種改進的FEM-DEM耦合異步長并行計算方法,分析不同步長比的串行和并行求解格式的計算效率。將系統模型進行區域分解,針對分區計算類型采用合適的積分時間步長[30]。將重疊邊界數據連續性傳遞引入到并行計算中,利用動態重疊邊界方法和進程間通信來傳輸相鄰子域之間的邊界信息。邊界數據傳輸過程中不涉及插值、截斷過程,在保證計算精度的同時有效提升數值計算效率。

1 FEM-DEM耦合模型建立

1.1 顯式有限元求解格式

在有限元計算中,含阻尼的線彈性結構動力學方程可以寫為[31]

(1)

式中M、C、K——質量矩陣、阻尼矩陣和剛度矩陣

F——外載荷向量

有限元求解中采用預測校正格式的Newmark積分格式,預測校正格式求解過程分為預測步、加速度求解以及校正步3個求解過程[32]。預測步中,第n+1時間步節點位移以及速度預測值可以由第n時間步節點的位移、速度以及加速度表示為

(2)

β、γ——Newmark時間積分控制參數,取β=0.25,γ=0.5

Δt——積分時間步長

此時方法具有二階精度,是無條件穩定的[33]。將式(2)代入到式(1)求出第n+1時間步節點加速度,表示為

(3)

式中Fn+1——第n+1時間步外載荷向量

(4)

1.2 顆粒離散元求解格式

離散元算法[34]中離散單元遵從牛頓第二定律,即

(5)

Fij——作用在單元i的相鄰單元j對它的作用力

Fi——作用在單元i的外載荷力

Ni——與單元i相鄰發生作用的所有單元數量

N——系統中參與計算的所有單元數量

為了使離散元法準確描述材料的連續介質力學特性,采用連接鍵[35]形式的離散元法來處理連續介質力學問題。將連續介質離散成若干個攜帶質量、位移、速度和加速度的離散單元(本文稱為顆粒),每個顆粒視為連續介質中的一個等效質點,通過相鄰顆粒之間設定的連接鍵來表示兩個顆粒之間的相互作用關系,實現連續介質的力學特性表達。

本文設定顆粒之間是彈性相互作用,通過連接兩個顆粒之間連接鍵的變形來表示顆粒的相對運動,所有顆粒的運動即為連續介質宏觀的特性。為了計算便捷,假設在整個計算過程中彈性系數等都保持不變,忽略加載歷史和顆粒變形等。

為了統一耦合模型的求解格式,采用velocity-Verlet[36]顯式積分算法來求解顆粒的運動方程,即

(6)

式中xi(n+1)——單元i第n+1個時間步位移

vi(n+1)——單元i第n+1個時間步速度

ai(n+1)——單元i第n+1個時間步加速度

fi(n+1)——單元i第n+1個時間步合力

xi(n)——單元i第n個時間步位移

vi(n)——單元i第n個時間步速度

ai(n)——單元i第n個時間步加速度

1.3 有限元離散元耦合求解格式

在有限元求解格式中,將式(1)寫成如下形式

(7)

其中

在離散元求解格式中,將式(5)轉換為

(8)

其中

比較式(7)和式(8),可以看出在兩種不同的計算方法中,運動方程在求解形式上保持一致,可以將兩種不同的方法耦合到一個計算程序中[5]。

2 交替格式的耦合異步長并行計算方法

耦合求解格式也適用于多個計算分區多核的求解過程[23],為了方便描述,以有限元和離散元兩個計算分區為例進行異步長并行計算,采用邊界重疊計算單元的形式實現一個系統時間步內有限元和離散元區域邊界數據信息交換。

2.1 基于重疊邊界的區域分割

將整個計算系統劃分為有限元和離散元兩個計算子分區,在有限元分區中,將參與計算的節點分為外部節點、邊界節點和內部節點[8]。在離散元分區中,將參與計算的顆粒分為外部顆粒、邊界顆粒和內部顆粒3類。其中有限元的邊界節點和外部節點與離散元的外部顆粒和邊界顆粒分別對應,構成重疊邊界區域,用于邊界數據傳遞,如圖1所示。在有限元異步長計算中,需要多重的外部節點來保證在同系統時間步時邊界節點數據可以完整求出[8]。相比于單純的有限元多重節點的劃分,本文提出的方法在傳遞邊界數據時只有各自的一層外部節點(顆粒)和邊界節點(顆粒),這是因為在離散元計算中,通常假定顆粒只與和它相鄰的顆粒發生接觸,即只有最近鄰顆粒會發生接觸,不考慮次近鄰和其他近鄰的顆粒。在有限元計算中,節點只與和它相鄰的節點才會對總體剛度矩陣起作用。與多重節點方法相比,算法更加簡潔高效,可操作性更強。

圖1 重疊邊界有限元-離散元分區示意圖Fig.1 FEM-DEM partitioned schematic with overlapping boundary

2.2 并行計算方法

目前并行計算已經成為模擬當中重要的手段,尤其是大規模多自由度計算,采用并行計算可以發揮多核計算機優勢提高計算效率[37-38]。由于有限元和離散元均采用顯式計算格式,顯式算法具有解耦特性(求解下一時刻的相關信息只需用到相鄰節點上一時步的相關信息),這種求解格式的解耦特性使得程序具有天然并行性,為并行化計算提供了便利[22]。通過MPI[39]進行并行計算,單個計算進程負責一個計算分區的任務處理,在單個計算區域內部同時包含該區域內部和區域邊界的計算信息,通過相鄰計算區域的信息傳遞來保障系統數據的準確性和連續性。

并行算法性能常用加速比[40]來進行衡量。考慮不同步長比時的加速比,計算式為

(9)

式中t1——同步長時(即步長比為1時)計算時間

tη——步長比為η時的計算時間

2.3 交錯格式異步長并行計算流程

一般有限元的積分時間步長條件比離散元條件寬松,因此本文離散元采用小步長,有限元采用大步長。一個完整的系統時間步包含一個有限元分區時間步和多個離散元分區時間步。計算時間步長分別為Δt和Δtr,滿足Δt=ηΔtr(本文η均為整數)。為了描述異步長并行算法的計算流程,以步長比η=2進行說明,有限元分區(分區1)和離散元分區(分區2)各自采用CPU的一個進程去進行計算。分區2完成一個時間步的計算時,分區1并沒有任何動作,當分區2又進行了一個時間步的計算后,此時分區1根據分區內部設定的計算時間步長完成一次計算,整個計算時間步也達到了一個系統時間步。圖2為步長比η=2時,分區1和分區2中子循環時間步和系統時間步的并行計算流程。非重疊形式的異步長方法,邊界位移、速度和加速度均假定為時間的線性插值函數,根據位移、速度、加速度與時間的關系,假定加速度成線性變化,通過積分速度為時間的二次函數,而位移則為更高的三次函數,在處理邊界信息時會引起數據的不一致性[32]。因此在處理邊界信息時通過MPI通信協議將兩個進程內分區1和分區2重疊邊界處的位移和速度信息進行傳遞和互換。

圖2 耦合形式的異步長并行求解流程Fig.2 Parallel computing process of coupled multi-time step

圖3 不同分區邊界信息傳遞示意圖Fig.3 Schematic of boundary information transmission in different subdomains

在有限元計算中,剛度矩陣往往是具有一定帶寬的稀疏對稱矩陣,求解加速度過程中涉及矩陣求逆。采用稀疏矩陣的行壓縮CSR存儲格式[41]有較好的存儲和求解效率。因此本文處理有限元分區計算過程采取稀疏矩陣求解的格式。而離散元求解直接根據節點自由度編號形成方程,采用坐標存儲對應節點的信息即可。

3 精度驗證

圖4 6自由度彈簧-質量系統Fig.4 Six-DOF spring-mass system

離散元分區采用的時間步長為1×10-4s,有限元分區采用的時間步長為1×10-4η。取質量塊3的位移作為參考。圖5為采用不同步長比時與解析解的對比結果。從圖5可以看出,節點位移曲線基本重疊,幾種不同的方法都可以準確地計算出質量塊3位移,表明本文開發的有限元-離散元耦合異步長算法能夠滿足計算精度。

圖5 不同方法計算出的質量塊3位移曲線Fig.5 Mass-3 displacement calculated by different methods

采用相同步長比(η=5)下串行和并行求解進行對比。質量塊3的位移誤差(不同方法減去理論值結果)如圖6所示。從圖6可以看出,有限元-離散元耦合形式的異步長算法在相同計算時間步長、步長比的條件下并行計算誤差和串行計算誤差相比沒有發生變化,反映本文提出的分區并行數據傳遞策略的有效性。

圖6 質量塊3位移誤差曲線Fig.6 Displacement error curves of mass-3

4 數值算例和有效性評價

為了進一步驗證算法在精度和并行效率方面的特性,采用超材料抵抗沖擊載荷作為數值算例。基于質量-彈簧微結構超材料模型[42],研究超材料對于沖擊載荷起到的衰減作用。超材料系統區域分解及施加的脈沖載荷如圖7所示。超材料微結構由外殼和內部的諧振子構成,質量分別為ms和mv;諧振子與外殼之間的諧振子彈簧、外殼之間彈簧剛度分別為ks和kn。諧振子、外殼及諧振子彈簧共同組成一個超材料微結構,當滿足一定條件時,等效成為負質量實體[43],可以對沖擊載荷起到衰減作用。

圖7 超材料系統區域分割及載荷歷程Fig.7 Metamaterial system domain partition and load history

超材料微結構的動力學方程為

(1≤j≤N,i=j+M)

(10)

將整個計算模型劃分為有限元(分區1)和離散元(分區2)兩個計算分區,系統仿真計算時間設定為75 ms。一個超材料微結構自由度為2個,為了使各分區參與計算的自由度數量相等,因此M=100,N=200,Q=500,計算參數設定為:ms=0.05 kg,mv=0.05 kg,kn=7.8×106N/m,ks=2.49×106N/m。離散元分區計算中的積分時間步長設定為1×10-4s。當系統中沒有超材料結構時,系統由M+N+Q個實心顆粒組成。圖8為采用不同方法獲得的離散元分區中第250號顆粒的時域波形(位移與最大位移之比)。從圖8可以看出,超材料顯著降低了結構位移響應幅值。當步長比η=5時,所提出的有限元-離散元耦合異步長并行計算方法與Newmark[44]方法獲得的位移響應一致,驗證了本文所提出方法的有效性和精度。

圖8 離散元分區中第250號顆粒的時域波形圖Fig.8 Time domain waveform of 250th particle in DEM domain

對于異步長算法而言,分區步長比、串行及并行格式都會影響計算效率。為了研究算法在不同步長比條件下的計算效率,統計采用多組步長比時消耗的時間,結果見表1。從表1可以看出,隨著步長比的增加,計算總時間逐漸減少,算法的加速比增加,步長比達到20時的計算時間縮短為步長比為1(同步長)時的65.6%。由于在計算過程中,采用異步計算模式,在未達到一個系統時間步內,分區2進行多次計算;當達到一個系統時間步內,分區1進行一次迭代計算,提高步長比意味著在分區2步長不變的情況下,分區1采用更大的計算時間步長,從而降低了分區1的迭代次數,因此分區1計算時間有所縮減,進而提高了整體計算效率。

表1 不同步長比下串行計算時間Tab.1 Computing time with different time-step ratios

圖9為相同步長比條件下采用串行格式與并行格式所消耗的求解時間。可以看出,隨著步長比的增加,采用并行格式的求解時間逐漸降低,這與串行求解格式下得到的結論一致。在相同步長比時,采用并行求解格式的計算時間低于串行求解格式的計算時間,也就是說采用異步長并行求解格式,在串行異步長算法基礎上,使得計算效率又得到了進一步的提升。從圖9也可看出,當步長比達到5以上時,繼續增加步長比,并行計算消耗的時間并沒有明顯降低,而是維持在一個穩定狀態。

圖9 不同格式下的計算時間Fig.9 Calculation time in different formats

從計算結果來看,隨著步長比的增加,無論是串行格式還是并行格式的計算時間與時間比率的下降幅度都在逐漸減小,計算時間并不會隨著步長比增加而一直保持相同的速率下降,而是趨于一個穩定值。這是因為在一個異步通信周期內,進程之間等待時間過長,出現了明顯的分區負載不均衡現象。分區計算負載不平衡在一定程度上阻礙了計算效率的提升,因此研究適用于有限元和離散元耦合的負載均衡算法,采用更多核數的并行計算可以進一步提升多時空多尺度異步長算法性能。

5 結論

(1)所有分區均采用相類似的顯式積分算法求解格式,通過數據切分和MPI通信協議充分發揮多核計算機性能優勢,邊界信息傳遞只與相鄰的分區有關,提高了并行計算過程信息傳遞效率。

(2)傳遞過程中不涉及位移和速度插值、截斷等,有效地提高了計算精度。相比于單純的有限元多重節點的耦合策略,重疊區域范圍小,數據傳遞量小,傳遞效率高。

(3)重疊區域只有兩層節點(顆粒),各分區根據自己的單元特征選擇適合的時間步長,突破了以往有限元-離散元耦合計算采用多層過渡區域的限制。