胡立軍 袁海專(zhuān) 杜玉龍

1) (衡陽(yáng)師范學(xué)院數(shù)學(xué)與統(tǒng)計(jì)學(xué)院,衡陽(yáng) 421002)

2) (湘潭大學(xué)數(shù)學(xué)與計(jì)算科學(xué)學(xué)院,湘潭 411105)

3) (北京航空航天大學(xué)數(shù)學(xué)科學(xué)學(xué)院,北京 100191)

使用低耗散激波捕捉格式對(duì)高超聲速流動(dòng)問(wèn)題進(jìn)行數(shù)值模擬時(shí)經(jīng)常會(huì)遭受不同形式的激波不穩(wěn)定性.本文基于二維無(wú)黏可壓縮Euler方程,對(duì)低耗散HLLEM格式進(jìn)行激波穩(wěn)定性分析.結(jié)果表明: 激波面橫向通量中切向速度的擾動(dòng)增長(zhǎng)誘發(fā)了格式的不穩(wěn)定性.通過(guò)增加耗散來(lái)治愈HLLEM格式的激波不穩(wěn)定性.為了避免引入過(guò)多的耗散進(jìn)而影響剪切層的分辨率,定義激波探測(cè)函數(shù)和亞聲速區(qū)探測(cè)函數(shù),使得只有在計(jì)算激波層亞聲速區(qū)的橫向數(shù)值通量時(shí)才增加耗散,其余地方的數(shù)值通量依然采用低耗散的HLLEM格式來(lái)計(jì)算.穩(wěn)定性分析和數(shù)值模擬的結(jié)果表明,改進(jìn)的HLLEM格式不僅保留了原格式高分辨率的優(yōu)點(diǎn),還大大提高了格式的魯棒性,在計(jì)算強(qiáng)激波問(wèn)題時(shí)能夠有效地抑制不穩(wěn)定現(xiàn)象的發(fā)生.

1 引 言

盡管流體力學(xué)數(shù)值方法的研究在最近幾十年取得了巨大的進(jìn)步,然而依然存在一些挑戰(zhàn)性的問(wèn)題,涉及激波、剪切層、激波/湍流邊界層干擾等現(xiàn)象的高超聲速流動(dòng)問(wèn)題的數(shù)值模擬就是其中之一[1,2].構(gòu)造精確、高效、魯棒性好的激波捕捉格式是精確模擬高速流動(dòng)問(wèn)題的關(guān)鍵.近似黎曼解法器是一類(lèi)流行的激波捕捉方法.一個(gè)好的黎曼解法器不僅要能夠精確地分辨各種間斷,還要具有很好的魯棒性,能夠抑制激波不穩(wěn)定現(xiàn)象的發(fā)生.然而,能夠精確分辨接觸間斷和剪切波的低耗散黎曼解法器(例如Roe,HLLC和HLLEM)在計(jì)算強(qiáng)激波問(wèn)題時(shí)都會(huì)產(chǎn)生嚴(yán)重的不穩(wěn)定現(xiàn)象,從而影響數(shù)值模擬的可靠性.因此,分析低耗散數(shù)值格式激波不穩(wěn)定性產(chǎn)生的根源并且在保留高分辨率優(yōu)點(diǎn)的同時(shí)去治愈格式的激波不穩(wěn)定性成為了計(jì)算流體力學(xué)數(shù)值方法研究的一個(gè)重要問(wèn)題.

