吳宗鐸，趙 勇，嚴(yán) 謹(jǐn)，宗 智，高 云

（1.廣東海洋大學(xué)海洋工程學(xué)院，廣東 湛江 524088；2.大連海事大學(xué)交通運(yùn)輸工程博士后流動(dòng)站，遼寧 大連 116026；3.大連理工大學(xué)船舶工程學(xué)院，遼寧 大連 116024；4.西南石油大學(xué)油氣藏地質(zhì)及開發(fā)工程國(guó)際重點(diǎn)實(shí)驗(yàn)室，四川 成都 610500）

多介質(zhì)的激波間斷問題一直以來都是計(jì)算流體力學(xué)方面的一個(gè)備受關(guān)注的問題，目前廣泛應(yīng)用的一個(gè)處理方法是利用Level Set 函數(shù)追蹤氣水界面，同時(shí)結(jié)合近似Riemann 解下的虛擬流體質(zhì)點(diǎn)完成對(duì)沖擊波參數(shù)的計(jì)算[1]。雖然得到的結(jié)果較理想，但是計(jì)算過程復(fù)雜，且計(jì)算依賴近似Riemann 解，無法在JWL (Jones-Wilkins-Lee)及Mie-Grüneisen 等復(fù)雜形式狀態(tài)方程下展開更深入的研究[2-3]。

相比之下，利用多介質(zhì)混合模型，可以將流場(chǎng)視為一個(gè)整體并用體積分?jǐn)?shù)或質(zhì)量分?jǐn)?shù)來區(qū)分流體，從而在不直接求近似Riemann 解的前提下完成計(jì)算。多介質(zhì)混合物模型早期由Abgrall 提出[4-5]，通過非守恒方程來計(jì)算由氣體或低壓縮度的固液體組成的流場(chǎng)[5-6]。隨后，Abgrall 將該混合物模型擴(kuò)展成一個(gè)七方程模型，能適應(yīng)一般形式狀態(tài)方程的流場(chǎng)[7]。柏勁松等[8]也利用此類多介質(zhì)混合物方程來模擬兩種以上介質(zhì)的沖擊過程，且配合等效方程可以適應(yīng)一些復(fù)雜形式的狀態(tài)方程[9]。梁珊等[10]和Liang 等[11]將Abgrall 的七方程模型與多GPU 的并行計(jì)算相結(jié)合來提高計(jì)算精確度。為了讓算法更好的適應(yīng)并行計(jì)算，劉娜等[12]將譜體積方法運(yùn)用到多介質(zhì)混合物模型中，提升了計(jì)算效率。

雖然該多介質(zhì)混合物模型在采用直角坐標(biāo)系處理平面波的問題時(shí)，取得不錯(cuò)的進(jìn)展，但是對(duì)于常用于處理球面波的球坐標(biāo)系下的多介質(zhì)問題，仍然具備一定的局限性。第一，多介質(zhì)混合物模型建立在速度和壓力在界面處維持平衡的基礎(chǔ)，僅對(duì)平面波適用，但是球面波的壓力和速度都會(huì)隨著傳播半徑的增加而衰減；第二，球坐標(biāo)系下，原點(diǎn)位置處為沖擊波運(yùn)動(dòng)的奇點(diǎn)位置，必須作出合適的處理；第三，為了增加工程實(shí)用性，計(jì)算模型需要能適應(yīng)復(fù)雜形式的狀態(tài)方程，而不僅僅是理想氣體方程等形式簡(jiǎn)單的狀態(tài)方程。

考慮到以上因素，本文將直角坐標(biāo)系下的基于Mie-Grüneisen 狀態(tài)方程下的多介質(zhì)混合模型延伸到球坐標(biāo)系下。為了壓力和速度不平衡問題以及坐標(biāo)原點(diǎn)的奇點(diǎn)問題，在原計(jì)算模型下作了參數(shù)修正并用質(zhì)量分?jǐn)?shù)替換體積分?jǐn)?shù)，并對(duì)奇點(diǎn)處數(shù)值差分作了特殊處理—利用相鄰網(wǎng)格點(diǎn)參數(shù)賦值給奇點(diǎn)參數(shù)但奇點(diǎn)處速度設(shè)定為0。最后，將通過數(shù)值算例驗(yàn)證其數(shù)值穩(wěn)定性和計(jì)算準(zhǔn)確性。

1 直角坐標(biāo)系下的計(jì)算模型

