基于Riemann問題的可壓縮多介質流體數值模擬及在二維爆炸沖擊波傳播問題中的應用

姚成寶，付梅艷,2，閆 凱,2，韓 峰

(1.西北核技術研究院，陜西 西安 710024; 2.北京大學 數學科學學院，北京 100871)

引 言

可壓縮多介質流動問題的數值模擬越來越受到人們的廣泛關注，諸如爆炸沖擊波、高速沖擊等，都屬于可壓縮多介質流動問題的范疇。在多介質流動問題的數值模擬中，計算區域內可能會出現多種介質相互作用。界面兩側的介質通常具有截然不同的狀態方程，且介質的初始狀態也會存在較大差異，給數值模擬特別是界面附近的數值計算帶來很大的困難。很多對單介質問題比較有效的數值方法在處理多介質問題時往往無法得到理想的結果，并且在界面處會出現非物理振蕩，降低計算精度，甚至使計算難以進行。如何準確描述界面的位置以及界面兩側介質之間的相互作用，是求解多介質流動問題的關鍵。

在Euler坐標系下，目前應用較多、用于描述和追蹤物質界面的方法包括兩類：彌散界面方法和清晰界面方法。其中，彌散界面方法[1-2]將包含界面在內的幾層網格上的物理量進行平滑磨光，形成具有一定寬度的平衡層。清晰界面方法將界面視為嚴格的接觸間斷，并清晰捕捉界面的位置，例如VOF方法[3]、Level Set方法[4]和Front Tracking方法[5]等。當確定物質界面的位置后，需進一步準確處理不同介質間的相互作用，常見的處理方法包括虛擬流體類方法[6-7]和切割單元法(Cut cell method)[8]。無論采用何種方法，如何確定界面上的狀態參數是一個關鍵環節。由于物質界面實際上是一個接觸間斷，求解界面上的多介質Riemann問題精確解[9]是求解該類問題的最直接和有效的方法，其能夠精確捕捉激波和接觸間斷的位置。

Riemann問題的求解和介質的狀態方程密切相關，不同狀態方程的Riemann問題求解也各不相同。對于形式比較簡單的狀態方程，Riemann問題可以表示為解析形式的代數方程并通過迭代得到任意精度的解，但工程應用中常見的狀態方程通常具有比較復雜的形式，尚不存在Riemann問題精確解的通用計算方法，大多采用近似解來進行計算。例如，Banks[10]對JWL 狀態方程的 Riemann 問題進行了分析，并利用二分法迭代求解接觸間斷上的壓力p滿足的代數方程。Kamm[11]在 Banks的基礎上，針對滿足凸性的一般狀態方程建立了類似的數值求解方法。Rallu[12]、FarHat[13]和Lee[14]分別針對JWL 狀態方程和Mie-Grüneisen形式的狀態方程，建立了一個2×2的代數方程組，并通過迭代求解該方程組，獲得了 Riemann 問題的近似解。

本研究在上述基礎上，提出了一種針對JWL、多項式等高度非線性狀態方程的多介質Riemann問題求解方法，能夠對激波、稀疏波作用到不同物質界面時的復雜波系結構進行高效、準確的求解。結合前期工作[15-16]，建立了一套能夠模擬極端條件下、具有復雜狀態方程的多種介質間相互作用的二維數值計算體系。最后，對Riemann問題的健壯性和精度進行了測試，并對水下、空氣自由場以及密閉空間下的TNT爆炸沖擊波傳播過程進行了數值模擬，驗證了數值方法的實際應用能力。

1 物理模型和數值方法

1.1 控制方程和狀態方程

考慮二維計算區域上的多介質流體問題，基于Euler坐標系描述流體的動力學行為，控制方程組如式(1)所示：

(1)

式中：ρ為密度；u、v為x、y方向的速度；E為單位體積總能量；p為壓力，下同。

TNT爆轟產物采用JWL狀態方程，如式(2)所示：

(2)