利用線(xiàn)性擾動(dòng)分析方法,Quirk[3]指出壓力和密度場(chǎng)的相互作用誘發(fā)了某些格式的不穩(wěn)定性,并且在強(qiáng)激波附近使用耗散的HLLE格式來(lái)消除Roe格式在數(shù)值模擬中的不穩(wěn)定現(xiàn)象.Gressier和Moschetta[4]利用特殊的擾動(dòng)分析得出結(jié)論: 能夠精確捕捉接觸間斷的數(shù)值格式通常都會(huì)遭受激波不穩(wěn)定性的困擾.Dumbser等[5]設(shè)計(jì)了一種矩陣穩(wěn)定性分析方法來(lái)分析格式激波不穩(wěn)定性產(chǎn)生的根源并且發(fā)現(xiàn)在計(jì)算二維穩(wěn)態(tài)激波時(shí)數(shù)值格式的不穩(wěn)定性與數(shù)值激波的結(jié)構(gòu)密切相關(guān).謝文佳等[6]使用奇偶失聯(lián)線(xiàn)性方法對(duì)某些迎風(fēng)格式進(jìn)行穩(wěn)定性分析,結(jié)果表明: 具有充足剪切黏性的數(shù)值格式可以避免一類(lèi)激波不穩(wěn)定現(xiàn)象(“紅玉”現(xiàn)象)的發(fā)生,通過(guò)分析,Shen 等[7]也得到了同樣的結(jié)論.通過(guò)增加數(shù)值耗散來(lái)治愈格式的不穩(wěn)定性已經(jīng)成為了一種主流方法.一般來(lái)說(shuō),有以下三種方式來(lái)增加數(shù)值耗散: 將耗散格式和低耗散格式進(jìn)行混合,增加數(shù)值格式的人工黏性,控制數(shù)值格式的耗散機(jī)制.Kim等[8]構(gòu)造了一種HLLC-HLLE混合格式,通過(guò)定義開(kāi)關(guān)函數(shù)將強(qiáng)激波附近數(shù)值通量的計(jì)算由不穩(wěn)定的HLLC格式切換成穩(wěn)定的HLLE 格 式.Wu 等[9]和 Shen 等[10]構(gòu) 造 了一 種Roe-HLLE混合格式來(lái)消除Roe格式的不穩(wěn)定性.與全面混合不同的是,他們僅僅在質(zhì)量方程和切向動(dòng)量方程將這兩種格式進(jìn)行混合.Ren[11]提出了一種旋轉(zhuǎn)Roe黎曼解法器,通過(guò)自動(dòng)地引入人工黏性來(lái)抑制原始Roe格式不穩(wěn)定現(xiàn)象的發(fā)生.Rodionov[12]用人工黏性系數(shù)代替分子黏性系數(shù)來(lái)治愈Godunov-型數(shù)值格式的“紅玉”現(xiàn)象并且在文獻(xiàn)[13,14]中分別討論了高階Godunov-型格式和三維情形下的激波不穩(wěn)定性.Chen等[15]發(fā)現(xiàn)Roe格式和HLLE格式在亞聲速區(qū)剪切速度的差異可以解釋這兩種數(shù)值格式不同的激波行為,并且在激波附近亞聲速區(qū)增加Roe格式的剪切黏性來(lái)抑制不穩(wěn)定現(xiàn)象的發(fā)生.Simon和Mandal[16,17]使用反擴(kuò)散控制和選波修正兩種不同的策略來(lái)控制HLLC格式的耗散機(jī)制進(jìn)而抑制格式的不穩(wěn)定性,并且使用同樣的策略來(lái)構(gòu)造激波穩(wěn)定的HLLEM格式.Xie等[18]利用相鄰網(wǎng)格的壓力比定義開(kāi)關(guān)函數(shù)來(lái)控制HLLEM格式的耗散水平,進(jìn)而治愈其不穩(wěn)定性.與增加數(shù)值耗散來(lái)抑制激波不穩(wěn)定性的流行方法不同的是,Fleischmann等[19]通過(guò)修改Roe格式特征值的計(jì)算進(jìn)而得到一種耗散更小、魯棒性更好的Roe-M格式.

盡管眾多研究人員對(duì)于數(shù)值格式的激波不穩(wěn)定性進(jìn)行了大量的研究,但是關(guān)于激波不穩(wěn)定性產(chǎn)生的根源仍然存在爭(zhēng)議,因此很多治愈方法都存在一定的局限性.本文以低耗散HLLEM格式為例,采用線(xiàn)性穩(wěn)定性分析方法來(lái)推導(dǎo)兩種特殊擾動(dòng)情形下的擾動(dòng)發(fā)展方程,并且結(jié)合相應(yīng)的數(shù)值實(shí)驗(yàn)來(lái)探究激波不穩(wěn)定性產(chǎn)生的根源.根據(jù)穩(wěn)定性分析的結(jié)論提出一種改進(jìn)的HLLEM格式來(lái)治愈原格式的不穩(wěn)定性并且進(jìn)行數(shù)值實(shí)驗(yàn)來(lái)驗(yàn)證改進(jìn)格式的魯棒性和精度.

2 控制方程組和HLLEM黎曼解法器

2.1 二維歐拉方程組

