一種新型的高分辨率通量差分裂格式

胡立軍， 趙昆磊， 袁海專

(1.衡陽師范學(xué)院 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，衡陽 421002；2.中國科學(xué)院數(shù)學(xué)與系統(tǒng)科學(xué)研究院 計(jì)算數(shù)學(xué)與科學(xué)工程計(jì)算研究所，北京 100190；3.湘潭大學(xué) 數(shù)學(xué)與計(jì)算科學(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格式計(jì)算代價小并且能夠精確捕捉各類間斷。但是，由于Roe格式不滿足熵條件，故在計(jì)算中會產(chǎn)生非物理的膨脹激波。在計(jì)算強(qiáng)激波問題時，低耗散的Roe，Osher，HLLEM 和HLLC格式還會遭遇紅玉(carbuncle)現(xiàn)象的困擾，這極大地影響了對于高超聲速流動問題的精確模擬。消除低耗散通量差分裂格式激波異?，F(xiàn)象的一些工作主要包括，與耗散格式進(jìn)行混合[6-8]、對Roe格式進(jìn)行熵修正[9]、構(gòu)造旋轉(zhuǎn)格式[10]、反擴(kuò)散控制方法[11]、多維耗散方法[12]、添加人工粘性[13]以及增加剪切粘性[14]等。這些方法通過局部增加格式的數(shù)值耗散成功地抑制了激波不穩(wěn)定現(xiàn)象，但是在計(jì)算復(fù)雜流動問題時，不合適的數(shù)值耗散可能會影響接觸面和剪切層的精度。

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

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

2 預(yù)備知識

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

(1)

式中

(2)

(3)

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

(4)

(5)

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

(6)

(7)

(8)

(9)

Roe格式能夠精確捕捉各類間斷，但是不滿足熵條件，在計(jì)算中會產(chǎn)生非物理的膨脹激波。在計(jì)算強(qiáng)激波問題時還會遭遇紅玉現(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)的通量差分裂格式來計(jì)算壓力數(shù)值通量,

(17)

(18)

(19)

3.2.2 對流數(shù)值通量

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

(20)

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

(21)

從式(13)可知，對流子系統(tǒng)沒有一組完備的線性無關(guān)的特征向量，因此無法直接將ΔU分解成對流子系統(tǒng)特征向量組的線性組合。采用文獻(xiàn)[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í)驗(yàn)表明，耗散控制參數(shù)β取1.1～2.0，都可以得到理想的結(jié)果，β越大耗散越小。后面的數(shù)值實(shí)驗(yàn)β均取1.5。

(26)

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

(27)

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

3.3 穩(wěn)定性分析

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

將計(jì)算區(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的特征值的實(shí)部。

圖1展示了穩(wěn)定性矩陣S的特征值在復(fù)平面上的分布?？梢钥闯觯疚臉?gòu)造的FDS -ZB -BVD格式是穩(wěn)定的，而Roe格式是不穩(wěn)定的。在計(jì)算強(qiáng)激波問題時，F(xiàn)DS -ZB -BVD格式會比Roe格式表現(xiàn)出更好的魯棒性。第4節(jié)的數(shù)值實(shí)驗(yàn)也證明了這一點(diǎn)。

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

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

4.1 一維算例

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

4.1.1 孤立的接觸間斷

計(jì)算馬赫數(shù)M0=0.1的孤立接觸間斷。將計(jì)算區(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問題

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

(37)

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

圖3 Lax問題的密度分布

4.1.3 慢行激波問題

將計(jì)算區(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 一維爆轟波問題

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

(39)

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

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

將計(jì)算區(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 二維算例

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

4.2.1 二維黎曼問題

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

(41)

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

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

4.2.2 移動接觸面問題

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

(42)

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

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

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

將計(jì)算區(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)在計(jì)算強(qiáng)激波問題時會遭遇紅玉現(xiàn)象[11-14]。計(jì)算經(jīng)典的奇偶失聯(lián)問題來檢驗(yàn)FDS -ZB -BVD格式的魯棒性。初始時位于x=5處馬赫數(shù)為10的平面激波沿x-方向傳播。將計(jì)算區(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格式有效地抑制了激波異常現(xiàn)象，得到了清晰和準(zhǔn)確的激波面，表現(xiàn)出了更好的魯棒性。

4.2.5 二維Sedov爆轟波問題

該問題涉及大壓力比和強(qiáng)激波，常用來檢驗(yàn)低耗散數(shù)值格式的魯棒性。將計(jì)算區(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時的壓力等值線?？梢钥闯?，F(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)的通量差分裂格式來進(jìn)行計(jì)算。對于弱雙曲的對流子系統(tǒng)，通過添加廣義特征向量構(gòu)造通量差分裂格式進(jìn)行求解。為了提高接觸間斷的分辨率，利用界面變差下降(BVD)算法來減少對流通量耗散項(xiàng)中的密度差。激波穩(wěn)定性分析表明新格式可以有效抑制數(shù)值誤差的增長。一系列數(shù)值實(shí)驗(yàn)證明了本文構(gòu)造的新型通量差分裂格式比傳統(tǒng)的Roe格式具有更高的分辨率和更好的魯棒性。因此，本文方法是一種精確和魯棒的數(shù)值方法，可以廣泛應(yīng)用于可壓縮流的數(shù)值模擬。