胡立軍， 趙昆磊， 袁海專

(1.衡陽師范學(xué)院 數(shù)學(xué)與統(tǒng)計學(xué)院，衡陽 421002；2.中國科學(xué)院數(shù)學(xué)與系統(tǒng)科學(xué)研究院 計算數(shù)學(xué)與科學(xué)工程計算研究所，北京 100190；3.湘潭大學(xué) 數(shù)學(xué)與計算科學(xué)學(xué)院,湘潭 411105)

1 引 言

數(shù)值模擬高馬赫數(shù)流動問題需要精確魯棒的數(shù)值格式。無論在理論研究還是在工程應(yīng)用中，迎風(fēng)格式都是最受歡迎的激波捕捉方法。根據(jù)反映控制方程組物理特性的不同方式，通?？蓪⑵浞譃橥坎罘至?FDS)格式和通量向量分裂(FVS)格式。流行的通量差分裂格式主要包括Roe格式[1]、Osher格式[2]和HLL-類格式[3-5]。其中，利用系統(tǒng)的特征值和特征向量構(gòu)造的Roe格式計算代價小并且能夠精確捕捉各類間斷。但是，由于Roe格式不滿足熵條件，故在計算中會產(chǎn)生非物理的膨脹激波。在計算強激波問題時，低耗散的Roe，Osher，HLLEM 和HLLC格式還會遭遇紅玉(carbuncle)現(xiàn)象的困擾，這極大地影響了對于高超聲速流動問題的精確模擬。消除低耗散通量差分裂格式激波異常現(xiàn)象的一些工作主要包括，與耗散格式進行混合[6-8]、對Roe格式進行熵修正[9]、構(gòu)造旋轉(zhuǎn)格式[10]、反擴散控制方法[11]、多維耗散方法[12]、添加人工粘性[13]以及增加剪切粘性[14]等。這些方法通過局部增加格式的數(shù)值耗散成功地抑制了激波不穩(wěn)定現(xiàn)象，但是在計算復(fù)雜流動問題時，不合適的數(shù)值耗散可能會影響接觸面和剪切層的精度。

以Steger-Warming格式[15]和van Leer格式[16]為代表的通量向量分裂方法在計算中具有很好的魯棒性，但是其主要的缺點是不能精確分辨接觸間斷和剪切波。此外，基于對流-壓力通量分裂方法構(gòu)造的混合迎風(fēng)格式也吸引了許多研究人員的關(guān)注。Liou等[17]將動量方程中的壓力項分裂出來，并且構(gòu)造了一種可以精確分辨接觸間斷的對流迎風(fēng)分裂方法(AUSM)。目前，AUSM格式已經(jīng)發(fā)展成為具有廣泛應(yīng)用前景的AUSM-類格式[18-20]。Zha等[21]將動量方程中的壓力項以及能量方程中速度與壓力的乘積項分裂出來進而構(gòu)造了一種混合迎風(fēng)格式。此外，Toro等[22]提出了一種新的對流-壓力通量分裂方法，其特點是對流通量中完全不含壓力項。基于對流-壓力通量分裂方法構(gòu)造的混合迎風(fēng)格式形式簡單并且能夠精確捕捉接觸間斷，但是在計算某些問題時仍然會出現(xiàn)紅玉現(xiàn)象和激波波后振蕩[23]。

本文基于Zha-Bilgen對流-壓力通量分裂方法來構(gòu)造一種精確魯棒的新型通量差分裂格式。壓力子系統(tǒng)具有一組完備的線性無關(guān)的特征向量，因此可以構(gòu)造傳統(tǒng)的通量差分裂格式來進行計算。弱雙曲的對流子系統(tǒng)采用文獻[24]添加廣義特征向量的方法構(gòu)造通量差分裂格式來進行計算。為了提高接觸間斷的分辨率，利用界面變差下降(BVD)算法[25]來減小對流通量耗散項的密度差。穩(wěn)定性分析和數(shù)值實驗證明了新的通量差分裂格式比Roe格式具有更好的魯棒性和更高的分辨率。

2 預(yù)備知識

二維無粘可壓縮流的控制方程組可以表示為

(1)

式中

(2)

(3)

