胡立軍， 杜玉龍

(1.衡陽師范學(xué)院 數(shù)學(xué)與統(tǒng)計學(xué)院，衡陽 421002；2.北京航空航天大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，北京 100191)

1 引 言

基于Godunov型數(shù)值格式[1]的有限體積法由于其良好的守恒性和網(wǎng)格適應(yīng)性，已經(jīng)成為求解雙曲守恒律系統(tǒng)最具代表性的數(shù)值方法。通過求解局部黎曼問題來得到網(wǎng)格界面的數(shù)值通量是實現(xiàn)有限體積方法的關(guān)鍵步驟。從物理角度來看，采用迭代方法求解黎曼問題的精確解似乎是最合理的，但是由于其效率低且實現(xiàn)難度大，在實際問題的計算中精確黎曼求解器很少使用。因此構(gòu)造性能良好的近似黎曼求解器一直都是計算流體力學(xué)研究的熱點問題。

過去幾十年，研究人員構(gòu)造了許多不同性能的近似黎曼求解器。根據(jù)捕捉接觸間斷的能力，可以將其分為兩類。非全波求解器，如Rusanov格式[2]、HLLE格式[3]和HLL-CPS格式[4]，過高的耗散行為在計算中不能精確分辨接觸波或者剪切波；全波求解器，如Roe格式[5]、Osher格式[6]和HLLC格式[7]，在計算中能夠精確捕捉接觸間斷和剪切波。但是，在計算強激波問題時，這些低耗散的求解器會遭遇嚴(yán)重的不穩(wěn)定現(xiàn)象，這大大限制了它們在高超聲速流動問題中的應(yīng)用。

研究人員嘗試在保留全波求解器精確分辨接觸間斷優(yōu)點的同時來消除它們的激波不穩(wěn)定性，其中最流行的方法是根據(jù)當(dāng)?shù)亓鲌鲈诤纳⒏袷胶偷秃纳⒏袷街g進(jìn)行切換的混合方法[8-12]。盡管混合格式可以成功地抑制激波失穩(wěn)現(xiàn)象，但是對于復(fù)雜的流動問題，不恰當(dāng)?shù)拈_關(guān)函數(shù)可能會影響接觸間斷的分辨率，并且兩種格式間的突然切換也可能會影響到格式的收斂速度[13]。此外，增加多維耗散也是抑制全波黎曼求解器激波失穩(wěn)現(xiàn)象的一種常用方法，主要通過增加格式的剪切粘性[14,15]、構(gòu)造旋轉(zhuǎn)黎曼求解器[16,17]或者法向速度重建技術(shù)[18]來實現(xiàn)。旋轉(zhuǎn)格式在每個網(wǎng)格界面處需要計算兩次數(shù)值通量，因此計算效率低。而單純地增加剪切粘性在某些情形下并不能完全消除激波異常現(xiàn)象[19]。最近，Chen等[20]通過在動量通量中引入剪切粘性構(gòu)造了一種具有良好激波穩(wěn)定性的HLLC+格式。

此外提高非全波黎曼求解器的接觸捕捉精度也是構(gòu)造性能良好的數(shù)值格式的一種策略。本文利用雙曲正切函數(shù)[21]來重構(gòu)界面兩側(cè)的密度值，并且結(jié)合邊界變差遞減算法[22]來減小單波Rusanov格式耗散項中的密度差，從而提高格式分辨接觸間斷的能力。數(shù)值結(jié)果表明,本文構(gòu)造的黎曼求解器比全波的HLLC格式具有更高的接觸分辨率和更好的激波穩(wěn)定性。

2 背景知識

守恒形式的二維歐拉方程組為

(1)

式中

(2)

(3)

本文的比熱比γ取1.4。

采用有限體積法對式(1)進(jìn)行空間離散，可以得到半離散方程為

(4)

式中Ui,j為守恒向量U的單元平均值，F(xiàn)i +1/2,j和Gi,j +1/2為數(shù)值通量。式(4)可以改寫成ODE系統(tǒng)為

dU/dt=L(U)

(5)

2.1 時間離散

采用優(yōu)化的三階TVD龍格-庫塔格式[23]來進(jìn)行時間離散,