考慮多介質(zhì)的可壓縮流體組成的混合流場(chǎng)，其運(yùn)動(dòng)形式可以用歐拉方程來描述：

式中：ρ、p、u、E 分別為密度、壓力、速度和單位質(zhì)量能量[13]。單位質(zhì)量的總能量為內(nèi)能和動(dòng)能的總和，且E=e+u2/2。

流場(chǎng)中介質(zhì)的物理特性，可以用Mie-Grüneisen 狀態(tài)方程來表述：式中：Γ 為Mie-Grüneisen 系數(shù)，pref和eref為參考點(diǎn)的壓力和單位質(zhì)量?jī)?nèi)能，與密度ρ 有關(guān)。

1.1 界面處熱力學(xué)關(guān)系

界面處，流體的壓力p 和速度u 均保持平衡[4]，此時(shí)p 與u 在界面處不發(fā)生變化，而密度ρ 與內(nèi)能e 則出現(xiàn)數(shù)值間斷。因此?p/?x 和?u/?x 約為0，式(1)中的質(zhì)量守恒與能量守恒可以分別表示為：

將式(2)代入式(4)，有：

由于Γ、pref、eref都只和密度ρ 有關(guān)的函數(shù)，且式(5)對(duì)于任意的表達(dá)式pref(ρ)和eref(ρ)均成立。這里首先假設(shè)pref、eref都為0，那么有[2,9]

同理，可依次類推得到關(guān)于pref和eref的非守恒形式的方程：

1.2 波面處熱力學(xué)關(guān)系

由于式(6)～(8)的推導(dǎo)過程需要基于p 和u 在界面處連續(xù)這一前提條件，但是在左右波面處，p 和u 均出現(xiàn)間斷面，式(3)和(4)不再成立。此時(shí)，對(duì)(6)作如下考慮：

因此，在左右波面處雖然p 和u 均間斷，但由于在同類介質(zhì)中1/Γ 等參數(shù)的導(dǎo)數(shù)保持連續(xù)，因此，式(6)依然成立。同理，式(7)和(8)在左右波面處也成立。

1.3 輸送方程

對(duì)于多介質(zhì)流場(chǎng)，必須利用輸送方程來捕捉界面的運(yùn)動(dòng)。流場(chǎng)中，某種介質(zhì)的界面隨時(shí)間變化的運(yùn)動(dòng)方程，稱為輸送方程。一般利用體積分?jǐn)?shù)zi來建立介質(zhì)的輸送方程：

式中：i 為介質(zhì)的種類序號(hào)。非守恒變量1/Γ、pref/Γ、ρeref均為獨(dú)立變量，但它們對(duì)的ρ 偏導(dǎo)數(shù)可表達(dá)成關(guān)于體積分?jǐn)?shù)zi的加權(quán)疊加。這樣，綜合式(1)、(6)～(8)和(10)，可完成直角坐標(biāo)系下的多介質(zhì)流場(chǎng)的求解。

2 球坐標(biāo)系下的計(jì)算模型

球坐標(biāo)系下，則需要考慮更多的問題。

第一，在采用直角坐標(biāo)系下計(jì)算平面波問題時(shí)，p 與u 在界面處維持不變且偏導(dǎo)數(shù)?p/?x 與?u/?x 基本為0。該數(shù)值特性一直從界面處向左右延伸，直到波面處才形成間斷，如圖1 所示。但是采用球坐標(biāo)系處理球面波問題時(shí)，當(dāng)沖擊波沿著半徑r 向外擴(kuò)散時(shí)，傳播的形式則不再是平面沖擊波，而是球面沖擊波。此時(shí)，雖然p 與u 依然在界面處連續(xù)，但是球面波的擴(kuò)散過程中p 與u 會(huì)隨著半徑增加而發(fā)生衰減，并非保持不變。因此，球面沖擊波在界面處的偏導(dǎo)數(shù)?p/?r 與?u/?r 不再是0，而是一個(gè)穩(wěn)定的負(fù)數(shù)。這種情況下，式(3)～(4)非守恒性方程的建立的理論基礎(chǔ)并不十分準(zhǔn)確。

第二，球坐標(biāo)系下的內(nèi)外介質(zhì)產(chǎn)生相互作用時(shí)，內(nèi)部的稀疏波(或激波)匯聚到一起后又反向朝著外部擴(kuò)散，原點(diǎn)處便形成了一個(gè)奇點(diǎn)。奇點(diǎn)處的稀疏波形成了獨(dú)特的反射機(jī)制。因此奇點(diǎn)的物理參數(shù)變化需仔細(xì)考慮。

