面向爆轟沖擊的分離式流固耦合數(shù)值模擬

張 森，郭曉威，*，甘新標(biāo)，龔春葉，楊文祥，李 超

(1.國(guó)防科技大學(xué)計(jì)算機(jī)學(xué)院，量子信息研究所兼高性能計(jì)算國(guó)家重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙 410073；2.中國(guó)空氣動(dòng)力研究與發(fā)展中心，綿陽 621000)

0 引言

爆轟波的沖擊響應(yīng)通常涉及高溫高壓下的流體動(dòng)力學(xué)，以及固體的變形和損傷等。這一復(fù)雜過程固有的多物理和多尺度特性導(dǎo)致了實(shí)驗(yàn)的高成本和高難度。隨著高性能計(jì)算技術(shù)和數(shù)值模擬方法的發(fā)展，對(duì)爆轟沖擊下固體響應(yīng)進(jìn)行數(shù)值模擬已成為一種有效的工程研究方法[1-3]。

爆轟傳播是一個(gè)典型的兩相或多相流問題?？紤]到不同域（如炸藥和空氣）之間跨界面的跳躍條件，該問題的主要數(shù)值難點(diǎn)是精確計(jì)算界面運(yùn)動(dòng)，可能的解決方法包括界面跟蹤法[4]、基于界面重構(gòu)過程的VOF（volume of fluid）方法[5]、水平集（level set,LS）方法[6]、擴(kuò)散界面法[7-9]等。其中擴(kuò)散界面法不需要計(jì)算界面的精確位置，而需要為過渡層中的混合物引入與擴(kuò)散界面相對(duì)應(yīng)的混合物模型或人工狀態(tài)方程。利用這一方法，Larrouturou提出了適用于理想氣體狀態(tài)方程的γ模型和質(zhì)量分?jǐn)?shù)模型[7]。然而，這些模型可能導(dǎo)致壓力或速度在界面上的振蕩。為了解決這個(gè)問題，Abgrall和Karni將γ模型的非保守方程應(yīng)用于該問題[8]。Shyue通過添加更多與狀態(tài)方程中所用參數(shù)相對(duì)應(yīng)的輸運(yùn)方程，將該模型擴(kuò)展到其他類型的狀態(tài)方程[9]。此外，離散玻爾茲曼方法（DBM）應(yīng)用在爆轟燃燒領(lǐng)域并取得進(jìn)展[10]，張玉東等基于離散玻爾茲曼模型，模擬和研究了四種不同反應(yīng)速率的爆轟現(xiàn)象，從水動(dòng)力量、非平衡量和熵產(chǎn)三方面研究了四種起爆方式的差異[11]。blastFoam[12-13]是Synthetik applicated Technologies基 于 開 源 軟 件OpenFOAM開發(fā)的適用于高爆轟、爆炸安全、氣流和一般可壓縮流動(dòng)的單相多相可壓縮流動(dòng)庫，包括13個(gè)狀態(tài)方程，允許在極端條件下對(duì)各材料進(jìn)行建模。blastFoam提供了一個(gè)比較完善的功能來模擬爆轟波的傳播過程，但是無法模擬爆轟沖擊引起的固體變形和損傷。固體的損傷變形可通過固體力學(xué)進(jìn)行建模計(jì)算，采用有限元方法對(duì)數(shù)學(xué)模型進(jìn)行離散。開源有限元框架deal.II[14]是一個(gè)支持創(chuàng)建有限元代碼的C++程序庫，該程序庫包含網(wǎng)格處理和細(xì)化、自由度處理、網(wǎng)格輸入和圖形格式的結(jié)果輸出等細(xì)節(jié)。

流固耦合（fluid structureinteraction,FSI）的主要理論模型可以從流體動(dòng)力學(xué)和固體力學(xué)的建模中得到。但是，除了單一的物理模型外，還需要考慮不同物理模型之間的相互作用。流固耦合數(shù)值模擬一般采用整體式或分離式方法。整體式方法將耦合系統(tǒng)作為整體來推導(dǎo)全局方程，然后用統(tǒng)一的方法對(duì)其進(jìn)行離散求解。該方法是魯棒且高效的，但其模塊化不足，在解決新問題時(shí)無法重用已有的代碼，需要進(jìn)行重構(gòu)、離散化、編碼。分離式方法固有的靈活性，可以重用現(xiàn)有相對(duì)成熟的單物理場(chǎng)求解器，快速開發(fā)代碼，并能在有限的時(shí)間內(nèi)快速解決問題。