U(1)=Un+ΔtL(Un)

(6)

式中Δt為時間步長。

2.2 五階WENO格式

采用五階WENO格式[24]重構(gòu)得到界面左右兩側(cè)的狀態(tài)值為

(7)

(8)

非線性權(quán)重系數(shù)定義為

(9)

(10)

式(7)的線性加權(quán)系數(shù)為

(11)

而式(8)的線性加權(quán)系數(shù)為

(12)

為防止分母為0，取ε=10-6。

利用WENO重構(gòu)得到界面左右兩側(cè)的狀態(tài)值，然后利用近似黎曼求解器求解局部黎曼問題，從而得到每個網(wǎng)格界面處的數(shù)值通量。

2.3 Rusanov求解器

最簡單的單波Rusanov黎曼求解器[2]的數(shù)值通量函數(shù)為

(13)

波速S定義為

S=max(|uL-aL|,|uR-aR|,|uL+aL|,|uR+aR|)

(14)

2.4 HLLC求解器

全波的HLLC黎曼求解器的數(shù)值通量函數(shù)為

(15)

式中

(16)

波速估計為

(17)

3 一種精確分辨接觸間斷的Rusanov黎曼求解器

Rusanov格式的通量函數(shù)式(13)可以改寫為

(18)

式中右端第一項為中心差分項，D1/2為耗散項，其表達(dá)式為

(19)

(20)

(21)

為了計算方便，本文直接給出網(wǎng)格界面左右兩側(cè)密度重構(gòu)值的計算公式為

(22)

式中

(23)

式中β為控制密度跳躍大小的參數(shù)，本文取β=1.6。

(24)

將由式(24)確定的界面兩側(cè)的密度值代入式(19)計算耗散項D1/2，從而得到一種低耗散的Rusanov格式(命名為LD -Rusanov，Low Diffusion Rusanov)。此外,本文為了提高非線性波的分辨率，采用Roe平均值來計算波速，

(25)

在計算多維問題時，采用算子分裂方法來逐維計算數(shù)值通量。因此本文所述的一維算法可以直接應(yīng)用于每個坐標(biāo)方向。雙曲正切重構(gòu)作為一種代數(shù)方法，在進(jìn)行多維計算時不會涉及到任何的幾何重構(gòu)，因此實施起來比較簡單。此外，與HLLC格式通過修改HLL格式的波系結(jié)構(gòu)來提高接觸波的分辨率不同的是，LD-Rusanov格式采用代數(shù)方法而不是修改原Rusanov格式的波系結(jié)構(gòu)來提高接觸間斷的捕捉能力，因此其很好地保留了原Rusanov格式的激波穩(wěn)定性。接下來，通過一系列的數(shù)值實驗來證明LD-Rusanov格式的高分辨率和強魯棒性。

4 數(shù)值結(jié)果

通過一些典型的數(shù)值算例來驗證LD -Rusanov格式的表現(xiàn)。

4.1 孤立的接觸間斷

首先考慮一個孤立的慢行接觸波，計算區(qū)域為[0,1]，其初始條件為

(26)

計算中使用網(wǎng)格數(shù)為100的均勻網(wǎng)格。從 圖1 可以看出，本文構(gòu)造的LD -Rusanov格式不僅可以更加精確地捕捉該接觸間斷，而且還消除了原Rusanov格式在間斷附近出現(xiàn)的過沖。圖2展示了取不同β值的LD -Rusanov格式計算該算例時的計算結(jié)果。可以看出，β=1.3,1.6和1,9都可以得到精確的結(jié)果。為了定量分析不同β值對于結(jié)果的影響，采用式(27)來計算間斷的厚度，

(27)

圖1 孤立的接觸間斷在t =3時的密度分布

圖2 取不同β值的LD -Rusanov格式的密度分布

由表1可知，β值越大，間斷的厚度越小，其分辨率也會越高，但是過大的β會卷起平行于速度方向的界面[21]。本文發(fā)現(xiàn)1.2～2.0的β值可以得到令人滿意的結(jié)果。

4.2 一維爆轟波問題

計算區(qū)域為[0,1]，初始條件為

(28)