考慮二維理想氣體流動(dòng)的守恒律方程組

式中,守恒量U,矢通量F(U) 和G(U) 的定義分別為:

其中ρ為密度,u和v分別為全局坐標(biāo)系中x方向和y方向的流體速度,p為靜壓力,E為總能.使用理想氣體狀態(tài)方程來(lái)使守恒律系統(tǒng)封閉:

其中比熱比γ=1.4.求解方程組(1)的守恒型數(shù)值方法可以表示成

其中 ?t為時(shí)間步長(zhǎng),?x和 ?y為空間步長(zhǎng),Fi+1/2,j和Gi,j+1/2為x方向和y方向的數(shù)值通量.

2.2 HLLEM黎曼解法器

Einfeldt等[20]在HLLE通量的基礎(chǔ)上添加與線(xiàn)性退化場(chǎng)相對(duì)應(yīng)的反擴(kuò)散項(xiàng),從而得到了一種可以精確分辨接觸間斷和剪切波的HLLEM格式,其通量表達(dá)式為

其中,SL和SR表示左、右波速,R2,3表示通量雅克比矩陣與線(xiàn)性退化場(chǎng)相對(duì)應(yīng)的右特征向量,α2,3為對(duì)應(yīng)的波強(qiáng),δ2,3表示控制耗散大小的系數(shù),他們的計(jì)算公式分別為:

與其他可以捕捉接觸間斷和剪切波的低耗散格式一樣,使用HLLEM格式計(jì)算多維強(qiáng)激波問(wèn)題時(shí)也會(huì)產(chǎn)生不同形式的激波不穩(wěn)定現(xiàn)象,這大大限制了HLLEM格式在高超聲速流動(dòng)問(wèn)題數(shù)值模擬中的應(yīng)用.接下來(lái)首先探究HLLEM格式激波不穩(wěn)定性產(chǎn)生的根源,進(jìn)而采用針對(duì)性的方法來(lái)治愈其不穩(wěn)定性.

3 HLLEM 格式激波不穩(wěn)定性的機(jī)理分析

3.1 線(xiàn)性穩(wěn)定性分析

流體介質(zhì)以速度u0=M0a0沿著x正方向(縱向) 運(yùn)動(dòng),沿y方向 (橫向) 的速度v0=0 ,其中M0為馬赫數(shù),為當(dāng)?shù)芈曀?為了弄清楚是哪個(gè)方向的擾動(dòng)增長(zhǎng)造成了HLLEM格式的不穩(wěn)定性,考慮在初始分布上增加兩類(lèi)特殊的擾動(dòng).

首先,在x方向上增加奇偶對(duì)稱(chēng)擾動(dòng),擾動(dòng)后的流場(chǎng)狀態(tài)為:

在超聲速(M0>1 )的情況下,擾動(dòng)發(fā)展方程為:

其中σx=?t/?x.(11)式的特征值為:

當(dāng)σx<1/(u0+a0) 時(shí) ,此時(shí)所有擾動(dòng)量都能有效衰減.選取兩組不同的初始擾動(dòng)量和初值:

經(jīng)過(guò)時(shí)間迭代20步,每個(gè)擾動(dòng)量的發(fā)展趨勢(shì)如圖1所示,其中σx取 0.2.

從圖1 可以看到,當(dāng)x方向存在擾動(dòng)時(shí),在超聲速情形下,所有擾動(dòng)量都會(huì)衰減到0,此時(shí)數(shù)值格式是穩(wěn)定的.

圖1 超聲速情形下擾動(dòng)量的發(fā)展曲線(xiàn) (a)初始擾動(dòng)量和初值(I);(b)初始擾動(dòng)量和初值(II)Fig.1.Evolutionary curves of perturbation under supersonic condition: (a) Initial perturbation values and initial conditions (I);(b) initial perturbation values and initial conditions (II).