圖 1 平面沖擊波界面處p 的數(shù)值特征Fig.1 The numerical property of p at the interfacefor plane shock wave

2.1 球坐標(biāo)系下控制方程

這里將球坐標(biāo)系下的控制方程表示成守恒變量U 與通量F 的關(guān)系為[13-14]:

其中：

式中：系數(shù)α=2。式(11)可以改寫成：

其中：

此時(shí)，將非守恒形式的方程組添加到式(11)中，即可得到球坐標(biāo)系下的控制方程。

2.2 熱力學(xué)參數(shù)修正

需要注意的是，非守恒形式的方程式(6)～(8)的推導(dǎo)，是以界面處于平衡狀態(tài)為前提的。在球坐標(biāo)系下，由于沖擊波是沿著球面向四周擴(kuò)散。隨著擴(kuò)散范圍的增加，幾乎所有的沖擊波參數(shù)都會(huì)隨著半徑增加而降低。此時(shí)，p 與u 在界面處將呈現(xiàn)出持續(xù)的減小趨勢(shì)。這樣，式(6)～(8)無法直接使用到式(12)中。在直角坐標(biāo)系和球坐標(biāo)系下，界面間斷的關(guān)系分別為

考慮到三個(gè)非守恒變量1/Γ、pref/Γ、ρeref均為密度ρ 的函數(shù)，這里每隔一個(gè)時(shí)間步對(duì)它們做修正:

式中：ρi為介質(zhì)i 的密度。而的其他參數(shù)如ρ、p、u、E 的計(jì)算，仍然借助當(dāng)前時(shí)間步下的Γ、pref、eref。而式(6)～(8)則保留在系統(tǒng)方程中，并且每隔一個(gè)時(shí)間步，非守恒變量1/Γ、pref/Γ、ρeref都會(huì)得到修正。當(dāng)時(shí)間步足夠密時(shí)，非守恒形式在球坐標(biāo)系下帶來的誤差可以控制在一定的范圍內(nèi)。

2.3 輸送方程變換

由于每隔一個(gè)時(shí)間步，都會(huì)需要利用ρi進(jìn)行修正。這里用質(zhì)量分?jǐn)?shù)代替體積分?jǐn)?shù)來構(gòu)造輸送方程。利用式(12)中的質(zhì)量守恒關(guān)系可得：

根據(jù)式(14)得到的質(zhì)量分?jǐn)?shù)yi可以很容易得到各介質(zhì)的密度ρi(i=1, 2, ···)。而非守恒方程(6)～(8)中的偏導(dǎo)數(shù)，則轉(zhuǎn)換成與質(zhì)量分?jǐn)?shù)相關(guān)的加權(quán)疊加：

式中：φ、φ、ψ 為非守恒變量1/Γ、pref/Γ、ρeref偏導(dǎo)數(shù)，下角標(biāo)1、2 代表流場(chǎng)中介質(zhì)的序號(hào)。

相比式(10)給出的基于體積分?jǐn)?shù)的輸送方程，基于質(zhì)量分?jǐn)?shù)的輸送方程能考慮到每一種介質(zhì)的密度ρi，從而得到更為合理的偏導(dǎo)數(shù)φ、φ、ψ。

2.4 奇點(diǎn)處理

對(duì)于球坐標(biāo)系來說，奇點(diǎn)是一個(gè)無法避免的問題。參考文獻(xiàn)[15]中對(duì)奇點(diǎn)的處理方式，奇點(diǎn)處質(zhì)點(diǎn)的運(yùn)動(dòng)速度為0，即u=0，但是奇點(diǎn)處的質(zhì)量和能量守恒關(guān)系依然成立：

展開后，將u=0 代入，有：

由于在奇點(diǎn)處主要涉及沖擊波的匯聚和發(fā)散，不涉及界面的運(yùn)動(dòng)，因此，可不考慮輸送方程。參數(shù)Γ、pref、eref由密度ρ 直接計(jì)算得到。

文獻(xiàn)[15]中，對(duì)偏導(dǎo)數(shù)?u/?r 采用了如下的計(jì)算方式：

式中：下角標(biāo)1 和2 代表網(wǎng)格點(diǎn)序號(hào)。但考慮到本文的方程體系非守恒性較強(qiáng)，將u1和u2都替換為奇點(diǎn)處的速度0，這樣偏導(dǎo)數(shù)?u/?r=0。再根據(jù)式(17)計(jì)算得到奇點(diǎn)處的其他參數(shù)。