式中：ρ0為初始密度；e為單位質量比內能，且滿足e=E/ρ-(u2+v2)/2；A、B、R1、R2和ω為JWL狀態方程的參數，且A=3.712×1011Pa，B=3.23×109Pa，R1=4.15，R2=0.95，ω=0.3。

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

p=(γ-1)ρe

(3)

式中：γ=1.4為空氣的絕熱指數。

水采用多項式狀態方程來描述，如式(4)所示：

(4)

式中：μ=ρ/ρ0-1；A1、A2、A3、B0、B1、Γ1、Γ2為多項式狀態方程的系數，且A1=2.2×109Pa，A2=9.54×109Pa，A3=1.457×1010Pa，B0=0.28，B1=0.28，Γ1=2.2×109Pa，Γ2=0Pa。

為后續進行統一分析，將JWL狀態方程、理想氣體狀態方程和多項式狀態方程寫成統一形式：

(5)

1.2 多介質流體數值方法

圖1 多介質流體計算模型示意圖Fig. 1 Schematic diagram of multi-media fluid calculation model

1.2.1 單元邊界數值通量

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

(6)

1.2.2 物質界面數值通量

物質界面數值通量是物質界面兩側不同流體之間的數值通量，且滿足：

(7)

1.2.3 守恒量更新

當計算得到Ki,n的單元邊界數值通量和物質界面數值通量后，可對該單元中每種流體的守恒量進行式(8)所示的更新：

(8)

2 多介質Riemann問題數值求解

2.1 多介質Riemann問題

由于多維Riemann問題的復雜性，相關的理論研究十分困難，目前常見的做法是將多維Riemann問題簡化為界面法向上的局部一維Riemann問題[9]，利用界面兩側流體的壓力和法向速度的連續條件來求解一維Riemann問題，此時Riemann問題可表示為求解式(9)所示的一維多介質Euler方程組的初值問題：

Uτ+F(U)ξ=0

(9)

式中：τ和ξ分別表示時間和空間坐標；

U=[ρ,ρu,E]Τ,F=[ρu,ρu2+p,(E+p)u]Τ,

且滿足初值條件：

圖2 一維多介質Riemann問題的波系結構Fig.2 Typical wave structure of one-dimensional Riemann problem

f(p)=fl(p,Ul)+fr(p,Ur)+ur-ul=0

(10)

其中，

(11)

式中：c表示聲速，k=l,r。

利用牛頓迭代法對式(10)進行迭代求解，可得到界面上的最終壓力p，此時界面上的速度u可表示為：

2.2 多介質Riemann問題的數值求解

對于JWL、多項式等形式復雜的狀態方程，fl(p,Ul)和fr(p,Ur)具有高度非線性，解析求解式(10)非常困難。本研究采用數值方法來近似求解fl(p,Ul)和fr(p,Ur)，并通過在數值求解過程中進行誤差控制來保證近似解收斂到式(10)的精確解。

2.2.1 激波

當p>pk時，該非線性波的類型為激波。根據內能e的Rankine-Hugoniot關系：

(12)

(13)

利用牛頓迭代法對式(13)進行迭代求解：

(14)

2.2.2 稀疏波

當p≤pk時，該非線性波的類型為稀疏波，p和ρ滿足等熵關系dρ/dp=1/c2。結合式(11)的稀疏波關系，可得：

(15)

利用自適應步長的Runger-Kutta-Fehlberg方法[17]對式(15)進行聯立求解，可完成稀疏波分支fk(p,Uk)的求解。其中c表示聲速，且滿足：

2.2.3 Riemann問題的求解步驟

綜上所述，本研究提出的Riemann問題的整體求解步驟如下：

(1)提供物質界面壓力的初始估計p0和整體誤差閾值ε0，其中

(2)假設第n步迭代的值pn已知，確定左、右非線性波的類型：

①如果pn>max{pl,pr}，兩側非線性波均為激波；

②如果min{pl,pr}

③如果pn≤min{pl,pr}，兩側非線性波均為稀疏波。

(3)根據非線性波的類型和當前迭代步的局部求值誤差εn，嚴格控制激波分支的非線性代數方程(13)的迭代殘差和稀疏波分支的常微分方程組(15)的局部求值誤差，計算得到當前步的fl(pn,Ul)和fr(pn,Ur)；

