歐拉坐標系下具有銳利相界面的可壓縮多介質流動數值方法研究

姚成寶 付梅艷 韓 峰 閆 凱

(西北核技術研究院，西安 710024)

(北京大學數學科學學院，北京 100871)

引言

多介質大變形問題廣泛存在于軍事和民用領域，例如爆炸力學效應、高速沖擊問題、武器設計以及自然災害的防護等.上述問題通常存在多種介質間強烈的相互作用，具有強耦合、大變形、高度非線性等特點，物理現象極其復雜.由于各介質之間的物理特性、初始狀態差異極大，很多對單介質問題比較有效的數值方法在處理多介質問題時，通常會引起較大的數值震蕩和誤差[1-2]，甚至導致計算無法繼續，因此需要研究合適的多介質數值計算方法，準確描述不同介質間的相互作用.

多介質流動數值方法按其采用的坐標系可以分為Lagrange 方法[3]、Euler 方法[4]、ALE 方法[5-6]以及Euler-Lagrange 耦合方法[7]等.由于Euler 方法的坐標系固定在空間中，計算網格不隨介質的運動而運動，不發生網格畸變，因此適合處理大變形問題.但由于介質在網格間流動，如果計算區域內存在多種介質，則有可能出現同一網格中包含多種介質的情況.如何描述相界面的位置以及混合網格中各介質之間的相互作用，是Euler 方法求解多介質問題的難點.目前比較流行的、用于描述和追蹤自由面和相界面運動軌跡的方法包括：VOF(volume of fluid)方法[8-9]、MOF(moment of fluid)方法[10]、Front Tracking 方法[11-12]以及Level Set 方法[13-14]等.此外，當網格中包含多種介質時，如果不同介質的密度和狀態方程差異很大，在相界面附近容易產生速度和壓力的非物理震蕩，因此在Euler坐標系下求解多介質問題時，應盡可能避免直接利用相界面兩側介質的狀態參數來求解控制方程，而是需要對混合網格進行特殊處理，盡可能抑制和消除非物理震蕩.目前主要的處理方法包括兩類.第一類方法是虛擬流體類方法(包括GFM[15]、MGFM[16-20]和RGFM[21-22]等)，其假設在相界面的另一側存在相同介質的虛擬流體，介質穿過相界面時滿足等熵關系，且壓力和法向速度連續.第二類方法是切割單元法(cut cell method)，該方法根據求解精度顯式地重構出相界面，并將含有多種介質的單元切割成幾個部分.切割單元法的優點是能夠清晰的描述相界面的位置和拓撲結構，而且能處理比較復雜的相界面以及拓撲結構發生變化時的情形，但缺點是單元會被相界面切割成幾個部分，人為地形成碎片單元.當碎片單元的尺寸過小時，會嚴重影響數值計算的穩定性和效率，此時需要采用額外的技術來克服上述困難.例如，Nourgaliev 等[23]采用了兩套網格來存儲數據，其中一套網格為非結構網格，另一套為結構網格.數值計算時，相界面對結構網格進行切割，并通過對切割過程中產生的碎片單元進行網格融合，得到另一套非結構網格，在一定程度上克服了CFL 條件對時間步長的限制.Liu 等[24]引入了一個混合方法，通過對碎片單元相鄰網格的守恒量進行再分配，來達到數值穩定的目的.白曉等[25]在切割單元過程中引入了虛擬流體方法，每種介質切割單元的狀態可以覆蓋到一個完整的單元，實現了碎片單元的狀態量更新，并克服了數值不穩定性.Hu 等[26-27]在二維交錯型結構網格中提出了一種守恒型的混合技術，對體積分數小于0.5 的碎片單元，根據界面的法向確定碎片單元在x，y方向以及對角線方向的同介質相鄰單元，并以界面的法向分量和各相鄰單元的體積分數作為權重來計算這3 個單元與目標碎片單元之間的守恒量交換，從而對碎片單元的守恒量進行了修正，提高了數值計算的穩定性.Meyer 等[28]在Hu 的研究基礎上進行了改進，并將算法拓展了三維情形.在對碎片單元的守恒量進行修正時，碎片單元的所有相鄰單元都被用于進行混合操作，并將體積分數的冪函數形式加入到權重因子中，使得體積分數越大的的單元所占的權重越大，進一步提高了混合技術的穩定性.Lauer 等[29]將Meyer 的混合技術擴展到了非均勻的結構網格中，利用混合技術對體積分數小于0.6 的碎片單元進行了滿足守恒性質的修正處理.