在亞聲速( 0

(14)式的特征值為:

當(dāng)σx<1/(a0+u0) 時(shí) ,1(i=1,2,3,4) ,此時(shí)所有擾動(dòng)量均能有效衰減.同樣選取兩組不同的初始擾動(dòng)量和初值:

經(jīng)過(guò)時(shí)間迭代20步,每個(gè)擾動(dòng)量的發(fā)展趨勢(shì)如圖2所示,其中σx取 0.4.

從圖2 可以看到,當(dāng)x方向存在擾動(dòng)時(shí),在亞聲速情形下,所有擾動(dòng)量也會(huì)衰減到0.因此當(dāng)流體沿著x方向運(yùn)動(dòng)時(shí),該方向的擾動(dòng)不會(huì)誘發(fā)HLLEM格式的不穩(wěn)定性.

接下來(lái),考慮在y方向增加奇偶對(duì)稱(chēng)擾動(dòng),x方向上流場(chǎng)無(wú)擾動(dòng).此時(shí)的流場(chǎng)狀態(tài)為

由于在x方向沒(méi)有擾動(dòng),因此數(shù)值通量Fi?1/2,j=Fi+1/2,j=0.使用 HLLEM格式計(jì)算數(shù)值通量Gi,j?1/2和Gi,j+1/2,進(jìn)而得到擾動(dòng)量的發(fā)展方程為:

(18)式的特征值為:

圖2 亞聲速情形下擾動(dòng)量的發(fā)展曲線(xiàn) (a)初始擾動(dòng)量和初值(III);(b)初始擾動(dòng)量和初值(IV)Fig.2.Evolutionary curves of perturbation under subsonic condition: (a) Initial perturbation values and initial conditions (III);(b) initial perturbation values and initial conditions (IV).

其中σy=?t/?y.當(dāng)σy<1/a0時(shí),和是衰減的,而密度擾動(dòng)和切向速度擾動(dòng)不會(huì)衰減.設(shè)初值ρ0=1.4,u0=1 ,p0=1 ,選取以下四組不同的初始擾動(dòng)量來(lái)觀(guān)察在該擾動(dòng)下各個(gè)物理量的擾動(dòng)發(fā)展趨勢(shì)

3.2 激波不穩(wěn)定性的多維特性

關(guān)于激波不穩(wěn)定性是一種一維激波的異?，F(xiàn)象還是一種多維現(xiàn)象仍然是有爭(zhēng)議的.Chauvat等[21]和Kitamura等[22]發(fā)現(xiàn)激波不穩(wěn)定性與內(nèi)部激波結(jié)構(gòu)有關(guān)并且主張它是一種一維激波的異常現(xiàn)象.然而,純粹的一維算例不會(huì)遭受激波不穩(wěn)定現(xiàn)象的困擾,并且3.1節(jié)的線(xiàn)性穩(wěn)定性分析的結(jié)論表明:y方向的擾動(dòng)會(huì)造成x方向的激波不穩(wěn)定性,反之亦然.因此,本文認(rèn)為它是一種多維現(xiàn)象.

圖3 橫向擾動(dòng)下擾動(dòng)量的發(fā)展曲線(xiàn) (a)初始擾動(dòng)量 (V);(b)初始擾動(dòng)量 (VI);(c)初始擾動(dòng)量 (VII);(d) 初始擾動(dòng)量 (VIII)Fig.3.Evolutionary curves of perturbation in transverse direction: (a) Initial perturbation values (V);(b) initial perturbation values (VI);(c) initial perturbation values (VII);(d) initial perturbation values (VIII).

圖4 二維 Sedov 爆轟波問(wèn)題的密度等值線(xiàn) (a) HLLEM;(b) x-HLLE+y-HLLEM (c) x-HLLEM+y-HLLEFig.4.Density contours for 2D Sedov blast wave problem: (a) HLLEM;(b) x-HLLE+y-HLLEM (c) x-HLLEM+y-HLLE.

接下來(lái),計(jì)算文獻(xiàn)[19]中的二維Sedov爆轟波問(wèn)題來(lái)證明激波不穩(wěn)定性的多維特性.初始時(shí),區(qū)域 [ 0,2.4]×[0,2.4]中心的壓力pin=3.5×105,其余地方的壓力pout=10?10,整個(gè)區(qū)域的密度為1,速度u0和v0均為 0.計(jì)算中使用 4 80×480 的笛卡爾網(wǎng)格,上、下、左、右都使用反射邊界條件.圖4為計(jì)算得到的二維Sedov爆轟波問(wèn)題的密度等值線(xiàn),可以看到,采用 HLLEM 格式進(jìn)行計(jì)算時(shí),在激波面和坐標(biāo)軸平行的位置出現(xiàn)了明顯的“紅玉”現(xiàn)象.當(dāng)x方向的通量采用耗散的HLLE格式計(jì)算時(shí),y方向的“紅玉”現(xiàn)象消失了,反之亦然.該數(shù)值實(shí)驗(yàn)的結(jié)果與線(xiàn)性穩(wěn)定性分析的結(jié)論相一致,從而進(jìn)一步證明了激波不穩(wěn)定現(xiàn)象的多維特性.

4 一種改進(jìn)的 HLLEM 格式

激波面橫向擾動(dòng)的發(fā)展方程(18)式表明: 密度和切向速度的擾動(dòng)不會(huì)衰減.由于激波不穩(wěn)定性是一種多維現(xiàn)象,并且一維問(wèn)題存在熵波但是沒(méi)有剪切波.因此我們認(rèn)為: 激波面橫向通量中切向速度的擾動(dòng)增長(zhǎng)造成了HLLEM格式的激波不穩(wěn)定性.通過(guò)增加橫向通量的耗散來(lái)抑制不穩(wěn)定現(xiàn)象的發(fā)生,其具體表達(dá)式為

在耗散項(xiàng)的作用下,橫向擾動(dòng)的發(fā)展方程式為:

(22)式表明橫向通量中的切向速度的擾動(dòng)也會(huì)衰減.

4.1 激波探測(cè)函數(shù)

在實(shí)際的數(shù)值模擬中,只需要在激波附近增加耗散就可以抑制不穩(wěn)定現(xiàn)象的發(fā)生,為了避免在不適當(dāng)?shù)奈恢迷黾雍纳⑦M(jìn)而影響剪切層的分辨率,利用網(wǎng)格界面的壓力比來(lái)探測(cè)激波的位置:

顯然,h∈(0,1) ,且在激波附近壓力差較大,從而h→0.由前面的分析可知: 在x方向增加耗散可以抑制y方向的不穩(wěn)定性,反之亦然.因此,計(jì)算x方向的數(shù)值通量Fi+1/2,j時(shí),只需檢測(cè)與之相鄰的y方向的網(wǎng)格界面.計(jì)算y方向的數(shù)值通量Gi,j+1/2時(shí),只需檢測(cè)與之相鄰的x方向的網(wǎng)格界面.即

由于不穩(wěn)定性的嚴(yán)重程度會(huì)隨著激波強(qiáng)度的增加而增加,我們采用文獻(xiàn)[15]中的余弦函數(shù)來(lái)定義激波探測(cè)函數(shù):

此外還給出另外兩種激波探測(cè)函數(shù):

圖5畫(huà)出了這三種激波探測(cè)函數(shù)的曲線(xiàn)圖.從中可以看到三個(gè)函數(shù)的取值都會(huì)隨著激波強(qiáng)度的減弱而減小,且對(duì)于固定的h,滿(mǎn)足g1>g>g2.

圖5 三種激波探測(cè)函數(shù)的函數(shù)曲線(xiàn)Fig.5.Curves for three different shock-detecting functions.

4.2 M-HLLEM格式

Xu和Li[23]在分析定常激波問(wèn)題時(shí)發(fā)現(xiàn)激波不穩(wěn)定性是由激波層亞聲速區(qū)的擾動(dòng)增長(zhǎng)導(dǎo)致的,這一結(jié)論得到了眾多研究人員的支持[6,15,19].因此定義亞聲速區(qū)探測(cè)函數(shù)

式中M表示馬赫數(shù).在亞聲速區(qū) 0<φ<1 ,在超聲速區(qū)φ=0.

因此,在一般的結(jié)構(gòu)化四邊形網(wǎng)格下,最終的耗散項(xiàng)可以表示為

其中n=(nx,ny) 為網(wǎng)格界面單位外法向量,g和φ分別為激波探測(cè)函數(shù)和亞聲速區(qū)探測(cè)函數(shù).因此,一種改進(jìn)的HLLEM格式(M-HLLEM,Modified HLLEM)數(shù)值通量的表達(dá)式為

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

本節(jié)計(jì)算一些典型的算例來(lái)檢驗(yàn)本文構(gòu)造的M-HLLEM格式的魯棒性和精度.

5.1 隨機(jī)數(shù)值噪聲問(wèn)題