(4)更新得到n+1迭代步時，物質界面上的壓力：

(5)當壓力的相對變化達到指定的整體誤差ε0時，迭代終止，將得到的收斂解pn近似作為物質界面的壓力p*，否則返回步驟(2)。

(6)計算物質界面上的法向速度：

3 健壯性與精度測試

利用本研究數值方法分別對強稀疏波Riemann問題和氣-水多介質Riemann問題進行數值模擬，驗證數值方法的精度和健壯性。另外，若無特殊說明，本研究所有算例均采用kg-m-s單位制，其中密度ρ、速度u和壓力p的單位分別為kg/m3、m/s和Pa。

3.1 Riemann問題的健壯性測試

利用一個強稀疏波問題[9]對本研究設計的Riemann問題求解方法的健壯性進行測試。計算區域為[0,1]m，初始界面位于x=0.5m，界面兩側的介質均為理想氣體，且初值條件如下：

將初值條件進行無量綱縮放，測試數值方法在不同的數量級下能否始終保持與解析解一致。其中，無粘、可壓縮流體的無量綱關系滿足p∝u3和p∝ρ3，即當u和ρ進行相同的量級變化后，p相應的變化3個量級。表1給出了在不同量級的初值條件下，數值模擬結果與理論結果的對比，可知本研究的 Riemann 問題求解器在不同量級的初值條件下均具有較好的一致性，驗證了方法的健壯性。

表1 Riemann求解器的健壯性測試

3.2 氣-水多介質Riemann問題

計算一個氣-水多介質Riemann問題，驗證本研究方法求解多介質Riemann問題的精度。其中，氣體采用JWL狀態方程來描述，水采用多項式狀態方程來描述。計算區域為[0,1]m×[0, 5×10-3]m，網格尺寸為2.5×10-3m。初始條件如下：

計算時間為t=8.0×10-5s。數值計算結果與參考解的對比如圖3所示，對比結果表明本研究的數值方法能夠準確地捕捉Riemann問題解的波系結構，且在激波和物質界面附近沒有產生非物理震蕩，說明本研究方法能夠有效處理JWL和多項式等具有高度非線性的狀態方程。

圖3 JWL-多項式Riemann問題Fig.3 JWL-Polynomial Riemann problem

4 在爆炸沖擊波計算中的應用

進一步利用數值方法分別對TNT水下爆炸、TNT空中爆炸自由場、柱形容器內沖擊波反射等問題進行數值模擬，并將計算結果與實測數據進行比對，驗證數值方法在爆炸沖擊波問題中的應用能力。

4.1 水下爆炸沖擊波數值模擬

首先計算100kg TNT在水下自由場中爆炸產生的沖擊波傳播過程，分別采用JWL和多項式狀態方程來描述TNT爆炸產物和水介質。整個計算區域為[0,15] ×[0,4]m。為提高計算效率，利用H型網格自適應技術來捕捉沖擊波的波陣面。其中，原始的背景網格尺寸為30cm，初始時刻經過網格自適應加密后，TNT爆炸產物和水界面附近的最小網格尺寸約為1mm。整個計算區域的網格總量控制在(6～10)×104范圍內，從而在捕捉激波峰值的同時，能夠有效減少計算量，提高計算效率。

計算得到水下爆炸不同距離處的沖擊波峰值超壓和沖量，并與實驗結果[18]進行了對比，如圖4所示，兩者符合較好。

圖4 不同距離處的沖擊波峰值超壓和沖量Fig.4 Peak overpressures and impulses at different radii

下面計算一個淺水爆炸問題。計算區域為[1, 1]×[1, 1]m，將初始半徑為0.0527m的球形TNT放置在水面下方0.2m處，在初始時刻引爆TNT炸藥，計算爆炸產生的沖擊波在水中和水-氣界面中的傳播過程。計算模型包含3種介質：TNT 爆炸產物、水和空氣，同時也包含兩個物質界面，即爆炸產物與空氣的界面以及水與空氣之間的界面。其中，TNT、水和空氣分別采用 JWL、多項式和理想氣體狀態方程來描述。初值條件如下：