本文提出了一種基于Euler 坐標系的非結構網格、具有銳利界面的二維和三維守恒型多介質流動數值方法，可用于模擬可壓縮流體和彈塑性固體在極端物理條件下的大變形動力學行為.利用分片線性的Level Set 函數清晰重構出分段線性的相界面，建立了單純形網格內具有多種介質的相界面拓撲結構，理論上可以處理全局任意種介質、局部3 種介質的多介質流動問題.提出了一種守恒型的相界面兩側介質相互作用的處理方法，能夠保證不同介質間的通量守恒.提出了一種非結構網格上的單元聚合算法，通過對每個碎片單元建立至少包含一個完整單元的同介質單元片，并對該單元片內的守恒量進行基于單元體積的加權平均，擴大了碎片單元的特征尺度，消除了由于網格被相界面分割成較小碎片、違反CFL 條件，進而可能帶來數值不穩定的問題.最后對激波-氣泡相互作用、淺埋爆炸、空中強爆炸沖擊波和典型坑道內沖擊波傳播問題開展了數值模擬，驗證了數值方法的可靠性和計算效率.

1 物理模型

1.1 控制方程

考慮二、三維計算區域上的多介質流動問題，基于Euler 坐標系描述可壓縮流體和彈塑性固體的大變形動力學行為，其控制方程組如式(1)所示

其中，ρ 表示密度，u表示速度矢量，E表示單位體積總能量，p表示壓力，S表示固體的偏應力張量，對于流體S設為0，以下同.

1.2 狀態方程

除了描述物質(流體或固體)的質量、動量和能量守恒的控制方程外，物質動力學行為的完整描述還需要提供描述物質自身熱力學狀態的本構方程(狀態方程).本文涉及到的狀態方程包括JWL 狀態方程、理想氣體狀態方程以及剛性氣體狀態方程等.其中，JWL 狀態方程用于描述TNT 爆轟產物，如式(2)所示

其中，ρ0表示初始密度，e表示質量比內能，且滿足

A1,B1,R1,R2,為JWL 狀態方程的參數，其值分別為3.712,0.0323,4.15,0.95 和0.3.

空氣采用理想氣體狀態方程，如式(3)所示

其中，γ 表示空氣的絕熱指數.

剛性氣體狀態方程(4)可用于描述典型的固體等在受到爆炸或沖擊時的動態響應問題

和流體不同的是，固體材料除了能夠承受壓縮或膨脹帶來的體積變化外，還能夠承受一定程度的剪切變形(形狀改變).流體彈塑性模型將固體的力學響應過程分解為兩個獨立部分：體積變形和剪切變形，其中體積變形仍然由狀態方程來確定，而剪切變形在受到強沖擊或者爆炸載荷作用時可能經歷3 個過程：彈性變形階段、塑性變形階段以及流體動力學階段，如式(5)所示

其中，sij表示偏應力分量，表示Von Mises 等效應力，表示偏應變率，μe和μp分別表示彈性和塑形剪切模量，Ye和Yp分別表示彈性和塑形屈服應力.

2 多介質流動數值方法

采用基于歐拉坐標系、具有互不相溶特征的二、三維多介質流動計算模型來進行數值離散.在任意時刻，整個計算區域? 可以分成兩個子區域?+(t)和??(t)，且滿足

其中，Γ(t)是兩種互不相溶介質之間的相界面.

利用有限體積方法求解每種介質的控制方程組，采用Level Set 方法捕捉不同介質之間的相界面，并通過在界面法向上求解局部一維多介質Riemann 問題的精確解來計算相界面上的數值通量.整個計算流程如下.

2.1 相界面追蹤

Level Set 方法把隨時間運動的相界面Γ(t)定義為距離函數的零等值面[30].本文利用特征線方法求解Level Set 函數的演化方程

根據特征線方法的思想，只要知道tn+1時刻網格頂點Xi上的流體在tn時刻的位置x(tn+1)，就可以得到新的Level Set 函數在Xi上的值.對每個頂點Xi，利用特征線方法可得

