陳凌蕙，徐 偉

（南昌航空大學(xué) 數(shù)學(xué)與信息科學(xué)學(xué)院，南昌 330063）

引 言

在對(duì)多介質(zhì)流進(jìn)行數(shù)值模擬時(shí)，有時(shí)會(huì)遇到流體的密度、壓力或內(nèi)能非常小的情況，如爆炸、高速射流及慣性約束核聚變等問題。此時(shí)用傳統(tǒng)的高精度格式，如WENO，DG等進(jìn)行求解，常會(huì)由于數(shù)值誤差而得到負(fù)的密度或壓力值，而這顯然是不符合實(shí)際情況的，同時(shí)這也會(huì)使得算法不穩(wěn)定，從而導(dǎo)致計(jì)算失敗。針對(duì)這類問題，目前已有一些研究工作，發(fā)展出了一些具有正保護(hù)功能的數(shù)值算法，但大多是針對(duì)單介質(zhì)流體的，而對(duì)于多介質(zhì)流的相關(guān)工作則相對(duì)較少。這是因?yàn)?，多介質(zhì)流體在介質(zhì)界面兩邊的密度與壓力等狀態(tài)常常會(huì)相差巨大，其所使用的狀態(tài)方程也不盡相同，甚至不同介質(zhì)使用的控制方程有時(shí)也不一樣，這些都使得構(gòu)造求解多介質(zhì)流動(dòng)問題的正保護(hù)數(shù)值算法要困難得多。解決這類問題的一個(gè)簡單的方法是強(qiáng)行將負(fù)的密度或壓力值等賦值為零，但這會(huì)使得算法的精度和守恒性都無法得到保證。Einfeldt在早期給出了一類求解可壓縮Euler方程的，具有正保護(hù)功能的Godunov格式[1]，但是只有一階精度。Perthame則分別給出了具有正保護(hù)功能的一階Lax-Friedrichs格式[2]和二階Boltzmann格式[3]。然而要想將這些格式中給出的正保護(hù)方法在不破壞精度的條件下，推廣運(yùn)用于任意高精度格式是非常困難的。Zhang和Shu在上述Lax-Friedrichs格式的基礎(chǔ)上給出了一種更易于推廣到高精度的正保護(hù)格式構(gòu)造方法[4]，并給出了判斷格式是否具有正保護(hù)性質(zhì)的充分條件，但這種方法只是適用于理想氣體。Cheng和Shu則結(jié)合HLLC方法和Lagrange方法，給出了一種針對(duì)多介質(zhì)流動(dòng)問題，且具有正保護(hù)功能的數(shù)值算法[5]。上述方法進(jìn)行正保護(hù)的大致思路都是通過構(gòu)造正保護(hù)限制器，進(jìn)而利用限制器對(duì)計(jì)算出現(xiàn)的負(fù)密度負(fù)壓力等值進(jìn)行修正，同時(shí)保證不破壞格式原有的精度及守恒性。

另外，在求解多介質(zhì)流問題時(shí)，若直接使用傳統(tǒng)的高精度格式，由于介質(zhì)界面兩邊流體狀態(tài)方程的不同，數(shù)值結(jié)果在界面附近會(huì)出現(xiàn)非物理振蕩或計(jì)算失敗。針對(duì)該問題，F(xiàn)edkiw等給出了一種GFM方法[6]，這種方法通過在界面附近定義虛擬流體，將多介質(zhì)問題轉(zhuǎn)化為多個(gè)單介質(zhì)問題，然后再使用已有的高精度方法進(jìn)行計(jì)算。最初的GFM方法在計(jì)算時(shí)，將虛擬流體的壓力與速度值取為當(dāng)?shù)財(cái)?shù)值，再由熵修正得到密度。但是在計(jì)算水-氣多介質(zhì)流時(shí)，該算法的穩(wěn)定性較差，從而產(chǎn)生了一些變形，如nGFM法[7]、mGFM法[8]以及rGFM法[9]等。其中mGFM法與rGFM法都是通過在介質(zhì)界面附近構(gòu)造黎曼問題，并利用黎曼問題的解來對(duì)界面附近的單元進(jìn)行特殊處理。不同的是，mGFM法只對(duì)虛擬流體中的單元進(jìn)行定義，而在真實(shí)流體中只對(duì)密度進(jìn)行等熵修正。這種做法實(shí)際上使得真實(shí)流體與虛擬流體在界面附近的單元沒有得到同步更新，從而在計(jì)算某些問題時(shí)還是會(huì)出現(xiàn)虛假振蕩。而rGFM法則利用黎曼問題的解同時(shí)對(duì)真實(shí)流體與虛擬流體進(jìn)行更新，能更真實(shí)的反映流體的特性。