馬赫數(shù)為6的平面激波沿著x方向運(yùn)動(dòng),初始時(shí)激波位于x=5 ,右側(cè)波前的初始條件為(ρ0,u0,v0,p0)=(1.4,0,0,1),左側(cè)的波后狀態(tài)可以通過(guò)激波和馬赫數(shù)之間的關(guān)系式來(lái)得到.計(jì)算區(qū)域 [ 0,1000]×[0,20]被劃分成 1 000×20 的笛卡爾網(wǎng)格.這是一個(gè)純粹的一維激波,為了誘發(fā)不穩(wěn)定性,在波前的初始分布上增加取值從?0.5×10?5到0.5×10?5的隨機(jī)數(shù)值擾動(dòng)

其中αk(k=1,2,3,4) 為?0.5 到 0.5 之間的隨機(jī)數(shù).

在該算例中,y方向速度的最大量值max(|v|)可以用來(lái)衡量不穩(wěn)定性的嚴(yán)重程度.圖6展示了t=120時(shí),原始的HLLEM格式以及與不同的激波探測(cè)函數(shù)結(jié)合的改進(jìn)格式的密度等值線(xiàn)圖.圖7展示了y方向速度的最大量值隨著時(shí)間變化的曲線(xiàn)圖.可以看到,HLLEM格式的激波結(jié)構(gòu)完全被破壞,且y方向速度的最大量值也從 1 0?5量級(jí)增長(zhǎng)到100量級(jí).而三種改進(jìn)格式都得到了清晰的激波面,且y方向速度的最大量值也始終保持在初始時(shí)的10?5量級(jí)附近.三種激波探測(cè)函數(shù)的數(shù)值結(jié)果之間幾乎沒(méi)有差異,因此后面的數(shù)值算例均選用流行的余弦函數(shù)(25)式來(lái)計(jì)算其函數(shù)值.

圖6 隨機(jī)數(shù)值噪聲問(wèn)題的密度等值線(xiàn) (a) HLLEM;(b)M-HLLEM-g;(c) M-HLLEM-g1;(d) M-HLLEM-g2Fig.6.Density contours of random numerical noise problem: (a) HLLEM;(b) M-HLLEM-g;(c) M-HLLEM-g1;(d)M-HLLEM-g2.

圖7 隨機(jī)數(shù)值噪聲問(wèn)題 y 方向速度的最大量值Fig.7.Maximum magnitude of velocity in y-direction of random numerical noise problem.

5.2 高超聲速繞柱流問(wèn)題

高超聲速繞柱流問(wèn)題經(jīng)常用來(lái)檢驗(yàn)數(shù)值格式是否會(huì)遭受一種嚴(yán)重的不穩(wěn)定現(xiàn)象—“紅玉”現(xiàn)象(carbuncle).馬赫數(shù)為20的自由來(lái)流流經(jīng)半徑為1的圓柱體,整個(gè)區(qū)域流場(chǎng)的初始條件(ρ0,u0,v0,p0)=(1.4,20,0,1),具體的計(jì)算區(qū)域、網(wǎng)格劃分和邊界條件可以參考文獻(xiàn)[24].本文采用320×40的結(jié)構(gòu)化貼體四邊形網(wǎng)格來(lái)計(jì)算.圖8畫(huà)出了時(shí)間t=4 時(shí),使用HLLEM和M-HLLEM格式計(jì)算得到的密度等值線(xiàn).可以清晰地看到,HLLEM格式出現(xiàn)了嚴(yán)重的“紅玉”現(xiàn)象,而M-HLLEM格式消除了不穩(wěn)定現(xiàn)象,得到了具有清晰激波面的穩(wěn)態(tài)弓形激波.

圖8 高超聲速繞柱流問(wèn)題的密度等值線(xiàn) (a) HLLEM;(b) M-HLLEMFig.8.Density contours of hypersonic flow over a cylinder:(a) HLLEM;(b) M-HLLEM.

5.3 雙馬赫反射問(wèn)題

Woodward和Colella[25]提出的雙馬赫反射問(wèn)題是檢驗(yàn)數(shù)值格式魯棒性的一個(gè)著名算例.馬赫數(shù)為10的斜激波向底部壁面運(yùn)動(dòng),激波面與壁面所成 角 度 為 60°,計(jì) 算 區(qū) 域 [ 0,4]×[0,1]被 劃 分 成480×120的笛卡爾網(wǎng)格.詳細(xì)的初始條件、邊界條件可參考文獻(xiàn) [25].圖9畫(huà)出了時(shí)間t=0.2 時(shí),HLLEM和M-HLLEM格式的計(jì)算結(jié)果.可以看到,HLLEM格式出現(xiàn)了彎曲的馬赫桿和一個(gè)三角點(diǎn),這是一種明顯的非物理現(xiàn)象.而M-HLLEM格式消除了不穩(wěn)定現(xiàn)象,得到了清晰的馬赫桿.