再利用代數平均，可得新的Level Set 函數

利用顯式正系數格式[31]求解如下所示的重新初始化方程

2.2 相界面幾何重構

當得到tn時刻的Level Set 函數后，根據其分片線性的特征重構出界面單元內分片線性的相界面，包括：二維兩相、二維三相、三維兩相和三維三相情形，具體如下.

2.2.1 二維兩相情形

圖1 給出了二維情形下單元內包含兩種介質時的相界面幾何結構.假設在三角形?ABC中，頂點A;B;C處的Level Set 函數值分別為.如果某個頂點的函數值與其他兩個頂點具有相反的符號，該單元一定為相界面單元，且該單元中存在兩種介質.

圖1 二維相界面拓撲結構(兩相情形)Fig.1 Two-dimensional topology of phase interface(two-phase case)

不失一般性，假設在三角形?ABC中，C點的相編號為I，A,B兩點的相編號為II，且存在相界面EF.根據水平集函數分片線性的特征，可計算E;F的坐標

且第I 相、第II 相所占的面積分別為

2.2.2 二維三相情形

采用向量型水平集函數來描述超過兩相情形時的相界面.對任一點x，向量型水平集函數在該點的取值為，如果其最大分量為，則該點的相編號為i.其中，n表示相數.

如果三角形ABC中3 個頂點的相編號各不相同，該單元中存在3 種介質，如圖2 所示.不失一般性，假設A;B;C三點的相編號分別為I,II,III.根據水平集函數的分片線性特征，可得單元3 條邊上的相界面分界點D;E;F的位置

圖2 二維相界面拓撲結構(兩相情形)Fig.2 Two-dimensional topology of phase interface(three-phase case)

下面確定三相分界點G，記

采用如下近似

來確定三相分界點的位置，然后計算各相所占的面積和兩兩之間的相界面.其中G的位置為

此時，第I 相所占的面積為A?ADE+A?DEG，第II 相所占的面積為A?BDF+A?DFG，第III 相所占的面積為A?CEF+A?EFG.I,II 兩相的相界面為線段DG，I,III 兩相的相界面為線段EG，II,III 兩相的相界面為線段FG.

2.2.3 三維兩相情形

根據水平集函數與三維四面體網格的相交情況，可以將三維兩相的相界面幾何結構分為兩種情形，分別如圖3(a)和圖3(b)所示.

圖3 三維相界面拓撲結構(兩相情形)Fig.3 Three-dimensional topology of phase interface(two-phase case)

(1)情形1

四面體ABCD中有一個頂點與其他3 個頂點的相編號不一樣，不妨設A點的相編號為I，B,C,D三點的相編號為II(如圖3(a)所示)，此時

根據水平集函數的分片線性特征，使用和二維兩相情形類似的算法可得四面體棱上的相分界點E,F,G，此時相界面為?EFG，且第I 相所占的體積為

第II 相所占的體積為

(2)情形2

四面體中4 個頂點按相編號平均分為兩組.不妨設A,B的相編號為I，C,D的相編號為II，如圖3(b)所示.使用和情形1 類似的算法，可以計算出四面體棱上的相分界點E,F,G,H，此時相界面為四邊形EFGH.記

則II 相所占的體積為

I 相所占的體積為

2.2.4 三維三相情形

圖4 給出了三維情形下單元內包含3 種介質時的相界面拓撲結構.根據水平集函數的分片線性性質，四面體網格的四個頂點中必然有兩個為一相，另外兩個頂點各占據一相.不失一般性，不妨設A的相編號為I，B的相編號為II，C,D的相編號為III.

圖4 三維相界面拓撲結構(三相情形)Fig.4 Three-dimensional topology of phase interface(three-phase case)

使用和前面類似的算法，可以計算出四面體棱上的相分界點，包括：相I 與相III 的分界點E,F，相II 與相III 的分界點G,H，以及相I 與相II 的分界點K，此時四面體單元ABCD內包含兩個相界面：三角形EFK和四邊形FGHK.整個四面體單元的體積為

I 相所占的體積為

II 相所占的體積為

III 相所占的體積為

2.3 數值通量

在Euler 坐標系下的多介質數值計算中，每個單元內守恒量的變化來自兩個部分：(1)單元內介質與單元外同相的介質在單元邊界上通過流入流出或相互作用力而帶來的守恒量變化，稱為單元邊界通量；(2)不同相的介質在相界面處由于相互作用力帶來的守恒量變化，稱為相界面通量.