分離式方法將FSI問題分解為多個(gè)域，分別求解每個(gè)子域，并通過不同域之間的接口進(jìn)行耦合。分離式FSI主要包括兩個(gè)任務(wù)，即：流體與固體方程之間的耦合格式和流體與固體在界面上的數(shù)據(jù)映射。在過去的幾十年中，擬牛頓耦合方法變得越來越重要。Michler等提出的interface-Newton-Krylov方法[15]是第一個(gè)擬牛頓耦合方法。之后，Vierendeels等在塊迭代系統(tǒng)的基礎(chǔ)上，提出了界面塊擬牛頓最小二乘（IBQN-LS）耦合方案[16]，Degroote等提出了界面擬牛頓逆最小二乘法（IQN-ILS）[17]，Bogaers等提出了多矢量擬牛頓（MVQN）格式[18]。對(duì)于數(shù)據(jù)映射方法，常用的有最近鄰法、最近投影法、徑向基函數(shù)映射法等。精確代碼交互耦合環(huán)境preCICE提供了多物理耦合所需的所有組件，包括通信方法、數(shù)據(jù)映射方法以及基于擬牛頓耦合方法的高效和健壯的耦合迭代技術(shù)[19-21]。

分離式方法具有良好的功能可擴(kuò)展性和靈活性，但其收斂效率和并行可擴(kuò)展性通常較差。De Nayer等采用基于CoMA的分離式方法模擬了帶有柔性分流板（PfS-1a）的圓柱繞流[22]，其效率和可擴(kuò)展性受順序耦合格式和集中式工具的限制，客戶端與服務(wù)器之間的通信存在瓶頸。為解決上述問題，我們基于點(diǎn)對(duì)點(diǎn)的通信布局來傳輸耦合數(shù)據(jù)，消除瓶頸。

綜上所述，本文針對(duì)爆轟沖擊下的固體響應(yīng)，基于開源軟件OpenFOAM（blastFoam[12]）、deal.Ⅱ[14]、preCICE[19]實(shí)現(xiàn)分離式流固耦合的求解系統(tǒng)，兼具靈活性和可擴(kuò)展性，同時(shí)采用具體實(shí)例，即三維豎直墻體在高爆轟作用下的運(yùn)動(dòng)過程，進(jìn)行數(shù)值模擬，并對(duì)模擬結(jié)果進(jìn)行分析和討論，以驗(yàn)證求解系統(tǒng)的正確性。

1 數(shù)值方法

分離式流固耦合數(shù)值求解部分包含三個(gè)模塊，即流體計(jì)算模塊、固體計(jì)算模塊和耦合計(jì)算信息交換模塊。各個(gè)模塊的數(shù)值方法也不盡相同。

1.1 流體計(jì)算模塊

爆轟波傳播是一種高度可壓縮的超音速流動(dòng)。基于OpenFOAM軟件框架和blastFoam庫，采用有限體積法對(duì)流體模型進(jìn)行離散求解。

作為流體動(dòng)力學(xué)的基本控制方程，Navier-Stokes（N-S）方程主要由連續(xù)方程、動(dòng)量方程和能量方程三個(gè)基本物理守恒定律組成。N-S方程可以表示為：

其中，U指由體積分?jǐn)?shù)、質(zhì)量、動(dòng)量和能量組成的矢量，F(xiàn) 指相應(yīng)的通量，S表示源項(xiàng)的向量。它們被定義為：

式中：u表 示混合物速度，I 表示單位矩陣，ρ指混合物密度，E 指總能量，p表示壓力，SM和SE是指動(dòng)量和能量源項(xiàng)， SMμ和SEμ指黏性項(xiàng)，αi和ρi表示每相的體積分?jǐn)?shù)和密度。所有體積分?jǐn)?shù)之和定義為1，混合物的密度ρ由各相的體積分?jǐn)?shù)和密度確定，表示為ρ=