除非特別說明，否則比熱比γ取1.4。將計算區(qū)域劃分成結(jié)構(gòu)化的四邊形網(wǎng)格并且采用有限體積方法對式(1)進行空間離散,

(4)

(5)

Roe格式[1]數(shù)值通量的表達(dá)式為

(6)

(7)

(8)

(9)

Roe格式能夠精確捕捉各類間斷，但是不滿足熵條件，在計算中會產(chǎn)生非物理的膨脹激波。在計算強激波問題時還會遭遇紅玉現(xiàn)象。

3 一種新型的通量差分裂格式

基于Zha-Bilgen分裂方法來構(gòu)造一種精確和魯棒的通量差分裂格式。

3.1 Zha-Bilgen分裂

Zha等[21]將歐拉方程的通量分裂成對流通量和壓力通量兩部分。

(10)

對流通量Fc的雅可比矩陣為

(11)

Ac的特征值為

(12)

其線性無關(guān)的特征向量為

(13)

壓力通量Fp的雅可比矩陣為

(14)

Ap的特征值為

(15)

其線性無關(guān)的特征向量為

(16)

從式(12，15)可以看出，對流子系統(tǒng)的四個特征值都等于流體在該方向的速度，從而表現(xiàn)出很好的對流特性。而壓力子系統(tǒng)的特征值等于聲速的常數(shù)倍，很好地反映了聲速波的特性。

3.2 FDS -ZB -BVD格式

3.2.1 壓力數(shù)值通量

從式(16)可以看出，壓力子系統(tǒng)有一組完備的線性無關(guān)的特征向量，因此可以構(gòu)造傳統(tǒng)的通量差分裂格式來計算壓力數(shù)值通量,

(17)

(18)

(19)

3.2.2 對流數(shù)值通量

對流子系統(tǒng)滿足微分關(guān)系式,

(20)

將式(20)改寫成有限差分形式,

(21)

從式(13)可知，對流子系統(tǒng)沒有一組完備的線性無關(guān)的特征向量，因此無法直接將ΔU分解成對流子系統(tǒng)特征向量組的線性組合。采用文獻[24]的約旦標(biāo)準(zhǔn)型理論來添加廣義特征向量，并且注意到對流子系統(tǒng)有四個相等的特征值u，可以得到：

(22)

因此，對流子系統(tǒng)的數(shù)值通量可以表示為

(23)

為了減少數(shù)值耗散，利用THINC函數(shù)來重構(gòu)界面兩側(cè)的密度[25]，

(24)

式中ρmin=min(ρi -1,ρi +1), Δρ=|ρi +1-ρi -1|

A=[B-cosh(β)]/sinh(β)

θ=sgn(ρi +1-ρi -1),ε=10-20

(25)

數(shù)值實驗表明，耗散控制參數(shù)β取1.1～2.0，都可以得到理想的結(jié)果，β越大耗散越小。后面的數(shù)值實驗β均取1.5。

(26)

因此，基于Zha-Bilgen分裂方法并且結(jié)合BVD算法構(gòu)造的通量差分裂格式(命名為FDS -ZB -BVD)數(shù)值通量的表達(dá)式為

(27)

本文將通過穩(wěn)定性分析和數(shù)值實驗來證明該格式比Roe格式具有更好的魯棒性和更高的分辨率。

3.3 穩(wěn)定性分析

低耗散的數(shù)值格式在計算強激波問題時會遭遇激波不穩(wěn)定現(xiàn)象。采用矩陣穩(wěn)定性分析方法[26]對Roe格式和FDS -ZB -BVD格式進行穩(wěn)定性分析。

將計算區(qū)域[0,1]×[0,1]劃分成20×20的笛卡爾網(wǎng)格，馬赫數(shù)M0=7的穩(wěn)態(tài)駐激波位于第10列和第11列網(wǎng)格單元共有的網(wǎng)格界面上。激波的上游狀態(tài)為

(28)

下游狀態(tài)通過Rankine -Hugoniot關(guān)系式來得到

(29)

本文研究不同格式的數(shù)值誤差隨著時間的發(fā)展趨勢。整個區(qū)域的流場可以分解為

(30)

(31)

將線性化的通量函數(shù)代入式(4)，可以得到

(32)

將式(32)改寫成矩陣形式，

(33)