以二維兩相情形為例(如圖5 所示)，假設單元Ki;n與單元Kj;n具有共邊Sij;n，且包含物質界 面將Ki;n和Sij;n依次分為2 個部分，此時單元包含2 種介質，每種介質的守恒量分別定義為

圖5 多介質流動計算模型示意圖Fig.5 Schematic diagram for the multi-media flow calculation model

(1)單元邊界數值通量

單元邊界的數值通量是指單元邊界上同種介質之間的數值通量，且滿足

(2)相界面數值通量

相界面數值通量是相界面兩側不同介質之間的數值通量，且滿足

2.4 守恒量更新

當計算得到Ki;n的單元邊界數值通量和相界面數值通量后，可對該單元中每種介質的守恒量進行如下更新

2.5 單元碎片聚合

利用式(7)來更新守恒量，等價于在原始背景網格Λ 和相界面耦合的網格上求解控制方程組.由于有些網格單元被相界面分割成一些碎片，這些碎片單元的尺寸比原始單元的尺寸小，導致CFL 條件被破壞，可能造成數值不穩定.為了消除不穩定性，需要對進行特殊處理.首先建立單元片，使得每個碎片單元都處于同介質的單元片中，并且該單元片中至少包含一個完整單元；然后修正守恒量，對單元片中每個單元上的介質守恒量進行單元體積的加權平均，以保證整體物理量的守恒.圖6 為網格聚合算法及碎片處理的示意圖.

圖6 網格聚合算法及碎片處理Fig.6 Cell aggregation and fragment process

2.5.1 建立單元片

提出如下的聚合算法來建立單元片.首先把包含碎片的單元(即界面單元)收集起來

再把與界面單元相鄰并且只含有當前所考慮介質的單元收集起來，作為種子單元

具體算法流程如下：

2.5.2 修正守恒量

可以證明，后處理方法(8)滿足

后處理過程(8)實際上消去了同一單元片內部單元間的數值通量，將過程(7)的更新改為只考慮單元片邊界上的通量.由于每個單元片中至少有一個完整單元，因而滿足CFL 條件，整個算法的數值穩定性得到了保證.

如果一個小的單元碎片孤立存在于其他介質單元中，有可能整個單元片中沒有一個完整的單元與它含有相同的介質.因此，有些在算法完成之后只有一個單元，將這些單元合起來作為一個單元片，稱之為“孤島”.如果孤島的體積小于預設的忍量(原始單元體積的10?8)，直接將孤島上該介質的守恒量設為缺省的極低壓狀態，計算過程中孤島會逐漸縮小直至消失，這種做法可以有效地消除過于復雜的相界面，因此總的計算量得到控制.

3 數值算例

本節給出一些數值算例來對數值方法進行測試，包括多介質Riemann 問題、激波與氣泡相互作用、觸地淺埋爆炸以及典型坑道內的沖擊波傳播等問題，驗證數值方法的可靠性和計算效率.

3.1 多介質Riemann 問題

首先計算一個JWL 狀態方程的單介質Riemann問題，來自文獻[34].計算區域為[0,1]×[0,5.0×10?3],網格尺寸為2.5×10?3.其中，JWL 狀態方程參數取A1=854:5 GPa，A2=20:5 GPa，=0:25，R1=4:6,R2=1:35，ρ0=1840 kg/m3.初值條件如下

計算時間為t=1:2×10?5.計算結果與文獻[34]參考解的對比如圖7 所示，對比表明本文結果與文獻中的高精度參考解符合較好，且在相界面處密度具有高分辨率.

再計算一個JWL--彈性固體的Riemann 問題，其中氣體相用JWL狀態方程來描述，參數為A1=854:5 GPa，A2=20:50 GPa，=0:25，R1=4:6,R2=1:35，ρ0=1840 kg/m3.固體相的靜水壓力用剛性氣體狀態方程描述，且γ=4:4，=600 MPa，固體的偏應力遵守Hooke 定律，且βe=1014Pa·kg/m3[35].初值條件如下

圖7 Shyue 激波管問題Fig.7 Shyue shock tube problem

