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

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

引 言

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

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

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

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

1 物理模型和數(shù)值方法

1.1 控制方程和狀態(tài)方程

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

(1)

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

TNT爆轟產(chǎn)物采用JWL狀態(tài)方程，如式(2)所示：

(2)

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

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

p=(γ-1)ρe

(3)

式中：γ=1.4為空氣的絕熱指數(shù)。

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

(4)

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

為后續(xù)進(jìn)行統(tǒng)一分析，將JWL狀態(tài)方程、理想氣體狀態(tài)方程和多項式狀態(tài)方程寫成統(tǒng)一形式：

(5)

1.2 多介質(zhì)流體數(shù)值方法

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

1.2.1 單元邊界數(shù)值通量

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

(6)

1.2.2 物質(zhì)界面數(shù)值通量

物質(zhì)界面數(shù)值通量是物質(zhì)界面兩側(cè)不同流體之間的數(shù)值通量，且滿足：

(7)

1.2.3 守恒量更新

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

(8)

2 多介質(zhì)Riemann問題數(shù)值求解

2.1 多介質(zhì)Riemann問題

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

Uτ+F(U)ξ=0

(9)

式中：τ和ξ分別表示時間和空間坐標(biāo)；

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

且滿足初值條件：

圖2 一維多介質(zhì)Riemann問題的波系結(jié)構(gòu)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)進(jìn)行迭代求解，可得到界面上的最終壓力p，此時界面上的速度u可表示為：

2.2 多介質(zhì)Riemann問題的數(shù)值求解

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

2.2.1 激波

當(dāng)p>pk時，該非線性波的類型為激波。根據(jù)內(nèi)能e的Rankine-Hugoniot關(guān)系：

(12)

(13)

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

(14)

2.2.2 稀疏波

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

(15)

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

2.2.3 Riemann問題的求解步驟

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

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

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

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

②如果min{pl,pr}

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

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

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

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

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

3 健壯性與精度測試

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

3.1 Riemann問題的健壯性測試

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

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

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

3.2 氣-水多介質(zhì)Riemann問題

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

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

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

4 在爆炸沖擊波計算中的應(yīng)用

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

4.1 水下爆炸沖擊波數(shù)值模擬

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

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

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

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

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

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

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

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

圖6 TNT空中爆炸的沖擊波壓力云圖和自適應(yīng)網(wǎng)格圖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 柱形容器內(nèi)沖擊波反射數(shù)值模擬

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

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

圖8 柱形容器內(nèi)反射沖擊波壓力云圖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 結(jié) 論

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

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

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

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