球坐標(biāo)系下多介質(zhì)混合物模型的數(shù)值模擬*

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

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

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

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

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

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

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

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

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

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

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

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

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

由于Γ、pref、eref都只和密度ρ 有關(guān)的函數(shù)，且式(5)對于任意的表達(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)不再成立。此時，對(6)作如下考慮：

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

1.3 輸送方程

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

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

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

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

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

第二，球坐標(biāo)系下的內(nèi)外介質(zhì)產(chǎn)生相互作用時，內(nèi)部的稀疏波(或激波)匯聚到一起后又反向朝著外部擴(kuò)散，原點處便形成了一個奇點。奇點處的稀疏波形成了獨特的反射機(jī)制。因此奇點的物理參數(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)可以改寫成：

其中：

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

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

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

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

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

2.3 輸送方程變換

由于每隔一個時間步，都會需要利用ρ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 代表流場中介質(zhì)的序號。

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

2.4 奇點處理

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

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

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

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

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

2.5 差分格式

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

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

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

2.6 時間空間步長

時間步長Δt，滿足收斂的條件為：

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

式中：H=E+p/ρ。

3 計算結(jié)果

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

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

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

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

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

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

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

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

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

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

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

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

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

在實際工程中，最受重視的是第一次的峰值。圖5 為沖擊波峰值壓力隨運動距離變化的曲線。作為對比，將給初始的體積分?jǐn)?shù)模型[2]和改進(jìn)后的質(zhì)量分?jǐn)?shù)模型同Zamyshlyayev 的經(jīng)驗公式[20]曲線作了對比。對比發(fā)現(xiàn)，改進(jìn)后的質(zhì)量分?jǐn)?shù)模型可以得到更為準(zhǔn)確的計算結(jié)果。在壓力隨運動距離發(fā)生急劇下降的過程中，利用體積分?jǐn)?shù)得到的計算結(jié)果，得到的結(jié)果比經(jīng)驗公式的結(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ù)計算模型分別得到的壓力時程曲線。從不同的計算結(jié)果中可以看出，利用改進(jìn)后的質(zhì)量分?jǐn)?shù)模型，得到的計算結(jié)果好于體積分?jǐn)?shù)模型。各個半徑處的壓力時程，質(zhì)量分?jǐn)?shù)模型均吻合得比較好。由于Zamyshlyayev 的經(jīng)驗公式僅考慮沖擊波壓力最大峰值未考慮壓力的二次峰值，本文也僅模擬了沖擊波壓力的最大峰值。

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

4 結(jié) 語

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