式中q=20×20為網(wǎng)格單元總數(shù)，S為穩(wěn)定性矩陣。求解式(33)，得到數(shù)值誤差的演化方程為

(34)

誤差發(fā)展有界的條件為

max{Re[λ(S)]}≤0

(35)

式中Re[λ(S)]表示穩(wěn)定性矩陣S的特征值的實部。

圖1展示了穩(wěn)定性矩陣S的特征值在復(fù)平面上的分布。可以看出，本文構(gòu)造的FDS -ZB -BVD格式是穩(wěn)定的，而Roe格式是不穩(wěn)定的。在計算強激波問題時，F(xiàn)DS -ZB -BVD格式會比Roe格式表現(xiàn)出更好的魯棒性。第4節(jié)的數(shù)值實驗也證明了這一點。

4 數(shù)值結(jié)果和分析

通過數(shù)值實驗來檢驗FDS -ZB -BVD格式的精度和魯棒性。采用二階MUSCL格式獲得空間的二階精度，時間離散采用二階Runge -Kutta格式。

4.1 一維算例

計算幾個經(jīng)典的算例來比較FDS -ZB -BVD格式和Roe格式捕捉間斷的能力。

4.1.1 孤立的接觸間斷

計算馬赫數(shù)M0=0.1的孤立接觸間斷。將計算區(qū)域[0,1]均勻劃分成100個網(wǎng)格，初始條件為

(36)

圖2展示了t=2.0時的密度分布?？梢钥闯?，F(xiàn)DS -ZB -BVD格式可以更加精確地捕捉該接觸間斷。

圖1 穩(wěn)定性矩陣S的特征值在復(fù)平面上的分布

圖2 孤立的接觸間斷的密度分布

4.1.2 Lax問題

將計算區(qū)域[0,1]均勻劃分成100個網(wǎng)格，初始條件為

(37)

該算例由一個左行的稀疏波、一個右行的接觸間斷和強激波構(gòu)成。圖3展示了t=0.15時的密度分布?？梢钥闯?，F(xiàn)DS -ZB -BVD格式在接觸間斷處比Roe格式的耗散更小，捕捉激波和稀疏波時與Roe格式具有相同的分辨率。

圖3 Lax問題的密度分布

4.1.3 慢行激波問題

將計算區(qū)域[0,1]均勻劃分成500個網(wǎng)格，初始條件為

(38)

該算例由一個孤立的右行慢激波構(gòu)成。圖4 展示了t=1.5時的密度分布。可以看出，Roe格式在波后出現(xiàn)了明顯的振蕩，而FDS -ZB -BVD格式有效地抑制了波后的振蕩，表現(xiàn)出了更好的魯棒性。

圖4 慢行激波問題的密度分布

4.1.4 一維爆轟波問題

將計算區(qū)域[0,1]均勻劃分成500個網(wǎng)格，初始條件為

(39)

該算例包含了復(fù)雜的間斷結(jié)構(gòu)和波系間的相互作用。圖5展示了t=0.38時的密度分布。可以看出，F(xiàn)DS -ZB -BVD格式對于各個波系具有更高的分辨率。

4.1.5 激波-熵波相互作用問題

將計算區(qū)域[-5,5]均勻劃分成1000個網(wǎng)格，初始條件為

(40)

該算例涉及到馬赫數(shù)為3的激波和密度擾動產(chǎn)生的余弦波相互作用，產(chǎn)生了包含光滑和間斷結(jié)構(gòu)的流場。圖6展示了t=1.8時的密度分布。可以看出，F(xiàn)DS -ZB -BVD格式比Roe格式具有更高的分辨率。

圖6 激波-熵波相互作用問題的密度分布

4.2 二維算例

通過算例檢驗FDS -ZB -BVD格式在計算二維問題時的精度和魯棒性。本文僅考慮用結(jié)構(gòu)化四邊形網(wǎng)格來劃分規(guī)則的矩形區(qū)域，對于不規(guī)則的計算區(qū)域，可以使用非結(jié)構(gòu)網(wǎng)格(如三角形網(wǎng)格)進行區(qū)域的劃分，然后利用本文構(gòu)造的數(shù)值格式結(jié)合旋轉(zhuǎn)不變性式(5)計算界面的數(shù)值通量，最后應(yīng)用有限體積方法來離散控制方程組。