壓力由一個(gè)特定的狀態(tài)方程來定義，在這個(gè)狀態(tài)方程中，混合物的內(nèi)能、密度和體積分?jǐn)?shù)被用來計(jì)算總壓力。狀態(tài)方程的Mie-Grüneisen形式使用理想氣體項(xiàng)和偏差項(xiàng)來計(jì)算壓力，即：

其中，Γi是Grüneisen系數(shù)，e是材料的內(nèi)部能量，Πi指與理想氣體的壓力偏差。同時(shí)，聲速表示為：

不同的材料對(duì)應(yīng)不同的狀態(tài)方程。對(duì)于理想氣體，滿足Γi=γi，Πi=0；對(duì)于硬化氣體，滿足Γi=γi，Πi=γiai，其中，γi表示比熱比，ai表示參考?jí)毫?。?duì)于高能炸藥，Jones Wilkins Lee（JWL）狀態(tài)方程是最常用的描述高能炸藥的狀態(tài)方程，系數(shù)滿足：

1.2 固體計(jì)算模塊

假設(shè)固體為線性彈性的，則可用彈性方程描述固體的運(yùn)動(dòng)。在deal.II軟件框架的基礎(chǔ)上，采用有限元方法對(duì)固體計(jì)算模塊進(jìn)行離散求解。

控制方程主要由運(yùn)動(dòng)方程、本構(gòu)方程和應(yīng)變位移方程組成，即：

其中，σ表示Cauchy應(yīng)力張量，f表示單位體積的體積力，ρ表示質(zhì)量密度，d 表示位移矢量，表示d對(duì)時(shí)間的二階導(dǎo)數(shù)。C表示四階剛度張量，ε表示無窮小應(yīng)變張量，運(yùn)算符“:”表示兩個(gè)張量的內(nèi)積。

在各向同性介質(zhì)中，本構(gòu)方程可簡(jiǎn)化為：

其中，λ表示Lamé第一參數(shù)，μ表示剪切模量，兩者都是Lamé參數(shù)。

在通過接口數(shù)據(jù)通信獲得應(yīng)力后，可將應(yīng)力代入本構(gòu)方程（9）求解應(yīng)變，然后代入應(yīng)變位移方程（10）求解位移。最后，將經(jīng)過一個(gè)時(shí)間步長(zhǎng)后的位移通過接口數(shù)據(jù)通信發(fā)送給流體計(jì)算模塊，完成耦合交互。

1.3 耦合計(jì)算信息交換模塊

對(duì)于分離式FSI，假設(shè)兩個(gè)求解器F 以及S，分別等價(jià)于映射和S。一個(gè)求解器將耦合接口上的輸入向量映射到輸出向量，并將其發(fā)送到另一個(gè)求解器。也就是說，以固體邊界上的位移D（ 或速度V）作為流體求解器的輸入，計(jì)算作用到固體上的應(yīng)力F作為輸出；固體求解器將這些應(yīng)力作為輸入，并計(jì)算位移（或速度）作為輸出。數(shù)學(xué)形式表示為：

分離式流固耦合數(shù)值模擬中的時(shí)間步可以在流體或固體求解器中通過少量固定數(shù)量的時(shí)間步來完成，也可以通過迭代過程直到收斂。前者對(duì)應(yīng)于顯式耦合格式，后者對(duì)應(yīng)于隱式耦合格式，兩者都有串行和并行格式。

串行顯式格式首先使用舊時(shí)間步的位移值dn作為輸入來運(yùn)行流體求解器，并獲得下一時(shí)間步的輸出應(yīng)力值 fn+1。然后使用輸出應(yīng)力值 fn+1作為輸入來運(yùn)行固體求解器以獲得下一時(shí)間步的位移值dn+1，即：

在串行顯式格式中，兩個(gè)求解器是交錯(cuò)運(yùn)行的，因此在負(fù)載平衡和最大化并行性方面不是最優(yōu)的。并行顯式格式使用舊時(shí)間步的位移和應(yīng)力值（dn,fn）作為兩個(gè)求解器的輸入，以獲得下一時(shí)間步的位移和應(yīng)力值（dn+1,fn+1）：