圖7 Shyue 激波管問題(續)Fig.7 Shyue shock tube problem(continued)

計算時間為t=10?4.圖8 給出了計算結果與參考解[35]的對比，從對比分析可知數值結果與參考解符合較好，且在接觸間斷和激波附近具有較高的分辨率.

圖8 JWL-彈性固體Riemann 問題Fig.8 JWL-elastic solid Riemann problem

3.2 激波與氣泡相互作用

激波與氣泡作用問題是可壓縮多介質流體問題的一個經典算例[36-38].初值條件如下：一個直徑為50 mm 的圓形氣泡放置在長311.5 mm、寬89 mm 的激波管的正中間，其余地方充滿空氣(如圖9 所示).初始時刻，氣泡右側5 mm 處有一個馬赫數為Ma=1:22、向左運動的平面激波.初始物理參數如表2 所示.計算模型的上、下邊界設為反射邊界條件，左、右邊界設為出流邊界條件，網格尺寸為1 mm.

圖9 激波與氣泡相互作用問題示意圖Fig.9 Computational model of gas-bubble interaction problem

表2 激波與氣泡問題的初始物理參數Table 2 Initial physical parameters of gas-bubble problem

圖10 給出了不同時刻的數值結果，并與Haas 和Sturtevant 的實驗結果[38]進行了比較.圖中第一列為實驗紋影圖，第二列為數值計算得到的密度等值線圖.

從圖10 可以看出，數值計算結果與實驗結果定性吻合.當激波從密度大的空氣進入密度小的氦氣時，氣泡的上游邊緣呈現出Richtmyer-Meshkov 不穩定性，其蘑菇形狀在445μs 時刻發生.從427μs 開始氣泡發生嚴重變形，到674μs 時空氣的上尖峰和下尖峰形成.隨后氣泡變得越來越薄，直至發生破裂.算例中的氣泡大變形問題對數值計算提出了很大的挑戰，但計算結果表明，本文的數值方法可以很好的追蹤相界面，同時能夠保持y方向上解的對稱性，證明了數值方法的正確性.

圖10 激波與氣泡相互作用問題.第一列從上到下依次為32μs,53μs,62μs,72μs,82μs,102μs,245μs,427μs,674μs和983μs;第二列從上到下依次為23μs,43μs,53μs,66μs,75μs,102μs,260μs,445μs,674μs和983μsFig.10 Gas-bubble interaction problem.First column:32μs,53μs,62μs,72μs,82μs,102μs,245μs,427μs,674μs and 983μs;Second column:23μs,43μs,53μs,66μs,75μs,102μs,260μs,445μs,674μs and 983μs

圖10 激波與氣泡相互作用問題.第一列從上到下依次為32μs,53μs,62μs,72μs,82μs,102μs,245μs,427μs,674μs和983μs;第二列從上到下依次為23μs,43μs,53μs,66μs,75μs,102μs,260μs,445μs,674μs和983μs(續)Fig.10 Gas-bubble interaction problem.First column:32μs,53μs,62μs,72μs,82μs,102μs,245μs,427μs,674μs and 983μs;Second column:23μs,43μs,53μs,66μs,75μs,102μs,260μs,445μs,674μs and 983μs(continued)

3.3 淺埋爆炸問題

計算一個二維三相問題--淺埋爆炸，用于測試本文方法在處理三相問題時的能力.計算區域內包含3 種介質：爆炸產物、空氣和地介質.爆炸產物初始半徑為0.1 m，初始時刻位于地下0.1 m 處.爆炸產物和空氣采用理想氣體狀態方程，絕熱指數均為1.4.地介質的靜水壓力部分采用剛性氣體狀態方程，偏應力部分采用理想彈塑性材料模型，且γ=4:4，=6 kPa，μe=853 kPa，μp=0，Ye=6:5 kPa.初值條件如下

圖11 和圖12 分別給出了淺埋爆炸條件下，典型時刻整個計算區域內的密度和壓力分布情況.在淺埋爆炸問題中，炸藥爆炸產生的沖擊波會劇烈壓縮地介質，并在地介質中形成應力波.當沖擊波突破地表面后，形成彈坑，并噴射到空氣中形成強沖擊波.從計算結果來看，本文的數值方法能夠正確反映淺埋爆炸時的實際物理過程，表明本文的三相數值方法是正確的.