圖5為計算得到的不同時刻水下爆炸沖擊波壓力云圖。TNT爆炸后在水中產生沖擊波，并同時在水中形成氣泡。由于炸藥距離水面較近，在t=20s左右時，沖擊波先到達水面，并與自由水面發生作用，同時產生向空氣中傳播的透射激波和向水中傳播的反射波。由于氣泡附近水壓很大，因此遇到反射波后會產生一個向上運動的壓縮波。隨著時間的演化，氣泡附近的水壓逐漸變得復雜，氣泡的形狀在膨脹過程中也逐漸變成橢球型。計算結果表明，本數值模擬方法能夠對淺水爆炸等涉及多種具有高度非線性狀態方程的介質間相互作用過程進行定性正確的數值模擬，為后續進一步的定量研究奠定了基礎。

圖5 淺水爆炸典型時刻的沖擊波壓力云圖Fig.5 Pressure contours of shallow water explosion at typical time

4.2 TNT空中爆炸自由場沖擊波數值模擬

計算1kg TNT在空氣自由場中爆炸產生的沖擊波傳播過程，分別采用JWL和理想氣體狀態方程來描述TNT爆炸產物和空氣。圖6(a)、(b) 為典型時刻沖擊波壓力云圖和相應的網格自適應圖，圖7為計算得到的距爆心不同距離處的沖擊波峰值超壓和沖量，并與相關的實測結果[19-20]進行對比，兩者基本符合，驗證了計算模型和計算方法的正確性。

圖6 TNT空中爆炸的沖擊波壓力云圖和自適應網格圖Fig.6 Pressure contours and adaptive meshes of TNT explosion in free air at typical time

圖7 不同距離處的沖擊波峰值超壓和沖量Fig.7 Peak overpressures and impulses at different radii

4.3 柱形容器內沖擊波反射數值模擬

在爆炸近區，爆炸產物和沖擊波尚未完全分離。由于爆轟產物的密度遠高于空氣，此時需要精確模擬爆炸產物與空氣之間的相互作用，才能準確計算該區域內的真實載荷狀態。利用本研究方法計算了1kg TNT爆炸產生的近區沖擊波在遇到剛性壁面時的反射過程，得到了不同時刻的沖擊波反射壓力云圖，如圖8所示。

由圖8可知，當爆炸產生的沖擊波傳播到剛性壁面時，首先在爆心截面發生正反射，然后沿壁面向兩側繼續傳播，并發生規則反射。當入射角增大到臨界角附近時，沖擊波進一步發生了馬赫反射，最終形成比較復雜的馬赫波。圖9給出了計算得到的爆心截面處的沖擊波正反射峰值超壓和沖量，與文獻[20-22]以及實測數據[23]的結果比較吻合。

圖8 柱形容器內反射沖擊波壓力云圖Fig.8 Reflective pressure contours of TNT explosion in cylindrical vessel

圖9 不同距離處的正反射沖擊波峰值超壓和沖量Fig.9 Peak overpressures and impulses of normal reflection at different radii

5 結 論

(1)提出了一種處理高度非線性狀態方程的多介質Riemann問題的通用求解方法，能夠準確、高效處理多種復雜介質間的相互作用。

(2)結合Euler坐標系下具有清晰界面的可壓縮多介質流體數值方法，建立了一套能夠模擬具有高密度比、高壓力比以及復雜狀態方程的二維多介質流體計算體系。

(3)利用強稀疏波Riemann問題和氣-水多介質Riemann問題，驗證了數值方法的健壯性和精度。

(4)結合并行計算和網格自適應技術，對TNT爆炸沖擊波在水下、空氣自由場和柱形容器內的傳播過程進行了數值模擬，得到了與實測數據一致的數值結果，驗證了數值方法在爆炸沖擊波問題中應用的可靠性，為爆炸沖擊波在復雜環境下的傳播規律研究奠定了基礎。