兩個(gè)解算器之間沒有依賴關(guān)系，因此可以在同一時(shí)間步中同時(shí)執(zhí)行。然而，這兩種方案都可能導(dǎo)致不穩(wěn)定，可以使用迭代方法迭代執(zhí)行兩個(gè)求解器，直到收斂。對(duì)于串行隱式格式，表達(dá)式為：

而對(duì)于并行隱式格式，表達(dá)式為：

當(dāng)?shù)諗繒r(shí)，則滿足：

由于求解器之間相互獨(dú)立，所使用的網(wǎng)格在界面上可能會(huì)出現(xiàn)不一致的情況，因此需要在不匹配的網(wǎng)格之間進(jìn)行數(shù)據(jù)映射。常見的數(shù)據(jù)映射方法有最近鄰映射、最近投影映射、徑向基函數(shù)（RBF）映射等。

求解器之間通信的主要任務(wù)是交換耦合數(shù)據(jù)。開源耦合庫preCICE為通信任務(wù)提供了解決方案，使用MPI或TCP/IP套接字建立通信。在耦合模擬之前，preCICE將在每個(gè)求解器進(jìn)程之間建立一個(gè)通信通道，每個(gè)并行求解器將選擇一個(gè)進(jìn)程作為“主進(jìn)程”來控制整個(gè)數(shù)值模擬過程。在通信過程中，采用了點(diǎn)對(duì)點(diǎn)的通信方式，不使用任何中央服務(wù)器。通過先前建立的通信通道，各求解器進(jìn)程之間的通信是本地和異步的[19]。

2 驗(yàn)證算例的幾何構(gòu)型與參數(shù)配置

圖1 測(cè)試算例的計(jì)算域示意圖及其三視圖Fig.1 Computational domain and itsthree views

采用三維豎直墻體在高爆轟作用下的運(yùn)動(dòng)過程作為測(cè)試算例驗(yàn)證和分析流固耦合數(shù)值模擬的求解系統(tǒng)。測(cè)試算例的計(jì)算域示意圖及其三視圖如圖1所示，計(jì)算域分為兩部分：流體部分和固體部分。固體部分為0.5a ×2.25a ×10a的豎直墻體，流體部分則是立方體空間內(nèi)去除豎直墻體剩下的區(qū)域。紅點(diǎn)表示爆轟中心，距固體和下邊界的距離均為0.3048 m。在中心截面處豎直墻體周圍設(shè)置了六個(gè)采樣點(diǎn)，采樣點(diǎn)的分布情況如圖2所示，紅點(diǎn)表示爆轟中心，青點(diǎn)表示各采樣點(diǎn)。使用兩種不同規(guī)模的網(wǎng)格（分別為457.6萬和764萬個(gè)單元）進(jìn)行網(wǎng)格無關(guān)性檢驗(yàn)，結(jié)果如圖3所示。由圖中可以看出采樣點(diǎn)1和2處的壓力變化與網(wǎng)格規(guī)模無關(guān)，因此數(shù)值模擬的結(jié)果與網(wǎng)格規(guī)模無關(guān)。下文中除并行測(cè)試外的數(shù)值模擬結(jié)果均采用以下結(jié)構(gòu)化網(wǎng)格，即流體網(wǎng)格包含466.9萬個(gè)點(diǎn)，1382萬個(gè)面和457.6萬個(gè)單元。爆轟中心附近和流體與固體界面處的網(wǎng)格較密集，遠(yuǎn)離爆轟中心處網(wǎng)格較稀疏。固體網(wǎng)格包含64.8萬個(gè)點(diǎn)、59萬個(gè)面和19.2萬個(gè)單元，自由度為618723。

圖2 中心截面處豎直墻體周圍的采樣點(diǎn)分布示意圖Fig. 2 Distribution of sampling pointsaround the central section of the vertical wall

圖3 網(wǎng)格無關(guān)性檢驗(yàn)(采樣點(diǎn)1和2處的壓力變化)Fig.3 Mesh independency test(pressure at sampling points1 and 2)

流體部分由兩相材料組成，即C4炸藥和空氣。各相材料的狀態(tài)方程類型及其系數(shù)值如表1所示。C4炸藥的狀態(tài)方程類型為JWL，而空氣的狀態(tài)方程類型為硬化氣體（Stiffened Gas）。對(duì)于固體部分，豎直墻體的材料由三個(gè)參數(shù)決定，即密度（ρ）和兩個(gè)Lamé參數(shù)（λ 、μ）?；炷翂w的密度為2440 kg/m3， λ為9.03×109，μ為13.54×109；鐵板墻體密度為7000 kg/m3，λ為102.42×109，μ為80.47×109。