計算中使用網(wǎng)格數(shù)為200的均勻網(wǎng)格。該算例的解由兩個相互作用的爆轟波形成的復(fù)雜間斷結(jié)構(gòu)構(gòu)成。從圖3可以看出，相比原始的Rusanov格式和全波的HLLC格式，本文構(gòu)造的LD -Rusanov格式能夠更加精確地捕捉各個波系。

圖3 一維爆轟波問題在t =0.038時的密度分布

4.3 Titarev-Toro問題

考慮Titarev等[24]提出的激波-熵波相互作用問題。計算區(qū)域[-5,5]均勻劃分成500個網(wǎng)格，初始條件為

(29)

該算例涉及到高頻振蕩的正弦波和激波的相互作用。從圖4可以看出，Rusanov格式的計算結(jié)果具有較大的數(shù)值耗散，而本文構(gòu)造的LD -Rusanov格式得到了更加精確的解，在捕捉高頻波時比全波的HLLC格式具有更高的分辨率。

圖4 Titarev-Toro問題在t =5時的密度分布

4.4 二維爆炸問題

二維爆炸問題是一維Sod問題的推廣，計算該算例來檢驗格式捕捉不同波系的能力。計算區(qū)域[-1,1]×[-1,1]均勻劃分成101×101的正方形網(wǎng)格，初始條件為

(30)

計算時間為t=0.25。此時，該算例的解由一個向外運動的圓形激波和接觸面以及向中心運動的稀疏波構(gòu)成。圖5展示了時間t=0.25時沿徑向y=0的密度分布?？梢钥闯?，本文構(gòu)造的 LD -Rusanov 格式能夠精確地捕捉各個波系，特別是對接觸間斷的分辨率明顯優(yōu)于原始的Rusanov格式和HLLC格式。

圖5 二維爆炸問題在t =0.25時的密度分布

4.5 移動的圓形接觸面問題

計算區(qū)域[0,1]×[0,1]劃分成100×100的均勻網(wǎng)格，其初始分布為

(31)

該算例的解由一個以勻速(u,v)=(1,1)運動的圓形接觸面構(gòu)成。在t=0.3時刻該問題的精確解為

(32)

從圖6可以看出，相比于原始的Rusanov格式和全波的HLLC格式，本文構(gòu)造的LD -Rusanov格式對于該圓形接觸面具有更高的分辨率。

圖6 移動的圓形接觸面問題的計算結(jié)果

4.6 奇偶失聯(lián)問題

本文數(shù)值實驗表明，在捕捉接觸間斷時LD -Rusanov格式比HLLC格式具有更高的分辨率。文獻(xiàn)[8-20]曾報道精確分辨接觸間斷的數(shù)值格式(如HLLC格式)在計算多維強激波問題時會出現(xiàn)嚴(yán)重的不穩(wěn)定現(xiàn)象。接下來，計算幾個典型的強激波問題來檢驗LD -Rusanov格式的魯棒性。

首先,計算Quirk[8]提出的奇偶失聯(lián)算例來評估格式的激波穩(wěn)定性。馬赫數(shù)為20的平面運動激波從左向右傳播，其初始位置為x=5。右側(cè)的波前狀態(tài)為(ρ,u,v,p)R=(1.4,0,0,1)，左側(cè)的波后狀態(tài)由Rankine-Hugoniot關(guān)系式計算得到。計算區(qū)域[0,600]×[0,20]分割成600×20的矩形網(wǎng)格，且y-方向的網(wǎng)格中心線上存在±10-3的奇偶小擾動。圖7展示了t=20時的密度等值圖，為了全面比較格式的激波穩(wěn)定性，本文還展示了兩種穩(wěn)定版本的HLLC格式(HLLCM[13]和HLLC+[20])的計算結(jié)果?？梢钥闯?，HLLC格式表現(xiàn)出了明顯的激波失穩(wěn)現(xiàn)象，而LD -Rusanov格式和另外兩種穩(wěn)定版本的HLLC格式均得到了穩(wěn)定的激波面。

圖7 奇偶失聯(lián)問題的密度等值圖

4.7 穩(wěn)態(tài)斜激波問題