2.5 差分格式

在有限體積差分下，第j 個(gè)網(wǎng)格點(diǎn)的變量Uj，可按如下格式計(jì)算：

式中：Δt 和Δr 分別為時(shí)間和空間步長(zhǎng)，上標(biāo)(n)代表第n 步時(shí)間序號(hào)。對(duì)于物理量第j 和j+1 個(gè)網(wǎng)格點(diǎn)之間的通量和 源項(xiàng)分別為：

式中：角標(biāo)j+1/2 表示通量中的參數(shù)，角標(biāo)j 表示網(wǎng)格點(diǎn)的參數(shù)；在球面沖擊波中，通量表示沿半徑r 方向通過j 點(diǎn)的物理變量。

2.6 時(shí)間空間步長(zhǎng)

時(shí)間步長(zhǎng)Δt，滿足收斂的條件為：

式中：CCFL是為了保證迎風(fēng)型計(jì)算格式收斂而設(shè)置的一個(gè)0～1 之間的系數(shù)，滿足收斂的條件為CCFL＜1,并且當(dāng)CCFL取得越小時(shí)間步長(zhǎng)越小，計(jì)算越精細(xì), 這里取CCFL=0.4； uj和cj表示第j 個(gè)網(wǎng)格點(diǎn)的質(zhì)點(diǎn)速度和聲速。聲速cj的表達(dá)式為：

式中：H=E+p/ρ。

3 計(jì)算結(jié)果

3.1 無因次參數(shù)的氣水作用算例

這里考慮一個(gè)無因次的氣水作用問題[16]。初始時(shí)刻，中心氣團(tuán)的物理狀態(tài)為：p=8.3×103，ρ=1.27，u=0，外部水的物理狀態(tài)為：p=1.0，ρ=1.27，u=0。氣體和水的狀態(tài)方程分別為

式中：γ 和B 為描述水和氣體在低壓縮度時(shí)的熱力學(xué)參數(shù)。氣團(tuán)初始半徑為R0。由于中心氣團(tuán)的壓力和密度高于周邊水的。因此，高壓的氣團(tuán)會(huì)在水中形成沖擊波并推動(dòng)界面向外部擴(kuò)散。對(duì)沖擊波波面和介質(zhì)界面的運(yùn)動(dòng)位置的變化情況進(jìn)行記錄（波面為水中壓力的最大值位置，界面位置為從外向內(nèi)氣體的質(zhì)量分?jǐn)?shù)剛超過0.5 的位置），并繪制了沖擊波和界面隨時(shí)間變化的曲線，如圖2 所示。從圖2 可以看出，本文所得到的沖擊波和界面的運(yùn)動(dòng)位置與Liu 等[16]和Flores 等[17]的結(jié)果都比較接近。三個(gè)數(shù)據(jù)的計(jì)算結(jié)果，在氣水作用的初期差別不大。但是隨著沖擊波和界面向外擴(kuò)散，本文的計(jì)算結(jié)果與Flores 等[17]的計(jì)算結(jié)果相比，誤差略偏大，但仍然在一個(gè)合理范圍內(nèi)。

圖 2 沖擊波與界面位置隨時(shí)間的變化Fig.2 Time evolutions of position of shock wave and interface

圖 3 沖擊波到達(dá)3 倍R0 時(shí)，壓力與速度的分布Fig.3 The distributions of pressure and velocity when shock wave reaches 3R0

圖3 為沖擊波到達(dá)3 倍氣團(tuán)初始半徑時(shí)候的壓力和速度曲線。此時(shí)，向內(nèi)收縮的稀疏波已達(dá)到中心原點(diǎn)處。由于氣團(tuán)已擴(kuò)散開，中心處的壓力開始逐漸下降，但中心處的速度仍然為0。速度曲線中，由中心向外的第一個(gè)拐點(diǎn)為稀疏波的端點(diǎn)，第二個(gè)拐點(diǎn)則為界面所在位置。另外，從圖3 中可以看出，處于數(shù)值奇點(diǎn)的原點(diǎn)處，數(shù)值穩(wěn)定性較好。證明本文所用的奇點(diǎn)處理方式能效果良好。

3.2 圓形藥球的水下爆炸算例

假設(shè)一TNT 的實(shí)心藥球質(zhì)量為50 kg，將空間步長(zhǎng)取為炸藥球半徑的1/100。利用改進(jìn)后的多介質(zhì)混合模型計(jì)算球坐標(biāo)系下的沖擊波運(yùn)動(dòng)情況。