4.2.1 二維黎曼問題

將計算區(qū)域[0,1]×[0,1]劃分成512×512的正方形網(wǎng)格，初始條件為

(41)

圖7展示了t=0.23時的密度等值線?？梢钥闯觯現(xiàn)DS -ZB -BVD格式比Roe格式具有更高的分辨率。該算例表明，本文構(gòu)造的新型通量差分裂格式在計算二維問題時，可以有效地改善接觸間斷和相關(guān)流場結(jié)構(gòu)的分辨率。

圖7 二維黎曼問題的密度等值線

4.2.2 移動接觸面問題

將計算區(qū)域[0,1]×[0,1]劃分成200×200的正方形網(wǎng)格，初始條件為

(42)

該算例涉及到以(0.5,0.5)為圓心，0.3為半徑的圓環(huán)接觸面以速度(u,v)=(1,1)勻速運動。圖8展示了t=0.15時的密度等值線，此時圓心的位置為(0.65,0.65)?？梢钥闯?，F(xiàn)DS -ZB -BVD格式大大改善了接觸面的分辨率。

圖8 移動接觸面問題的密度等值線

4.2.3 Rayleigh-Taylor不穩(wěn)定性問題

將計算區(qū)域[0,0.25]×[0,1]劃分成120×480的正方形網(wǎng)格。區(qū)域上下部分兩種流體的初始條件為

(43)

圖9 Rayleigh-Taylor不穩(wěn)定性問題的密度等值線

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

能夠精確捕捉接觸間斷的低耗散數(shù)值格式(如Roe和HLLC)在計算強激波問題時會遭遇紅玉現(xiàn)象[11-14]。計算經(jīng)典的奇偶失聯(lián)問題來檢驗FDS -ZB -BVD格式的魯棒性。初始時位于x=5處馬赫數(shù)為10的平面激波沿x-方向傳播。將計算區(qū)域[0,1500]×[0,20]劃分成1500×20的均勻笛卡爾網(wǎng)格，除了y-方向的網(wǎng)格中心線處存在如下奇偶對稱擾動,

(44)

激波右側(cè)流體的密度為1.4，壓力為1，兩個方向的速度均為0。圖10展示了t=120時的密度等值線。可以看出，Roe格式出現(xiàn)了明顯的紅玉現(xiàn)象，激波面完全破壞，而FDS -ZB -BVD格式有效地抑制了激波異?，F(xiàn)象，得到了清晰和準(zhǔn)確的激波面，表現(xiàn)出了更好的魯棒性。

4.2.5 二維Sedov爆轟波問題

該問題涉及大壓力比和強激波，常用來檢驗低耗散數(shù)值格式的魯棒性。將計算區(qū)域[0,2.4]×[0,2.4]劃分成480×480的正方形網(wǎng)格。初始時，區(qū)域中心高壓區(qū)的壓力為3.5×105，其余地方的壓力為10-10，整個區(qū)域流體的密度為1，兩個方向的速度均為0。圖11展示了t=0.1時的壓力等值線?？梢钥闯觯現(xiàn)DS -ZB -BVD格式消除了Roe格式非物理的紅玉現(xiàn)象，表現(xiàn)出了更好的魯棒性。

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

圖11 二維Sedov爆轟波問題的壓力等值線

5 結(jié) 論

基于Zha-Bilgen對流-壓力通量分裂方法構(gòu)造了一種新型的通量差分裂格式。由于壓力子系統(tǒng)具有一組完備的線性無關(guān)的特征向量，因此構(gòu)造傳統(tǒng)的通量差分裂格式來進行計算。對于弱雙曲的對流子系統(tǒng)，通過添加廣義特征向量構(gòu)造通量差分裂格式進行求解。為了提高接觸間斷的分辨率，利用界面變差下降(BVD)算法來減少對流通量耗散項中的密度差。激波穩(wěn)定性分析表明新格式可以有效抑制數(shù)值誤差的增長。一系列數(shù)值實驗證明了本文構(gòu)造的新型通量差分裂格式比傳統(tǒng)的Roe格式具有更高的分辨率和更好的魯棒性。因此，本文方法是一種精確和魯棒的數(shù)值方法，可以廣泛應(yīng)用于可壓縮流的數(shù)值模擬。