由上述可知，在計(jì)算多介質(zhì)流問題時(shí)，有兩處需要進(jìn)行正保護(hù)，一是在界面處計(jì)算黎曼問題的時(shí)候，一是對(duì)各單介質(zhì)進(jìn)行計(jì)算的時(shí)候。本文基于正保護(hù)限制器[4]與正保護(hù)HLLC方法[5]，給出了一種適用于氣-氣及水-氣多介質(zhì)問題的正保護(hù)數(shù)值算法，并通過多個(gè)算例驗(yàn)證了算法的有效性。

1 方 程

1.1 控制方程

一維無粘可壓縮Euler方程：

其中U=(ρ,ρu,E)T，F(xiàn)(U)=(ρu,ρu2+p,(E+p)u)T，E=ρe+ρu2/2，式中的ρ是密度，u是速度，p是壓強(qiáng)，e是單位質(zhì)量的內(nèi)能，E是單位體積的總能。

二維無粘可壓縮Euler方程：

其中

式中的ρ，e與E的含義與一維情況一樣，而u是v則分別 表示x和y方向上的速度。

1.2 狀態(tài)方程

為了求解Euler方程，需要給出流體的狀態(tài)方程。當(dāng)流體為理想氣體時(shí)，其方程為：

而水的狀態(tài)方程有多種，本文取

其中γ與B為熱力學(xué)常數(shù)，將在具體算例中給出。不難看出氣體與水的狀態(tài)方程可以表示成統(tǒng)一的形式。此外，式（4）也可寫為p+B=(γ?1)(ρe?B)。

2 正保護(hù)數(shù)值方法

2.1 正保護(hù)有限體積法

對(duì)于一維情況，首先取點(diǎn)集G：

利用Jensen不等式易證G為凸集。而所謂正保護(hù)格式，即要使得當(dāng)Un∈G時(shí)，能夠有Un+1∈G。

考慮有限體積法的一般形式：

其中，m=ρu。接著取單元Ij上的Gauss-Lobatto點(diǎn)集

且有2N?3≥k。

由文獻(xiàn)[4]可知，對(duì)于有限體積格式（5），若在tn時(shí)刻，?α,j，有以及。則在CFL條件：

可得tn+1時(shí)刻的∈G。

根據(jù)上述結(jié)論，下面將給出正保護(hù)限制器對(duì)單元Ij上的重構(gòu)多項(xiàng)式qj(x)進(jìn)行修正，使得修正后的滿足：

修正將分兩步進(jìn)行，具體如下：

首先，密度修正：在各單元Ij上，將密度函數(shù)ρj(x)修正為，即取

其中

于是

本文采用WENO方法進(jìn)行空間離散，但由于WENO格式與DG格式不同，它不給出重構(gòu)多項(xiàng)式的具體形式，而是直接得到點(diǎn)值與。因此為了使用上述限制器對(duì)重構(gòu)多項(xiàng)式進(jìn)行修正，可以按文獻(xiàn)[10]中的方式得到重構(gòu)多項(xiàng)式，再對(duì)其進(jìn)行修正，繼而利用修正后的多項(xiàng)式再計(jì)算出所需的點(diǎn)值與。而時(shí)間離散則使用三階Runge-Kutta時(shí)間離散，式（5）實(shí)際上就是其中的第一步。由于其后兩步都是凸組合的形式，因此格式仍然具有正保護(hù)性質(zhì)。

對(duì)于二維情況，不妨設(shè)區(qū)域被劃分為若干矩形網(wǎng)格單元Iij。在每個(gè)小矩形網(wǎng)格單元中，取如下點(diǎn)集：

與一維情況類似，若在tn時(shí)刻，?α,β,i,j，有。則在CFL條件：

可得tn+1時(shí)刻的。其中qij(x,y)為單元Iij上 各守恒量的k次近似多項(xiàng)式向量(ρij(x,y),mij(x,y),ni j(x,y),Eij(x,y))T。因此，二維情形下正保護(hù)限制器的構(gòu)造及修正方式也與一維情形均基本相同。

2.2 正保護(hù)界面處理方法

本文采用rGFM方法來重定義界面兩側(cè)的流體狀態(tài)。使用該方法需要在界面處求解黎曼問題，因此也需要一種具有正保護(hù)功能的算法。

首先考慮一維情況。假設(shè)界面位于單元j與j+1之間，如圖1所示。

圖1 一維界面邊界條件定義

于是構(gòu)造如下黎曼問題：

對(duì)于求解黎曼問題，已有許多成熟的求解器。本文采用基于雙激波近似的HLLC求解器，如圖2所示。

圖2 雙激波近似黎曼問題

圖中的u?是介質(zhì)界面的速度，sL與sR分別是左右激波的波速，是該黎曼問題的解。以左激波為例，由于流體在激波兩側(cè)的狀態(tài)滿足Rankine-Hugoniot條件，則有