表1 流體計(jì)算模塊各相材料的狀態(tài)方程及其系數(shù)值Table1 The EOSand coefficients of each phase in the fluid module

數(shù)值格式是數(shù)值模擬的重要部分。對(duì)于流體求解器，通量格式選擇HLLC格式；一階時(shí)間導(dǎo)數(shù)項(xiàng)格式選擇二階強(qiáng)穩(wěn)定Runge-Kutta（RK2SSP）格式；從單元中心到面中心的插值格式為cubic插值，標(biāo)量的重構(gòu)格式是vanAlbada，向量的重構(gòu)格式是vanAlbadaV。除此之外，梯度格式為faceMDLimited leastSquares 1.0、表面法向梯度格式為corrected、拉普拉斯格式為Gauss linear corrected。對(duì)于固體求解器，選擇具有二階精度的Crank-Nicolson格式。

對(duì)于耦合計(jì)算信息交換模塊，耦合交互的數(shù)據(jù)為位移和應(yīng)力，數(shù)據(jù)映射方法采用最近鄰映射。耦合格式采用串行隱式格式（方程（15）），加速格式采用IQN-ILS[17]方法。

3 結(jié)果與討論

算例是在高性能計(jì)算平臺(tái)上進(jìn)行測(cè)試的，該平臺(tái)包含數(shù)百個(gè)計(jì)算節(jié)點(diǎn)，每個(gè)節(jié)點(diǎn)包含兩個(gè)Xeon E5-2620處理器（6核，2.10 GHz）。

圖4是Beyer報(bào)告中的爆轟過程定性分析圖，顯示出高能炸藥爆炸后的爆轟波傳播過程。在高爆轟作用下，豎直墻體受到較大的沖擊力，導(dǎo)致壁面變形。圖5是三維豎直墻體在高爆轟作用下的運(yùn)動(dòng)過程的整體模擬圖。該圖包括豎直墻體的位移、流場(chǎng)中的速度流線和密度等值線以及下邊界的速度分布。圖中豎直墻體的材質(zhì)為鐵板，鐵板墻體的變形程度小于混凝土墻體。

圖4 Beyer報(bào)告[23]中的爆轟過程定性分析圖Fig.4 A sketch of denotation process provided by Ref.[23]

炸藥在爆炸點(diǎn)起爆后很快就擴(kuò)散到豎直墻體。在0.04 ms左右，爆轟波首先與豎直墻面接觸，但豎直墻面暫時(shí)沒有受到?jīng)_擊力的影響，尚未發(fā)生變形。在0.07 ms左右，豎直墻體所受到的壓力高達(dá)10.20 GPa，此時(shí)豎直墻體在爆轟接觸點(diǎn)處有輕微變形。從圖5可以看出，隨著時(shí)間的推移，豎直墻受到持續(xù)的應(yīng)力，變形繼續(xù)從初始接觸點(diǎn)向外延伸。當(dāng)爆轟波傳播到豎直墻體下邊界的角落處時(shí)，壓力和密度不斷積累。在0.2 ms時(shí)，墻角處的壓力為最大壓力，達(dá)到52.24 MPa，密度達(dá)到47.71 kg/m3。在0.6 ms時(shí)，豎直墻體的變形達(dá)到最大值，x方向偏移達(dá)到0.045 m。當(dāng)爆轟波撞擊豎直墻體及其下邊界時(shí)，爆轟波反彈并向后傳播。爆轟波向外傳播，對(duì)豎直墻體施加力的作用，引起豎直墻體的變形。同時(shí)，由于豎直墻體的變形，爆轟波傳播的計(jì)算域發(fā)生變化，從而影響爆轟波的傳播。它們相互耦合、相互作用。爆轟波向外傳播，與豎墻相互耦合，豎直墻體右上角出現(xiàn)漩渦。數(shù)值模擬結(jié)果展示的爆轟過程與Beyer[23]報(bào)告的爆轟傳播過程相一致。