炸藥的狀態(tài)方程為JWL 狀態(tài)方程[3]：

式中：A、B、R1、R2、ω 均為常數(shù)，由圓筒試驗(yàn)標(biāo)定得到[3]；θ=ρ/ρ0，ρ0為物質(zhì)的初始密度。水的多項(xiàng)式狀態(tài)方程為[18]：

式中：μ=θ-1；a1～a3, b0～b3為展開后的多項(xiàng)式系數(shù)。狀態(tài)方程(23)～(24)的參數(shù)取值見表1[19]。

表 1 狀態(tài)方程參數(shù)Table 1 Coefficients of EOS

本算例中，將計(jì)算范圍擴(kuò)大至20 倍以上的藥球半徑。這樣，可以得到?jīng)_擊波運(yùn)動(dòng)到不同位置處的壓力，如圖4 所示。從圖4 可以看出，由于在偏導(dǎo)數(shù)φ、φ、ψ 的處理方式差異，兩種計(jì)算模型得到的壓力分布有差別。但都具備較好的數(shù)值穩(wěn)定性。當(dāng)沖擊波運(yùn)動(dòng)到較遠(yuǎn)位置時(shí)，會(huì)出現(xiàn)多個(gè)壓力峰值。

在實(shí)際工程中，最受重視的是第一次的峰值。圖5 為沖擊波峰值壓力隨運(yùn)動(dòng)距離變化的曲線。作為對(duì)比，將給初始的體積分?jǐn)?shù)模型[2]和改進(jìn)后的質(zhì)量分?jǐn)?shù)模型同Zamyshlyayev 的經(jīng)驗(yàn)公式[20]曲線作了對(duì)比。對(duì)比發(fā)現(xiàn)，改進(jìn)后的質(zhì)量分?jǐn)?shù)模型可以得到更為準(zhǔn)確的計(jì)算結(jié)果。在壓力隨運(yùn)動(dòng)距離發(fā)生急劇下降的過程中，利用體積分?jǐn)?shù)得到的計(jì)算結(jié)果，得到的結(jié)果比經(jīng)驗(yàn)公式的結(jié)果略微偏大。

圖 4 沖擊波到達(dá)不同位置處的壓力曲線Fig.4 Pressure curves at the time instants when the shock wave reaches different positions

圖 5 沖擊波在不同位置處的壓力峰值變化Fig.5 Peak value variation of pressure of the shock wave at different positions

圖6 為利用體積分?jǐn)?shù)與質(zhì)量分?jǐn)?shù)計(jì)算模型分別得到的壓力時(shí)程曲線。從不同的計(jì)算結(jié)果中可以看出，利用改進(jìn)后的質(zhì)量分?jǐn)?shù)模型，得到的計(jì)算結(jié)果好于體積分?jǐn)?shù)模型。各個(gè)半徑處的壓力時(shí)程，質(zhì)量分?jǐn)?shù)模型均吻合得比較好。由于Zamyshlyayev 的經(jīng)驗(yàn)公式僅考慮沖擊波壓力最大峰值未考慮壓力的二次峰值，本文也僅模擬了沖擊波壓力的最大峰值。

圖 6 不同計(jì)算模型下壓力隨時(shí)間變化曲線Fig.6 Time evolutions of pressure for different numerical models

4 結(jié) 語(yǔ)

針對(duì)球坐標(biāo)下多介質(zhì)界面捕捉困難的問題，本文對(duì)原始的Mie-Grüneisen 多介質(zhì)混合物模型作了熱力學(xué)參數(shù)修正、輸送方程變換、特殊奇點(diǎn)處理等多項(xiàng)數(shù)值改進(jìn)。改進(jìn)后的基于質(zhì)量分?jǐn)?shù)的模型，可以在常規(guī)狀態(tài)方程下準(zhǔn)確捕捉界面和沖擊波的運(yùn)動(dòng)，且可以在適應(yīng)JWL 等形式復(fù)雜的狀態(tài)方程。而適當(dāng)?shù)钠纥c(diǎn)處理也保證了數(shù)值計(jì)算在奇點(diǎn)處的穩(wěn)定性。另外，數(shù)值算例證明，相比原始的體積分?jǐn)?shù)模型，改進(jìn)后的質(zhì)量分?jǐn)?shù)模型可以得到較為準(zhǔn)確的計(jì)算結(jié)果。