其中

而對(duì)于右側(cè)激波則可以得到類似的式子。又由于在介質(zhì)界面處的速度和壓力連續(xù)，即有。于是由上述等式可得黎曼問題的近似解：

可以證明，當(dāng)取

sL=min(uL?cL,?),sR=max(uR+cR,+),

而對(duì)于二維多介質(zhì)問題，則需要沿介質(zhì)界面的法向構(gòu)造并求解黎曼問題[9]。由于此時(shí)得到的黎曼問題仍是一維形式，因此可以使用上面給出的正保護(hù)方法進(jìn)行求解。與一維問題相比，二維問題中介質(zhì)界面的追蹤要復(fù)雜的多，本文采用Level Set方法[11-12]對(duì)界面進(jìn)行追蹤。

3 數(shù)值實(shí)驗(yàn)

下面將通過幾個(gè)數(shù)值算例來驗(yàn)證算法的有效性。這些算例由于初始條件比較極端，因此會(huì)出現(xiàn)低密度低內(nèi)能區(qū)域，傳統(tǒng)高精度方法容易算得負(fù)值，從而導(dǎo)致程序崩潰。

算例1 氣?氣雙稀疏波問題

該問題的計(jì)算區(qū)域?yàn)閇?1,1]，介質(zhì)界面位于區(qū)域中心。其左側(cè)流體的初始狀態(tài)為：

右側(cè)的初始狀態(tài)為：

在該問題中，介質(zhì)界面兩側(cè)的流體同時(shí)向兩側(cè)反向運(yùn)動(dòng)，從而產(chǎn)生兩個(gè)強(qiáng)稀疏波，使得中心區(qū)域的壓力大幅下降，幾乎形成真空。圖3給出了網(wǎng)格數(shù)為300，t=0.6時(shí)的數(shù)值結(jié)果。

圖3 t=0.6時(shí)的數(shù)值解

算例2 氣?液激波管問題

這是一個(gè)氣?液多介質(zhì)流問題，計(jì)算區(qū)域?yàn)閇0,1]，介質(zhì)界面位于x=0.3處。界面左邊氣體的初始狀態(tài)為：

右邊液體的狀態(tài)為：

該問題是一個(gè)極端條件問題，兩側(cè)流體的初始?jí)毫與狀態(tài)方程中的參數(shù)B相差非常大。隨著時(shí)間推進(jìn)，將形成一個(gè)向左的激波和一個(gè)向右的稀疏波，圖4給出了t=0.000 24時(shí)刻，網(wǎng)格數(shù)取為300的數(shù)值結(jié)果?？梢钥吹剑词箖煞N流體的物理特性及初始條件都相差較大，本文給出的正保護(hù)算法依然能準(zhǔn)確的捕捉到各種間斷，并且介質(zhì)界面非常清晰。

圖4 t=0.000 24時(shí)的數(shù)值解

算例3 高壓氣泡爆炸問題

這是一個(gè)二維氣?氣多介質(zhì)流問題。計(jì)算區(qū)域?yàn)閇?1,1]×[?1,1]，其中充滿靜止的理想氣體，中心處有一個(gè)半徑為0.3的靜止高壓氣泡。氣泡內(nèi)部的流體初始狀態(tài)為：

外部流體的初始狀態(tài)為：

該問題中，氣泡爆炸會(huì)產(chǎn)生一個(gè)向內(nèi)的稀疏波，使得內(nèi)部的密度與壓力值急劇減小。圖5、圖6給出了網(wǎng)格數(shù)為200×200，t=0.000 7時(shí)刻的數(shù)值結(jié)果?？梢钥吹?，本文算法保證了結(jié)果的正確性，且能準(zhǔn)確的追蹤到介質(zhì)界面。

圖5 t=0.000 7時(shí)的等值線圖

圖6 沿中軸線y=0的數(shù)值解

4 結(jié) 論

針對(duì)一些極端條件下的多介質(zhì)流問題，會(huì)因算法等誤差導(dǎo)致出現(xiàn)不合理負(fù)值的情況，本文給出了一種具有正保護(hù)功能的求解方法。該方法首先使用具有正保護(hù)功能的雙激波近似黎曼問題求解器求解黎曼問題，并用得到的近似解對(duì)界面附近的真實(shí)和虛擬流體狀態(tài)進(jìn)行重定義。接著在分別對(duì)各單介質(zhì)流體進(jìn)行求解時(shí)，將針對(duì)理想氣體的正保護(hù)WENO算法推廣運(yùn)用于液體，從而得到了一種高精度正保護(hù)數(shù)值算法。最后通過數(shù)值算例驗(yàn)證了該算法能有效保正結(jié)果的正確性，且分辨率高，能準(zhǔn)確捕捉到各類間斷。