圖6顯示了豎直墻體在高爆轟作用下的形變過程。豎直墻體的材料為混凝土。在初始爆轟接觸點(diǎn)處，首先形成凹面，然后變形向外擴(kuò)展。當(dāng)時(shí)間為0.71 ms時(shí)，垂直墻的位移達(dá)到最大值，x方向的最大偏移達(dá)到0.157 m。隨著時(shí)間的推移，雖然墻體中間位置的x方向偏移逐漸減小，但墻體的變形范圍仍在不斷擴(kuò)大。

圖7、圖8和圖9分別顯示了0.4、0.8、1.2、1.6 ms的速度流線圖、壓力等值線圖和密度等值線圖。此時(shí)，爆轟波經(jīng)過下邊界和豎直墻體的反彈后繼續(xù)向外傳播。0.3 ms后，豎直墻體右上方逐漸形成渦旋，流線規(guī)則地向外輻射，壓力和密度在波面處發(fā)生躍變（包括初始波和反彈波）。當(dāng)爆轟波向外傳播時(shí)，渦逐漸變大并逐步發(fā)展為梭形渦。圖7顯示了中心截面后半部分的速度流線圖，可以看出，在0.4 ms時(shí)，整個(gè)速度流線向外輻射，在豎直墻體的右上角開始出現(xiàn)一個(gè)小渦旋。

圖5 三維豎直鐵板墻體在高爆轟作用下運(yùn)動(dòng)過程的整體模擬圖（時(shí)間跨度0.1～0.9 ms，間隔0.1 ms）Fig.5 The movement of a 3D vertical iron wall under the condition of a high-explosive detonation(the time span is0.1～ 0.9 ms,the interval time is0.1 ms)

圖6 三維豎直混凝土墻體在高爆轟作用下的形變過程圖（時(shí)間跨度0.1～1.0 ms，間隔0.1 ms）Fig.6 The movement of a 3D vertical concretewall under the condition of a high-explosive detonation(the time span is0.1～ 1.0 ms,the interval time is0.1 ms)

圖7 中截面后半部分在0.4、0.8 、1.2、1.6 ms的速度流線圖Fig.7 The velocity streamlinesin the second part of the central cross section at 0.4, 0.8 ,1.2,and 1.6 ms

圖8 中截面后半部分在0.4、0.8 、1.2、1.6 ms的壓力等值線圖Fig.8 The pressure contours in the second part of the central crosssection at 0.4, 0.8, 1.2, and 1.6 ms

在0.8 ms時(shí)，爆轟波向外傳播，速度流線仍呈現(xiàn)向外輻射狀。此時(shí)渦旋已經(jīng)非常明顯，同時(shí)呈現(xiàn)出梭形。在1.2 ms時(shí)，渦旋持續(xù)變大，流線出現(xiàn)混亂，不是規(guī)則的向外輻射狀。在1.6 ms時(shí)，流線的混亂程度變大。圖8和圖9顯示了中心橫截面后半部分的壓力等值線和密度等值線?？梢钥闯?，波面處的壓力和密度發(fā)生躍變，波面處的壓力和密度較高，而爆轟中心附近出現(xiàn)低壓和低密度區(qū)。

各采樣點(diǎn)處的壓力-時(shí)間和密度-時(shí)間歷程曲線如圖10和圖11所示。壓力和密度在爆轟波傳播到達(dá)的時(shí)候瞬間跳躍式增長(zhǎng)，爆轟波傳播過后則迅速下降。通過壓力和密度變化的延遲情況可以判斷出爆轟波到達(dá)的順序。壓力和密度出現(xiàn)了多次增長(zhǎng)，表示初始波和反彈波經(jīng)過了采樣點(diǎn)。從圖中可以看出，爆轟波最先到達(dá)1號(hào)采樣點(diǎn)，壓力最高達(dá)1.61296×107Pa，密度最高達(dá)到25.9997 kg/m3。隨后到達(dá)2號(hào)、3號(hào)采樣點(diǎn)。其他采樣點(diǎn)處的密度最高到達(dá)5 kg/m3附近，之后不斷下降趨近于零。