圖11 淺埋爆炸典型時刻的密度云圖Fig.11 Density contours of sub-surface blast problem

3.4 空中強爆炸沖擊波傳播

計算一個當量為1 kt TNT、爆高為200 m 的空中強爆炸沖擊波傳播問題.強爆炸產物采用等溫等壓球模型，取爆炸總當量的85%作為強爆炸的力學初始能量[39].空氣采用理想氣體狀態方程，初始壓力為101.3 kPa，初始密度為1.29 kg/m3.計算區域在x=0 和y=0 設為固壁反射邊界條件，其余邊界設為無反射邊界條件.

圖12 淺埋爆炸典型時刻的壓力云圖Fig.12 Pressure contours of sub-surface blast problem

計算得到典型時刻的沖擊波壓力等值線如圖13 所示，當空中強爆炸產生的沖擊波傳播到地面時，最早在爆心投影點發生正反射，然后以逐漸增大的入射角發生斜反射，產生雙激波結構的規則反射(圖13(a)).隨著沖擊波向外繼續傳播，入射角不斷增大.當入射角增大到臨界角附近時將發生馬赫反射，此時反射波陣面和入射波陣面的交點離開地面，形成垂直于地面的馬赫桿(圖13(b)).隨著沖擊波的進一步傳播，馬赫桿不斷增長.

圖14 給出了數值計算得到的地面沖擊波載荷分布(峰值超壓和沖量)與實測結果擬合的經驗公式[40]的對比情況.由對比可知，計算得到的地面沖擊波參數與實測結果符合得較好，最大誤差不超過40%.

圖13 空中強爆炸典型時刻的壓力云圖Fig.13 Pressure contours of air blast at typical time

圖13 空中強爆炸典型時刻的壓力云圖(續)Fig.13 Pressure contours of air blast at typical time(continued)

圖14 地面不同距離處的沖擊波參數Fig.14 Blast wave parameters on the ground at different radii

3.5 三維坑道中沖擊波傳播

計算了1 kg TNT 爆炸產生的沖擊波在一個“十”字形三維密閉空間中的傳播過程.考慮空間結構的對稱性，將其簡化為一個“L”型的結構.其中，E1=[0;2]×[0;2]×[0;6]m,E2=[0;6]×[0;2]×[0;2]m.TNT 和空氣分別采用JWL狀態方程和理想氣體狀態方程來描述.由于三維的計算量很大，本文采用了網格自適應技術并結合基于區域分解的并行計算策略，有效提高了計算效率.圖15 給出了典型時刻的沖擊波流場壓力云圖，圖16 給出了t=9:2×10?4s時刻典型位置處的沖擊波壓力切面圖.計算結果表明，TNT 爆炸后將在空氣中產生強沖擊波，起初波陣面以球面波的形式向外傳播.當沖擊波遇到固壁時，將發生反射并沿壁面繼續向前傳播，最終在整個密閉空間中形成非常復雜的載荷分布.

圖15 典型時刻的三維沖擊波壓力云圖和網格自適應圖Fig.15 Pressure contours and adaptive meshes of three-dimensional blast wave propagation

圖15 典型時刻的三維沖擊波壓力云圖和網格自適應圖(續)Fig.15 Pressure contours and adaptive meshes of three-dimensional blast wave propagation(continued)

圖16 t=9:2×10?4 s 時的沖擊波壓力切面圖Fig.16 Typical pressure sections of blast wave at t=9:2×10?4 s

4 結論

本文提出了一種基于Euler 坐標系的非結構網格上互不相溶的、具有銳利相界面的守恒型多介質流動數值方法，建立了二維和三維單純形網格上的多介質流動數值模擬程序，可用于模擬可壓縮流體和彈塑性固體在爆炸與沖擊等極端物理條件下的大變形動力學行為.本文提出的多介質流動數值方法采用了向量型的水平集函數，理論上可以處理全局任意種介質、局部三種介質的多介質流動問題，且能夠保持銳利的相界面.此外，提出的基于非結構網格上的聚合算法能夠從本質上消除由于網格被相界面分割成較小碎片而帶來的數值不穩定性，有效提高了數值計算的效率和健壯性.最后結合具有代表性的數值算例和實際應用問題，驗證了數值方法在處理多介質大變形和相界面大變形問題中的能力.