迎風(fēng)格式在接觸間斷的數(shù)值耗散及其誘導(dǎo)誤差

韓 芳,魏雁昕,劉 君

(大連理工大學(xué) 航空航天學(xué)院, 遼寧 大連 116024)

在計(jì)算流體力學(xué)(computational fluid dynamics, CFD)中,數(shù)值耗散主要來(lái)自空間離散格式離散對(duì)流項(xiàng)和時(shí)間離散格式離散時(shí)間項(xiàng)所產(chǎn)生的截?cái)嗾`差,其大小與格式的計(jì)算精度及穩(wěn)定性息息相關(guān)。因此,CFD研究者們常常采用降低數(shù)值耗散的方法來(lái)構(gòu)造更高精度的計(jì)算格式[1-3]。空間離散格式包含差分格式及對(duì)流通量導(dǎo)數(shù)的求解兩部分,本文主要討論對(duì)流通量導(dǎo)數(shù)計(jì)算方法的數(shù)值耗散特性,而迎風(fēng)格式作為目前CFD求解對(duì)流通量的主流格式[4],成為本文的研究對(duì)象。

迎風(fēng)格式主要分為三類(lèi):矢通量分裂(flux vector splitting,FVS)格式、通量差分分裂(flux difference splitting,FDS)格式以及混合格式。現(xiàn)今CFD對(duì)迎風(fēng)格式耗散性的主要觀(guān)點(diǎn)[5]為:以Steger-Warming格式[6]、Van Leer格式[7]等為代表的FVS格式將對(duì)流通量分裂為正負(fù)兩部分,在其分裂過(guò)程中使用的簡(jiǎn)化策略誘導(dǎo)出過(guò)大的數(shù)值耗散,降低了格式對(duì)線(xiàn)性波、非線(xiàn)性波等的模擬精度;以Roe格式[8]、HLLC格式[9]為代表的FDS格式的基礎(chǔ)是求解局部黎曼問(wèn)題,其數(shù)值耗散較低,這使得FDS格式能夠準(zhǔn)確地捕捉到激波、接觸間斷等非線(xiàn)性波、線(xiàn)性波,但其在強(qiáng)激波處的計(jì)算穩(wěn)定性較差,例如Roe格式的“紅玉”現(xiàn)象問(wèn)題;混合格式作為以上兩種格式的結(jié)合,既具備了FVS格式在非線(xiàn)性波捕捉上的魯棒性,也具備了FDS格式在線(xiàn)性波捕捉上的高分辨率,其耗散較低且穩(wěn)定性?xún)?yōu)于FDS格式,例如AUSM+格式[10]、LDFSS格式[11]等。

迎風(fēng)格式的數(shù)值耗散大小可表現(xiàn)為其對(duì)間斷的捕捉能力大小。一般認(rèn)為,在相同的計(jì)算條件下,捕捉到間斷過(guò)渡區(qū)越窄,則格式的數(shù)值耗散越低,捕捉間斷的能力越強(qiáng),例如,Sod[12]、Lax[13]、Shu-Osher[14]等問(wèn)題。對(duì)存在間斷相互干擾的問(wèn)題,則以捕捉流場(chǎng)小尺度結(jié)構(gòu)的能力來(lái)衡量格式的數(shù)值耗散大小,在相同的計(jì)算條件下,捕捉到的小尺度結(jié)構(gòu)越多、越清晰,則表明格式的數(shù)值耗散越小,例如2-D Riemann問(wèn)題[15]等。不論是對(duì)間斷的捕捉,還是對(duì)流場(chǎng)小尺度結(jié)構(gòu)的捕捉,以上問(wèn)題對(duì)迎風(fēng)格式數(shù)值耗散的研究都集中于間斷本身或間斷附近,對(duì)在數(shù)值耗散作用下產(chǎn)生的激波間斷或接觸間斷等對(duì)下游流場(chǎng)影響的研究很少。文獻(xiàn)[16-18]在計(jì)算運(yùn)動(dòng)激波問(wèn)題時(shí)發(fā)現(xiàn),捕捉法計(jì)算運(yùn)動(dòng)激波會(huì)產(chǎn)生非物理波動(dòng),這一波動(dòng)會(huì)向流場(chǎng)下游進(jìn)行傳播,進(jìn)而污染激波下游區(qū)域,改變激波下游流場(chǎng)的結(jié)構(gòu)分布,這一現(xiàn)象在以上三種迎風(fēng)格式中都存在,但文章并未對(duì)非物理波動(dòng)的產(chǎn)生機(jī)理進(jìn)行分析。文獻(xiàn)[19-20]采用以上三種迎風(fēng)格式對(duì)馬赫數(shù)為3的運(yùn)動(dòng)正激波進(jìn)行了模擬,計(jì)算結(jié)果表明無(wú)論采用哪種格式,在激波由初始間斷形成數(shù)值過(guò)渡區(qū)的過(guò)程中,在激波下游都會(huì)產(chǎn)生一個(gè)等熵波和一個(gè)非等熵波。這兩個(gè)波在不同格式下大小不同,但始終存在,影響著激波下游流場(chǎng)結(jié)構(gòu)分布。文獻(xiàn)[21-22]在研究使用加權(quán)本質(zhì)無(wú)振蕩(weighted essentially non-oscillatory,WENO)格式計(jì)算含接觸間斷的可壓縮流時(shí)發(fā)現(xiàn),FVS格式特征值逐點(diǎn)分裂的不兼容性以及質(zhì)量方程、動(dòng)量方程、能量方程中各分量非線(xiàn)性離散的不一致性會(huì)導(dǎo)致流場(chǎng)中出現(xiàn)非物理振蕩,推薦使用全局Lax-Friedrichs(L-F)分裂格式,并在WENO格式實(shí)現(xiàn)過(guò)程中使用一組共同的權(quán)值來(lái)消除非物理振蕩。

為了進(jìn)一步分析不同迎風(fēng)格式數(shù)值耗散的產(chǎn)生機(jī)理,同時(shí)認(rèn)清上文所提到的非物理波動(dòng)是否與格式的數(shù)值耗散有關(guān),采用具有單一流場(chǎng)結(jié)構(gòu)的接觸間斷作為數(shù)值模擬對(duì)象開(kāi)展研究。綜合不同迎風(fēng)格式的數(shù)值耗散及其非物理波動(dòng)在接觸間斷問(wèn)題中的表現(xiàn),可以為不同的流場(chǎng)問(wèn)題選擇更加合適的計(jì)算格式。

1 控制方程與數(shù)值方法

二維笛卡爾坐標(biāo)系(x,y,t)下無(wú)量綱化的守恒型Euler方程為:

(1)

式中,U為守恒變量,F、G分別為x、y方向的對(duì)流通量,其表達(dá)式為:

(2)

其中,ρ為密度,u、v分別為x、y方向的速度,p為壓力,e為氣體總比內(nèi)能,其計(jì)算式為:

(3)

式中,γ為比熱比,本文使用量熱完全氣體模型,γ=1.4。

采用有限差分法進(jìn)行數(shù)值求解,其中時(shí)間離散采用一階顯式格式,因此方程(1)的離散形式可寫(xiě)成:

(4)

以上公式及變量皆為無(wú)量綱化的,下文數(shù)值實(shí)驗(yàn)中的所有變量和數(shù)據(jù)也是無(wú)量綱化的。

文中以迎風(fēng)格式作為研究對(duì)象,分別選用Van Leer、Roe、AUSM+格式作為三種迎風(fēng)格式的代表進(jìn)行數(shù)值計(jì)算。

2 接觸間斷問(wèn)題

2.1 一維接觸間斷問(wèn)題

使用二維均勻正交網(wǎng)格計(jì)算一維流場(chǎng),計(jì)算區(qū)域設(shè)定為[0,10]×[0,1],網(wǎng)格量為200×20。無(wú)量綱化的初始條件為:

(5)

其中,xd為間斷的初始位置。將xd左側(cè)的區(qū)域記為1區(qū),將xd右側(cè)的區(qū)域記為2區(qū),因此ρ1為間斷左側(cè)的初始密度參數(shù)。

對(duì)一維接觸間斷問(wèn)題,由之后的計(jì)算結(jié)果可以看出,在初始間斷形成數(shù)值過(guò)渡層的過(guò)程中出現(xiàn)了非物理波動(dòng),其傳播速度為特征速度,為避免擾動(dòng)傳出邊界,設(shè)置間斷的初始位置xd滿(mǎn)足式(6):

(6)

其中,c1為1區(qū)的聲速,tstop為計(jì)算終止時(shí)刻。根據(jù)初始間斷左右流場(chǎng)馬赫數(shù)的不同,即流場(chǎng)初始條件的不同,設(shè)置以下五個(gè)算例,其具體參數(shù)如表1所示。

表1 不同馬赫數(shù)條件下的接觸間斷初始條件

為計(jì)算方便,本節(jié)統(tǒng)一取計(jì)算終止時(shí)刻tstop=2.0,根據(jù)給定的CFL(Courant-Friedrichs-Lewy)數(shù)確定時(shí)間步長(zhǎng),此處CFL=0.5。此外,分別選取CFL∈{0.01、0.05、0.1、0.5}對(duì)初始條件5的流場(chǎng)進(jìn)行計(jì)算,可知在允許范圍內(nèi),時(shí)間步長(zhǎng)的大小對(duì)計(jì)算結(jié)果影響很小,可忽略不計(jì)。

為了消除差分格式及限制器等因素的影響,選擇一階迎風(fēng)格式進(jìn)行數(shù)值計(jì)算,圖1～5分別給出了五種初始條件下的流場(chǎng)參數(shù)分布曲線(xiàn),包括密度ρ、速度u及壓力p分布曲線(xiàn)。

從圖1可以看出,在靜止流場(chǎng)中,Roe格式及AUSM+格式可以完整地保持接觸間斷不變,其密度、速度及壓力參數(shù)分布曲線(xiàn)皆保持初始狀態(tài)。而Van Leer格式計(jì)算的流場(chǎng)參數(shù)分布曲線(xiàn)不僅有密度耗散,速度及壓力曲線(xiàn)也有波動(dòng)。

(a) 密度參數(shù)分布曲線(xiàn)(a) Density parameter distribution curves

(b) 速度參數(shù)分布曲線(xiàn)(b) Velocity parameter distribution curves

(c) 壓力參數(shù)分布曲線(xiàn)(c) Pressure parameter distribution curves

從圖2中可以看出,在流場(chǎng)為全場(chǎng)超聲速條件時(shí),三種格式對(duì)接觸間斷的數(shù)值模擬結(jié)果相同,都是僅有密度的耗散,速度及壓力參數(shù)無(wú)變化。

(a) 密度參數(shù)分布曲線(xiàn)(a) Density parameter distribution curves

(b) 速度參數(shù)分布曲線(xiàn)(b) Velocity parameter distribution curves

(c) 壓力參數(shù)分布曲線(xiàn)(c) Pressure parameter distribution curves

從圖3～5可以看出,在流場(chǎng)中存在亞聲速區(qū)域時(shí),Roe格式及AUSM+格式捕捉的接觸間斷僅有密度的變化,速度及壓力參數(shù)無(wú)變化。Van Leer格式捕捉的接觸間斷存在速度及壓力參數(shù)的變化,其波動(dòng)曲線(xiàn)與在靜止流場(chǎng)中的結(jié)果相似,不僅在間斷處有波動(dòng),在間斷兩側(cè)也產(chǎn)生了兩個(gè)非物理波動(dòng)。但是,盡管Van Leer格式計(jì)算的接觸間斷有速度及壓力的變化,其密度變化曲線(xiàn)與Roe格式及AUSM+格式的計(jì)算結(jié)果曲線(xiàn)幾乎重合。

根據(jù)以上算例的計(jì)算結(jié)果,對(duì)一維接觸間斷,Roe格式及AUSM+格式在靜止流場(chǎng)中沒(méi)有耗散,在非靜止流場(chǎng)中有密度耗散;Van Leer格式則在所有流場(chǎng)中都有密度的耗散。圖6給出了不同時(shí)刻靜止流場(chǎng)(初始條件1)中非物理波動(dòng)的位置變化曲線(xiàn),從圖中可以看出,在上文中出現(xiàn)的非物理波動(dòng)分別以特征速度u-c1和u+c2逐步遠(yuǎn)離間斷。將間斷左側(cè)波記為“波1”,間斷右側(cè)波記為“波2”,則圖中v1、v2分別代表兩個(gè)波動(dòng)的傳播速度,c1、c2分別為當(dāng)?shù)攸c(diǎn)聲速。

圖7給出了使用初始條件1計(jì)算的接觸間斷流場(chǎng)壓力波動(dòng)幅值隨時(shí)間變化曲線(xiàn),從圖中可以看出,隨著時(shí)間步的推進(jìn),誘導(dǎo)誤差在遠(yuǎn)離接觸間斷的同時(shí),其幅值大小也會(huì)逐漸降低,但不會(huì)消失。結(jié)合圖6和圖7可以看出,對(duì)定常問(wèn)題,在計(jì)算時(shí)間足夠長(zhǎng)的情況下,誤差會(huì)運(yùn)動(dòng)出有限的計(jì)算區(qū)域,流場(chǎng)結(jié)構(gòu)會(huì)恢復(fù)原本的狀態(tài)。但對(duì)非定常問(wèn)題,波動(dòng)對(duì)流場(chǎng)的影響無(wú)法消除,在參考文獻(xiàn)[16-18]的算例中可以看到波動(dòng)對(duì)間斷下游流場(chǎng)參數(shù)的污染。

此外,這逐步遠(yuǎn)離間斷的非物理波動(dòng)不應(yīng)歸于Van Leer格式的數(shù)值耗散,而是一種數(shù)值誤差,因?yàn)閿?shù)值耗散的作用是“抹平”間斷,而這兩個(gè)非物理波動(dòng)對(duì)密度曲線(xiàn)的過(guò)渡區(qū)并無(wú)影響。

(a) 密度參數(shù)分布曲線(xiàn)(a) Density parameter distribution curves

(b) 速度參數(shù)分布曲線(xiàn)(b) Velocity parameter distribution curves

(c) 壓力參數(shù)分布曲線(xiàn)(c) Pressure parameter distribution curves

(a) 密度參數(shù)分布曲線(xiàn)(a) Density parameter distribution curves

(b) 速度參數(shù)分布曲線(xiàn)(b) Velocity parameter distribution curves

(c) 壓力參數(shù)分布曲線(xiàn)(c) Pressure parameter distribution curves

(a) 密度參數(shù)分布曲線(xiàn)(a) Density parameter distribution curves

(b) 速度參數(shù)分布曲線(xiàn)(b) Velocity parameter distribution curves

(c) 壓力參數(shù)分布曲線(xiàn)(c) Pressure parameter distribution curves

圖6 波動(dòng)位置隨時(shí)間變化曲線(xiàn)(初始條件1)Fig.6 Curve of fluctuation position with time(initial condition 1)

綜上所述,本文認(rèn)為迎風(fēng)格式的數(shù)值耗散受流場(chǎng)參數(shù)影響較大,在不同馬赫數(shù)條件下數(shù)值耗散大小不同,不過(guò)文獻(xiàn)[23]的計(jì)算結(jié)果也顯示,在低馬赫數(shù)(0.1, 0.01, 0.001)條件下,Roe格式的耗散大小與馬赫數(shù)無(wú)關(guān),而 HLL格式的數(shù)值耗散隨著馬赫數(shù)的減小而變大。此外,Van Leer格式捕捉接觸間斷時(shí)產(chǎn)生的非物理波動(dòng)是一種數(shù)值誤差,它們的存在影響了流場(chǎng)的速度及壓力參數(shù)分布,但對(duì)Van Leer格式在接觸間斷問(wèn)題中的數(shù)值耗散大小并無(wú)影響。

2.2 對(duì)數(shù)值計(jì)算結(jié)果的理論推導(dǎo)及機(jī)理分析

滿(mǎn)足Euler方程的初始間斷在計(jì)算過(guò)程中會(huì)由數(shù)學(xué)間斷變成數(shù)值過(guò)渡區(qū),按照CFD理論,一維有限差分格式的數(shù)值解在過(guò)渡區(qū)滿(mǎn)足式(7):

(7)

因?yàn)橐粋€(gè)網(wǎng)格節(jié)點(diǎn)無(wú)法存儲(chǔ)兩組數(shù)據(jù),因此初始間斷必占據(jù)兩個(gè)網(wǎng)格節(jié)點(diǎn)。假設(shè)初始間斷的中心在半節(jié)點(diǎn)i+1/2處(又稱(chēng)為界面),那么i點(diǎn)的初始流場(chǎng)參數(shù)和i-1點(diǎn)的相等。

對(duì)于FVS格式,一階迎風(fēng)格式離散Euler方程寫(xiě)成:

(8)

式中,上標(biāo)“+”和“-”分別代表通量為正、負(fù)數(shù)值通量。

田志芳一個(gè)人在地窩子里，看看頭頂，看看腳下，一屁股坐在土臺(tái)上，嘆口氣，心想姆媽拼死阻攔都沒(méi)攔住她和哥哥，現(xiàn)在怪誰(shuí)呢，自己跳起腳要支邊。她垂下頭，把手中的沙棗花捧起來(lái)瞧，帶沙點(diǎn)的葉根處，確實(shí)有細(xì)小的花苞，同樣泛出密密麻麻的沙塵，形如青色的米粒，一粒一粒擠在一起，好似家鄉(xiāng)中秋的桂花。猜想，沙棗花開(kāi)了，是不是真有桂花那樣的千里香？

(9)

展開(kāi)有:

(10)

對(duì)于FDS格式,離散Euler方程后得到:

(11)

以Roe格式為例,采用存儲(chǔ)在網(wǎng)格節(jié)點(diǎn)上的流動(dòng)參數(shù)重構(gòu)界面通量,一階精度的界面通量為:

(12)

其中,廣義系數(shù)矩陣A*是根據(jù)Roe平均公式計(jì)算界面處左右變量以后得到。接觸間斷兩側(cè)參數(shù)滿(mǎn)足條件:ui+1=ui=u*,pi+1=pi=p*,可知界面密度變化量對(duì)動(dòng)量方程和能量方程沒(méi)有影響:

(13)

(14)

對(duì)于混合格式AUSM+,其界面通量可以寫(xiě)為:

Fi+1/2=ci+1/2Mai+1/2Φi+1/2+pi+1/2

(15)

考察i+1點(diǎn)流場(chǎng)參數(shù)更新,當(dāng)流場(chǎng)為靜止流場(chǎng)時(shí),半點(diǎn)馬赫數(shù)Mai+1/2=0,且全場(chǎng)壓力相等,因此有:

(16)

因此,在靜止流場(chǎng)中,使用AUSM+格式計(jì)算接觸間斷可保持流場(chǎng)參數(shù)不變。

(17)

(18)

(19)

綜上所述,Van Leer格式以局部馬赫數(shù)為依據(jù)將對(duì)流通量分為正負(fù)兩部分,在全場(chǎng)超聲速條件下,只有正通量參與流場(chǎng)參數(shù)更新,對(duì)接觸間斷問(wèn)題,雖然間斷兩側(cè)密度不等致使密度參數(shù)隨時(shí)間推進(jìn)發(fā)生變化,產(chǎn)生了密度數(shù)值耗散,但更新后的密度參數(shù)對(duì)速度及壓力參數(shù)無(wú)影響;在靜止流場(chǎng)或流場(chǎng)中存在亞聲速區(qū)域時(shí),由于在通量分裂表達(dá)式中聲速和壓力、密度之間是非線(xiàn)性的,將更新后的密度參數(shù)代入動(dòng)量方程及能量方程,流場(chǎng)的速度及壓力參數(shù)發(fā)生改變,從而產(chǎn)生了速度波動(dòng)及壓力波動(dòng)。

Roe格式計(jì)算接觸間斷時(shí)以密度波推動(dòng)流場(chǎng)參數(shù)更新,同時(shí)受流場(chǎng)速度影響,AUSM+格式計(jì)算接觸間斷時(shí)同樣以密度波推動(dòng)流場(chǎng)參數(shù)更新,同時(shí)受界面馬赫數(shù)(半點(diǎn)馬赫數(shù))影響,因此在靜止流場(chǎng)中使用Roe格式及AUSM+格式時(shí),流場(chǎng)參數(shù)無(wú)更新,自然無(wú)數(shù)值耗散產(chǎn)生;在非靜止流場(chǎng)中使用以上兩種格式時(shí),質(zhì)量方程參數(shù)更新,導(dǎo)致密度參數(shù)發(fā)生變換,產(chǎn)生密度數(shù)值耗散,但將更新后的密度參數(shù)代入動(dòng)量方程及能量方程,速度及壓力參數(shù)無(wú)變化。

同時(shí),由以上公式推導(dǎo)可以看出,當(dāng)初始間斷的中心位于半節(jié)點(diǎn)i+1/2處且全場(chǎng)為超聲速條件時(shí),點(diǎn)i位于間斷的上游區(qū)域,根據(jù)雙曲型方程擾動(dòng)不向上游傳播的特性,點(diǎn)i處流場(chǎng)參數(shù)不應(yīng)隨時(shí)間發(fā)生變化,因此一個(gè)時(shí)間步后,Van Leer格式及AUSM+格式計(jì)算的接觸間斷流場(chǎng),皆為點(diǎn)i+1處的密度參數(shù)發(fā)生變化,點(diǎn)i處流場(chǎng)參數(shù)無(wú)變化;而Roe格式計(jì)算的接觸間斷流場(chǎng),點(diǎn)i處的密度參數(shù)發(fā)生變化,這顯然不符合雙曲型方程擾動(dòng)不向上游傳播的特性,這或許與“紅玉”現(xiàn)象出現(xiàn)的原因有關(guān),有待進(jìn)一步驗(yàn)證。

2.3 二維接觸間斷問(wèn)題

上文給出了三種不同類(lèi)型的迎風(fēng)格式在一維接觸間斷中的數(shù)值計(jì)算結(jié)果,并通過(guò)公式推導(dǎo)對(duì)不同計(jì)算結(jié)果出現(xiàn)的原因進(jìn)行了分析。為了加深對(duì)不同迎風(fēng)格式下接觸間斷的認(rèn)識(shí),給出以下二維算例。

在二維問(wèn)題中,將初始接觸間斷放置于x方向?yàn)槌曀俚木鶆蛄鲌?chǎng)中,設(shè)置計(jì)算區(qū)域?yàn)閇0,2]×[-2,2],網(wǎng)格量為200×400。將間斷放置于y=0處,初始流動(dòng)參數(shù)為:

(20)

其中,ρ1=4,ρ2=1,u=2,v=0,p=1/1.4。左右邊界給定超聲速出入口邊界條件,上下邊界為一階外推。

在本節(jié)中,首先使用一階迎風(fēng)格式進(jìn)行計(jì)算,同樣選用Van Leer、Roe、AUSM+三種迎風(fēng)格式,時(shí)間步長(zhǎng)按CFL=0.5計(jì)算,計(jì)算終止時(shí)刻tstop=1.5。計(jì)算結(jié)果表明,當(dāng)使用Roe格式或AUSM+格式時(shí),流場(chǎng)參數(shù)都能保持初始值不變,因此本節(jié)只給出了使用Van Leer格式的計(jì)算結(jié)果。圖8是空間離散采用五階WENOZ格式[24]的計(jì)算結(jié)果,此時(shí)時(shí)間離散采用具有總變差減小(total variation diminishing,TVD)性質(zhì)的三階Runge-Kutta格式[25]。圖8給出了單接觸間斷下計(jì)算終止時(shí)流場(chǎng)的相對(duì)壓力誤差(δp=1.4Δp)分布云圖及渦量(wz=?v/?x-?u/?y)分布云圖。因?yàn)槌跏剂鲌?chǎng)速度為常數(shù),渦量為0,因此圖8給出的渦量分布云圖也代表了流場(chǎng)的誤差分布。從圖中可以看出,在計(jì)算過(guò)程中Van Leer迎風(fēng)格式不能完好地保持間斷,在間斷兩側(cè)會(huì)有誤差產(chǎn)生,誤差分別以當(dāng)?shù)芈曀傺貀方向向兩側(cè)進(jìn)行傳播,在x方向超聲速氣流的作用下形成以特征線(xiàn)為邊界的誤差分布范圍。

(a) 壓力誤差分布云圖(a) Cloud map of pressure error distribution

(b) 渦量分布云圖(b) Cloud map of vorticity distribution

為進(jìn)一步研究二維接觸間斷對(duì)流場(chǎng)結(jié)構(gòu)的影響,設(shè)計(jì)了同時(shí)存在兩個(gè)間斷的流場(chǎng)算例。初始流場(chǎng)參數(shù)設(shè)置為:

(21)

其中,ρ1=4.0,ρ2=1.0,u=2.0,v=0,p=1/1.4。

圖9給出了雙接觸間斷下,在計(jì)算時(shí)間t=1.5時(shí)的流場(chǎng)壓力誤差分布云圖及渦量分布云圖。從圖中可以看出兩個(gè)間斷所產(chǎn)生的誤差在相互干擾過(guò)程中會(huì)改變每個(gè)間斷的誤差分布范圍,且會(huì)產(chǎn)生復(fù)雜的小尺度結(jié)構(gòu)。

(a) 壓力誤差分布云圖(a) Cloud map of pressure error distribution

(b) 渦量分布云圖(b) Cloud map of vorticity distribution

以上算例可以看成是y方向的一維靜止接觸間斷與x方向超聲速自由流的組合問(wèn)題,因?yàn)閥方向速度為0,此時(shí)Roe格式和AUSM+格式對(duì)接觸間斷數(shù)值耗散為0,表現(xiàn)為流場(chǎng)參數(shù)保持初始值不變。而Van Leer格式在對(duì)接觸間斷有密度數(shù)值耗散的同時(shí),其在間斷從初始數(shù)學(xué)上的間斷變成有厚度的數(shù)值剪切層時(shí)誘導(dǎo)出的數(shù)值誤差相互干擾,會(huì)生成復(fù)雜的非物理小尺度結(jié)構(gòu),從而影響流場(chǎng)的結(jié)構(gòu)分布。

3 線(xiàn)性分布流場(chǎng)

第2節(jié)給出了使用不同迎風(fēng)格式數(shù)值模擬接觸間斷問(wèn)題所得到的計(jì)算結(jié)果,并通過(guò)公式推導(dǎo)對(duì)不同結(jié)果出現(xiàn)的原因進(jìn)行了理論分析。本節(jié)繼續(xù)給了以上三種迎風(fēng)格式在流場(chǎng)參數(shù)線(xiàn)性分布的流場(chǎng)中的計(jì)算結(jié)果。

設(shè)置無(wú)量綱計(jì)算區(qū)域?yàn)閇0,1]×[0,1],網(wǎng)格量為100×100。給定流場(chǎng)密度參數(shù)為線(xiàn)性分布,初始參數(shù)為:

(ρ,u,v,p)=(1+y,0,0,1/1.4)

(22)

邊界給初始值,即理論值。計(jì)算終止時(shí)刻tstop=2.0,時(shí)間步長(zhǎng)按CFL=0.5(通量無(wú)分裂時(shí)CFL=0.1)計(jì)算。

首先使用一階迎風(fēng)格式進(jìn)行計(jì)算,圖10、圖11在給出Van Leer、Roe、AUSM+三種迎風(fēng)格式計(jì)算結(jié)果的同時(shí),給出了對(duì)流通量無(wú)分裂時(shí)的計(jì)算結(jié)果。計(jì)算結(jié)果表明,在密度成線(xiàn)性分布的靜止流場(chǎng)中,Roe、AUSM+格式及通量無(wú)分裂計(jì)算方法都能保持流場(chǎng)不變,而Van Leer格式的計(jì)算結(jié)果有較大誤差產(chǎn)生,且該誤差隨計(jì)算時(shí)間增加而逐漸增大。

(a) Van Leer

(b) Roe

(c) AUSM+

(d) 通量無(wú)分裂 (d) Flux non splitting

圖11 中心點(diǎn)處密度絕對(duì)誤差隨時(shí)間變化曲線(xiàn)Fig.11 Time variation curves of density absolute error at central point

考慮到所用空間離散格式為一階迎風(fēng)格式,而此時(shí)流場(chǎng)為二階精度流場(chǎng),使用一階格式計(jì)算二階流場(chǎng)可能會(huì)產(chǎn)生誤差,因此進(jìn)一步使用二階迎風(fēng)MUSCL格式、五階WCNS格式[26]和五階WENO格式[27]對(duì)流場(chǎng)進(jìn)行了數(shù)值模擬,計(jì)算結(jié)果如圖12所示。

(a) 二階迎風(fēng)MUSCL格式(a) 2nd-order upwind MUSCL scheme

(b) 五階WCNS格式(b) 5th-order WCNS scheme

(c) 五階WENO格式(c) 5th-order WENO scheme

從以上計(jì)算結(jié)果可以看出,Roe格式和AUSM+格式能夠保持初始流場(chǎng)參數(shù)的線(xiàn)性分布不變,而Van Leer格式在空間重構(gòu)對(duì)象為原始變量時(shí)可以保持流場(chǎng)參數(shù)線(xiàn)性分布不變,在空間重構(gòu)對(duì)象為對(duì)流通量時(shí)則會(huì)有誤差產(chǎn)生,破壞流場(chǎng)的結(jié)構(gòu)分布。

4 結(jié)論

對(duì)三種迎風(fēng)格式在接觸間斷中的數(shù)值耗散問(wèn)題進(jìn)行了數(shù)值實(shí)驗(yàn),并通過(guò)公式推導(dǎo)對(duì)不同流場(chǎng)參數(shù)下數(shù)值耗散產(chǎn)生的機(jī)理進(jìn)行了分析。以Roe格式為代表的FDS格式及以AUSM+格式為代表的混合格式,在接觸間斷問(wèn)題中存在密度數(shù)值耗散,其耗散受密度差推動(dòng)產(chǎn)生,同時(shí)受流場(chǎng)速度(馬赫數(shù))影響,但因?yàn)榇藭r(shí)其質(zhì)量方程和動(dòng)量方程、能量方程為解耦關(guān)系,所以更新后的密度參數(shù)對(duì)速度及壓力參數(shù)無(wú)影響。以Van Leer格式為代表的FVS格式在計(jì)算接觸間斷問(wèn)題時(shí),不僅存在密度數(shù)值耗散,在流場(chǎng)靜止或流場(chǎng)內(nèi)存在亞聲速區(qū)域條件下,密度耗散的產(chǎn)生還會(huì)誘導(dǎo)出速度擾動(dòng)誤差及壓力擾動(dòng)誤差,該擾動(dòng)誤差對(duì)格式的數(shù)值耗散大小無(wú)影響,但對(duì)流場(chǎng)結(jié)構(gòu)的影響無(wú)法忽略。特別地,對(duì)L-F分裂格式,局部L-F分裂會(huì)出現(xiàn)以上非物理振蕩誤差,全局L-F分裂則不會(huì)產(chǎn)生非物理振蕩現(xiàn)象。因此,文章中的FVS格式不包括全局L-F分裂格式。

FVS格式在計(jì)算一維接觸間斷問(wèn)題時(shí)產(chǎn)生的誘導(dǎo)數(shù)值誤差在二維接觸間斷問(wèn)題中表現(xiàn)為復(fù)雜小尺度結(jié)構(gòu),若二維接觸間斷存在于復(fù)雜結(jié)構(gòu)流場(chǎng)中,復(fù)雜小尺度結(jié)構(gòu)的產(chǎn)生必會(huì)對(duì)流場(chǎng)結(jié)構(gòu)分析帶來(lái)困難。此外在空間重構(gòu)對(duì)象為對(duì)流通量時(shí),使用FVS格式計(jì)算密度線(xiàn)性分布的流場(chǎng),會(huì)破壞流場(chǎng)參數(shù)原有的梯度,產(chǎn)生較大的數(shù)值誤差,在高階格式條件下,誤差大大減小,但并未消除,使得流場(chǎng)的計(jì)算精度難以達(dá)到二階。在空間重構(gòu)對(duì)象為原始變量時(shí),使用FVS格式可以很好地保持流場(chǎng)梯度。