圖9 中截面后半部分在0.4、0.8 、1.2、1.6 ms的密度等值線圖Fig.9 The density contoursin the second part of the central cross section at 0.4, 0.8 ,1.2,and 1.6 ms

圖10 采樣點(diǎn)處的壓力-時(shí)間歷程曲線Fig.10 The time historiesof pressure

因?yàn)榛邳c(diǎn)對(duì)點(diǎn)通信模式和并行耦合格式，所以求解系統(tǒng)具有良好的并行可擴(kuò)展性。我們?cè)O(shè)置了兩種不同規(guī)模的網(wǎng)格（網(wǎng)格1：流體300萬單元，固體10萬單元；網(wǎng)格2：流體500萬單元，固體10萬單元），分別采用串行顯式（方程（13））和并行顯式（方程（14））耦合格式展開測(cè)試。圖12顯示了網(wǎng)格1和網(wǎng)格2在串行顯式和并行顯式耦合格式下的時(shí)間開銷（400個(gè)時(shí)間步）和對(duì)應(yīng)的加速比。對(duì)于同一網(wǎng)格，并行耦合格式的時(shí)間開銷相比于串行格式減少大約36%（網(wǎng)格1）和28%（網(wǎng)格2）。在256核四個(gè)測(cè)試的加速比達(dá)到最大，分別為122.3（網(wǎng)格1-串行耦合），141.2（網(wǎng)格1-并行耦合），147.0（網(wǎng)格2-串行耦合），178.0（網(wǎng)格1-并行耦合）。

圖11 采樣點(diǎn)處的密度-時(shí)間歷程曲線Fig.11 The time historiesof density

數(shù)值模擬給出了高能炸藥在三維豎直墻體后起爆并傳播的過程中各物理量的分布情況和高爆轟作用下豎直墻體的變形運(yùn)動(dòng)情況。模擬結(jié)果符合預(yù)期及爆轟過程的直觀印象，可為相關(guān)工程應(yīng)用的數(shù)字化設(shè)計(jì)和評(píng)價(jià)提供技術(shù)支持。

圖12 網(wǎng)格1和網(wǎng)格2在串行顯式和并行顯式耦合格式下的時(shí)間開銷（400個(gè)時(shí)間步）和對(duì)應(yīng)的加速比Fig. 12 Thetime costs (400 time-steps) and the speedup ratios of Mesh1 and Mesh2 with serial-explicit and parallel-explicit coupling schemes

4 結(jié)論

爆轟沖擊作用下固體響應(yīng)過程的研究在諸多工程應(yīng)用中，如損傷評(píng)估、建筑爆破、礦物開采等，具有重要意義。然而，在復(fù)雜的工程應(yīng)用中，準(zhǔn)確高效地模擬爆轟沖擊作用下固體響應(yīng)的過程存在許多困難和挑戰(zhàn)。

本文采用分離式流固耦合的方法，基于開源軟件OpenFOAM（blastFOAM）、preCICE和deal.Ⅱ?qū)崿F(xiàn)數(shù)值模擬的求解系統(tǒng)，并利用三維豎直墻壁在高爆轟作用下的運(yùn)動(dòng)算例開展數(shù)值模擬。對(duì)于數(shù)值模擬結(jié)果，我們從爆轟響應(yīng)整體過程、豎直墻體的形變過程、速度流線、壓力和密度等值線、采樣點(diǎn)處的壓力和密度變化等方面詳細(xì)展開分析。結(jié)果表明，該求解系統(tǒng)能夠準(zhǔn)確有效地模擬爆轟波傳播的各個(gè)階段以及高爆轟作用下豎直墻體的變形運(yùn)動(dòng)情況，模擬結(jié)果展示的爆轟過程與Beyer報(bào)告中的爆轟波傳播過程一致。求解系統(tǒng)具有良好的并行可擴(kuò)展性，在并行度為256核時(shí)的加速比達(dá)到141.2和178.0，并行效率分別達(dá)到55.2%和69.5%?？偠灾疚耐ㄟ^集成開放源代碼，實(shí)現(xiàn)了面向爆轟沖擊的分離式流固耦合數(shù)值模擬求解系統(tǒng)，該系統(tǒng)在諸多領(lǐng)域的工程應(yīng)用中都具有重要的現(xiàn)實(shí)意義。