圖9 雙馬赫反射問(wèn)題的密度等值線(xiàn) (a) HLLEM;(b) M-HLLEMFig.9.Density contours of double Mach reflection problem:(a) HLLEM;(b) M-HLLEM.

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

計(jì)算3.2節(jié)描述的二維Sedov爆轟波問(wèn)題,由于該算例涉及高壓力比的強(qiáng)激波問(wèn)題,因此可以用來(lái)作為檢驗(yàn)數(shù)值格式性能的一個(gè)挑戰(zhàn)性測(cè)試.圖10畫(huà)出了時(shí)間t=0.1 時(shí),HLLEM和 M-HLLEM格式的壓力等值線(xiàn).可以清晰地看到,HLLEM格式在激波面和網(wǎng)格線(xiàn)平行的位置出現(xiàn)了四個(gè)明顯的“紅玉”現(xiàn)象.而M-HLLEM格式有效地抑制了不穩(wěn)定現(xiàn)象的發(fā)生.

5.5 二維無(wú)黏接觸間斷問(wèn)題

數(shù)值格式可以用來(lái)計(jì)算黏性流的前提條件是能夠精確捕捉接觸間斷.該算例用來(lái)檢驗(yàn)改進(jìn)的HLLEM格式捕捉接觸間斷的能力.區(qū)域[0,1]×[0,1]被均勻劃分成 1 0×10 的正方形網(wǎng)格,上、下兩部分不同速度的兩種流體的初始條件為:

圖10 二維 Sedov 爆轟波問(wèn)題的壓力等值線(xiàn) (a) HLLEM;(b) M-HLLEMFig.10.Pressure contours of 2D Sedov blast wave problem:(a) HLLEM;(b) M-HLLEM.

圖11展示了計(jì)算迭代1000步后,x=0.5 處的密度分布.可以看到,HLLE格式得到了具有較大耗散的解,而HLLEM和M-HLLEM格式都能準(zhǔn)確捕捉該接觸間斷.這說(shuō)明M-HLLEM格式保留了原HLLEM格式精確分辨接觸間斷的優(yōu)點(diǎn).

圖11 二維無(wú)黏接觸間斷問(wèn)題的密度分布Fig.11.Density distribution of 2D inviscid contact discontinuity problem.

6 結(jié) 論

本文構(gòu)造了一種激波穩(wěn)定的HLLEM格式.線(xiàn)性穩(wěn)定性分析和相關(guān)的數(shù)值實(shí)驗(yàn)表明: 橫向通量中切向速度的擾動(dòng)增長(zhǎng)誘發(fā)了HLLEM格式的激波不穩(wěn)定性.通過(guò)增加耗散來(lái)治愈格式的不穩(wěn)定性,為了防止過(guò)多的耗散影響剪切層的分辨率,定義激波探測(cè)函數(shù)和亞聲速區(qū)探測(cè)函數(shù),使得耗散項(xiàng)僅僅添加在激波層亞聲速區(qū)域的橫向通量上,其余的數(shù)值通量依然采用低耗散的HLLEM格式來(lái)計(jì)算,從而在不犧牲精度的前提下,提高格式的魯棒性.數(shù)值模擬的結(jié)果表明: 改進(jìn)的HLLEM格式消除了激波不穩(wěn)定現(xiàn)象并且保留了原始HLLEM格式精確分辨接觸間斷的優(yōu)點(diǎn),因此該格式可以廣泛應(yīng)用于高超聲速可壓縮流的數(shù)值模擬中.采用類(lèi)似的方法來(lái)改進(jìn)其他的低耗散格式,結(jié)合高階精度重構(gòu)方法并且將其推廣到三維情形來(lái)計(jì)算不同的守恒律系統(tǒng)可以作為未來(lái)的研究工作.

感謝中國(guó)科學(xué)院計(jì)算數(shù)學(xué)研究所袁禮研究員的討論.