胡立軍， 趙昆磊， 袁海專(zhuān)

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

1 引 言

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

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

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

2 預(yù)備知識(shí)

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

(1)

式中

(2)

(3)

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

(4)

(5)

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

(6)

(7)

(8)

(9)

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

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

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

3.1 Zha-Bilgen分裂

Zha等[21]將歐拉方程的通量分裂成對(duì)流通量和壓力通量?jī)刹糠帧?/p>

(10)

對(duì)流通量Fc的雅可比矩陣為

(11)

Ac的特征值為

(12)

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

(13)

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

(14)

Ap的特征值為

(15)

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

(16)

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

3.2 FDS -ZB -BVD格式

3.2.1 壓力數(shù)值通量

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

(17)

(18)

(19)

3.2.2 對(duì)流數(shù)值通量

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

(20)

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

(21)

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

(22)

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

(23)

為了減少數(shù)值耗散，利用THINC函數(shù)來(lái)重構(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)

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

3.3 穩(wěn)定性分析

低耗散的數(shù)值格式在計(jì)算強(qiáng)激波問(wèn)題時(shí)會(huì)遭遇激波不穩(wěn)定現(xiàn)象。采用矩陣穩(wěn)定性分析方法[26]對(duì)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)通過(guò)Rankine -Hugoniot關(guān)系式來(lái)得到

(29)

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

(30)

(31)

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

(32)

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

(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)激波問(wèn)題時(shí)，F(xiàn)DS -ZB -BVD格式會(huì)比Roe格式表現(xiàn)出更好的魯棒性。第4節(jié)的數(shù)值實(shí)驗(yàn)也證明了這一點(diǎn)。

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

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

4.1 一維算例

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

4.1.1 孤立的接觸間斷

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

(36)

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

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

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

4.1.2 Lax問(wèn)題

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

(37)

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

圖3 Lax問(wèn)題的密度分布

4.1.3 慢行激波問(wèn)題

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

(38)

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

圖4 慢行激波問(wèn)題的密度分布

4.1.4 一維爆轟波問(wèn)題

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

(39)

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

4.1.5 激波-熵波相互作用問(wèn)題

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

(40)

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

圖6 激波-熵波相互作用問(wèn)題的密度分布

4.2 二維算例

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

4.2.1 二維黎曼問(wèn)題

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

(41)

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

圖7 二維黎曼問(wèn)題的密度等值線

4.2.2 移動(dòng)接觸面問(wèn)題

將計(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)動(dòng)。圖8展示了t=0.15時(shí)的密度等值線，此時(shí)圓心的位置為(0.65,0.65)。可以看出，F(xiàn)DS -ZB -BVD格式大大改善了接觸面的分辨率。

圖8 移動(dòng)接觸面問(wèn)題的密度等值線

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

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

(43)

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

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

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

(44)

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

4.2.5 二維Sedov爆轟波問(wèn)題

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

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

圖11 二維Sedov爆轟波問(wèn)題的壓力等值線

5 結(jié) 論

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