胡立軍， 杜玉龍

(1.衡陽師范學(xué)院 數(shù)學(xué)與統(tǒng)計(jì)學(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ù)值通量是實(shí)現(xiàn)有限體積方法的關(guān)鍵步驟。從物理角度來看，采用迭代方法求解黎曼問題的精確解似乎是最合理的，但是由于其效率低且實(shí)現(xiàn)難度大，在實(shí)際問題的計(jì)算中精確黎曼求解器很少使用。因此構(gòu)造性能良好的近似黎曼求解器一直都是計(jì)算流體力學(xué)研究的熱點(diǎn)問題。

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

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

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

2 背景知識(shí)

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

(1)

式中

(2)

(3)

本文的比熱比γ取1.4。

采用有限體積法對(duì)式(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 時(shí)間離散

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

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

(6)

式中Δt為時(shí)間步長(zhǎng)。

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)值，然后利用近似黎曼求解器求解局部黎曼問題，從而得到每個(gè)網(wǎng)格界面處的數(shù)值通量。

2.3 Rusanov求解器

最簡(jiǎn)單的單波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)

波速估計(jì)為

(17)

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

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

(18)

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

(19)

(20)

(21)

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

(22)

式中

(23)

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

(24)

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

(25)

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

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

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

4.1 孤立的接觸間斷

首先考慮一個(gè)孤立的慢行接觸波，計(jì)算區(qū)域?yàn)閇0,1]，其初始條件為

(26)

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

(27)

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

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

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

4.2 一維爆轟波問題

計(jì)算區(qū)域?yàn)閇0,1]，初始條件為

(28)

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

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

4.3 Titarev-Toro問題

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

(29)

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

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

4.4 二維爆炸問題

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

(30)

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

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

4.5 移動(dòng)的圓形接觸面問題

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

(31)

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

(32)

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

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

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

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

首先,計(jì)算Quirk[8]提出的奇偶失聯(lián)算例來評(píng)估格式的激波穩(wěn)定性。馬赫數(shù)為20的平面運(yùn)動(dòng)激波從左向右傳播，其初始位置為x=5。右側(cè)的波前狀態(tài)為(ρ,u,v,p)R=(1.4,0,0,1)，左側(cè)的波后狀態(tài)由Rankine-Hugoniot關(guān)系式計(jì)算得到。計(jì)算區(qū)域[0,600]×[0,20]分割成600×20的矩形網(wǎng)格，且y-方向的網(wǎng)格中心線上存在±10-3的奇偶小擾動(dòng)。圖7展示了t=20時(shí)的密度等值圖，為了全面比較格式的激波穩(wěn)定性，本文還展示了兩種穩(wěn)定版本的HLLC格式(HLLCM[13]和HLLC+[20])的計(jì)算結(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]。計(jì)算區(qū)域[0,50]×[0,30]劃分成50×30的正方形網(wǎng)格，沿著直線y=2(x-12)有一個(gè)穩(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展示了時(shí)間迭代1000步的計(jì)算結(jié)果。可以看出，全波的HLLC格式出現(xiàn)了明顯的激波不穩(wěn)定現(xiàn)象，而LD-Rusanov格式保留了原始Rusanov格式的穩(wěn)定性，得到了清晰的激波面。

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

4.8 圓柱繞流問題

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

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

4.9 正激波繞射問題

計(jì)算馬赫數(shù)為5.09的正激波繞射問題來驗(yàn)證格式的魯棒性。采用200×200的均勻網(wǎng)格來劃分計(jì)算區(qū)域[0,1]×[0,1]，其中角點(diǎn)位于(x,y)=(0.05,0.6)處。初始時(shí)正激波右側(cè)靜止流體的密度為1.4，壓力為1。角點(diǎn)上方左側(cè)的入口邊界條件可以由Rankine -Hugoniot關(guān)系式計(jì)算得到，區(qū)域底部和角點(diǎn)下方左側(cè)設(shè)置為反射邊界條件，在右邊界所有變量的梯度設(shè)置為零，區(qū)域頂部的邊界條件會(huì)隨著激波的運(yùn)動(dòng)而調(diào)整。圖10展示了t=0.15 時(shí)的密度等值圖?？梢钥闯觯琀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求解器耗散項(xiàng)中的密度差，進(jìn)而提高格式對(duì)于接觸間斷的分辨率。全波的HLLC格式通過修改HLL格式的波系結(jié)構(gòu)來提高接觸間斷的分辨率，因此無法保留HLL格式的激波穩(wěn)定性。本文構(gòu)造的LD -Rusanov格式使用一種代數(shù)方法來提高接觸間斷的捕捉能力，并沒有修改原Rusanov求解器的波系結(jié)構(gòu)，因此很好地保留了Rusanov的激波穩(wěn)定性。一系列數(shù)值實(shí)驗(yàn)表明LD -Rusanov格式比全波的HLLC格式具有更好的接觸間斷捕捉能力和魯棒性。本文構(gòu)造的低耗散數(shù)值格式展現(xiàn)出了良好的應(yīng)用前景，因此將其應(yīng)用于復(fù)雜流動(dòng)問題(如可壓縮湍流和化學(xué)反應(yīng)流)的數(shù)值模擬值得未來進(jìn)一步的研究。