接下來考慮穩(wěn)態(tài)斜激波問題[26]。計算區(qū)域[0,50]×[0,30]劃分成50×30的正方形網(wǎng)格，沿著直線y=2(x-12)有一個穩(wěn)態(tài)激波，其波前狀態(tài)為(ρ,u,v,p)L=(1,447.26,-223.5,3644.31)，波后狀態(tài)(ρ,u,v,p)R=(5.444,82.15,-41.05,207725.94)。圖8展示了時間迭代1000步的計算結(jié)果?？梢钥闯?，全波的HLLC格式出現(xiàn)了明顯的激波不穩(wěn)定現(xiàn)象，而LD-Rusanov格式保留了原始Rusanov格式的穩(wěn)定性，得到了清晰的激波面。

圖8 穩(wěn)態(tài)斜激波問題的密度等值線

4.8 圓柱繞流問題

計算馬赫數(shù)為20的無粘圓柱繞流問題來評估不同數(shù)值格式在高超聲速流動問題中的激波穩(wěn)定性。該算例具體的計算區(qū)域和初邊值條件的描述參考文獻(xiàn)[12]。本文采用40×80的非規(guī)則結(jié)構(gòu)四邊形網(wǎng)格來劃分計算區(qū)域。圖9展示了t=4時的密度等值圖?？梢钥闯?，HLLC格式表現(xiàn)出明顯的激波失穩(wěn)現(xiàn)象，出現(xiàn)了嚴(yán)重的紅玉現(xiàn)象，而LD -Rusanov格式和另外兩種穩(wěn)定版本的HLLC格式均得到了穩(wěn)定的數(shù)值解。

圖9 無粘圓柱繞流問題的密度等值圖

4.9 正激波繞射問題

計算馬赫數(shù)為5.09的正激波繞射問題來驗證格式的魯棒性。采用200×200的均勻網(wǎng)格來劃分計算區(qū)域[0,1]×[0,1]，其中角點位于(x,y)=(0.05,0.6)處。初始時正激波右側(cè)靜止流體的密度為1.4，壓力為1。角點上方左側(cè)的入口邊界條件可以由Rankine -Hugoniot關(guān)系式計算得到，區(qū)域底部和角點下方左側(cè)設(shè)置為反射邊界條件，在右邊界所有變量的梯度設(shè)置為零，區(qū)域頂部的邊界條件會隨著激波的運動而調(diào)整。圖10展示了t=0.15 時的密度等值圖?？梢钥闯觯琀LLC格式在正激波的上方出現(xiàn)了明顯的激波失穩(wěn)現(xiàn)象，而LD -Rusanov格式和另外兩種穩(wěn)定版本的HLLC格式均得到了穩(wěn)定的正激波。值得注意的是，在該算例中本文構(gòu)造的LD -Rusanov格式表現(xiàn)出了最好的穩(wěn)定性。

圖10 正激波繞射問題的密度等值圖

5 結(jié) 論

基于單波的Rusanov求解器構(gòu)造了一種低耗散的LD -Rusanov黎曼求解器。使用雙曲正切函數(shù)和五階WENO格式來重構(gòu)界面兩側(cè)的密度值。利用邊界變差下降算法從密度重構(gòu)值的不同組合中挑選出使界面兩側(cè)密度差最小的密度值，從而最大限度地減小Rusanov求解器耗散項中的密度差，進(jìn)而提高格式對于接觸間斷的分辨率。全波的HLLC格式通過修改HLL格式的波系結(jié)構(gòu)來提高接觸間斷的分辨率，因此無法保留HLL格式的激波穩(wěn)定性。本文構(gòu)造的LD -Rusanov格式使用一種代數(shù)方法來提高接觸間斷的捕捉能力，并沒有修改原Rusanov求解器的波系結(jié)構(gòu)，因此很好地保留了Rusanov的激波穩(wěn)定性。一系列數(shù)值實驗表明LD -Rusanov格式比全波的HLLC格式具有更好的接觸間斷捕捉能力和魯棒性。本文構(gòu)造的低耗散數(shù)值格式展現(xiàn)出了良好的應(yīng)用前景，因此將其應(yīng)用于復(fù)雜流動問題(如可壓縮湍流和化學(xué)反應(yīng)流)的數(shù)值模擬值得未來進(jìn)一步的研究。