迎風格式在接觸間斷的數值耗散及其誘導誤差

韓 芳,魏雁昕,劉 君

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

在計算流體力學(computational fluid dynamics, CFD)中,數值耗散主要來自空間離散格式離散對流項和時間離散格式離散時間項所產生的截斷誤差,其大小與格式的計算精度及穩定性息息相關。因此,CFD研究者們常常采用降低數值耗散的方法來構造更高精度的計算格式[1-3]。空間離散格式包含差分格式及對流通量導數的求解兩部分,本文主要討論對流通量導數計算方法的數值耗散特性,而迎風格式作為目前CFD求解對流通量的主流格式[4],成為本文的研究對象。

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

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

為了進一步分析不同迎風格式數值耗散的產生機理,同時認清上文所提到的非物理波動是否與格式的數值耗散有關,采用具有單一流場結構的接觸間斷作為數值模擬對象開展研究。綜合不同迎風格式的數值耗散及其非物理波動在接觸間斷問題中的表現,可以為不同的流場問題選擇更加合適的計算格式。

1 控制方程與數值方法

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

(1)

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

(2)

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

(3)

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

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

(4)

以上公式及變量皆為無量綱化的,下文數值實驗中的所有變量和數據也是無量綱化的。

文中以迎風格式作為研究對象,分別選用Van Leer、Roe、AUSM+格式作為三種迎風格式的代表進行數值計算。

2 接觸間斷問題

2.1 一維接觸間斷問題

使用二維均勻正交網格計算一維流場,計算區域設定為[0,10]×[0,1],網格量為200×20。無量綱化的初始條件為:

(5)

其中,xd為間斷的初始位置。將xd左側的區域記為1區,將xd右側的區域記為2區,因此ρ1為間斷左側的初始密度參數。

對一維接觸間斷問題,由之后的計算結果可以看出,在初始間斷形成數值過渡層的過程中出現了非物理波動,其傳播速度為特征速度,為避免擾動傳出邊界,設置間斷的初始位置xd滿足式(6):

(6)

其中,c1為1區的聲速,tstop為計算終止時刻。根據初始間斷左右流場馬赫數的不同,即流場初始條件的不同,設置以下五個算例,其具體參數如表1所示。

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

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

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

從圖1可以看出,在靜止流場中,Roe格式及AUSM+格式可以完整地保持接觸間斷不變,其密度、速度及壓力參數分布曲線皆保持初始狀態。而Van Leer格式計算的流場參數分布曲線不僅有密度耗散,速度及壓力曲線也有波動。

(a) 密度參數分布曲線(a) Density parameter distribution curves

(b) 速度參數分布曲線(b) Velocity parameter distribution curves

(c) 壓力參數分布曲線(c) Pressure parameter distribution curves

從圖2中可以看出,在流場為全場超聲速條件時,三種格式對接觸間斷的數值模擬結果相同,都是僅有密度的耗散,速度及壓力參數無變化。

(a) 密度參數分布曲線(a) Density parameter distribution curves

(b) 速度參數分布曲線(b) Velocity parameter distribution curves

(c) 壓力參數分布曲線(c) Pressure parameter distribution curves

從圖3～5可以看出,在流場中存在亞聲速區域時,Roe格式及AUSM+格式捕捉的接觸間斷僅有密度的變化,速度及壓力參數無變化。Van Leer格式捕捉的接觸間斷存在速度及壓力參數的變化,其波動曲線與在靜止流場中的結果相似,不僅在間斷處有波動,在間斷兩側也產生了兩個非物理波動。但是,盡管Van Leer格式計算的接觸間斷有速度及壓力的變化,其密度變化曲線與Roe格式及AUSM+格式的計算結果曲線幾乎重合。

根據以上算例的計算結果,對一維接觸間斷,Roe格式及AUSM+格式在靜止流場中沒有耗散,在非靜止流場中有密度耗散;Van Leer格式則在所有流場中都有密度的耗散。圖6給出了不同時刻靜止流場(初始條件1)中非物理波動的位置變化曲線,從圖中可以看出,在上文中出現的非物理波動分別以特征速度u-c1和u+c2逐步遠離間斷。將間斷左側波記為“波1”,間斷右側波記為“波2”,則圖中v1、v2分別代表兩個波動的傳播速度,c1、c2分別為當地點聲速。

圖7給出了使用初始條件1計算的接觸間斷流場壓力波動幅值隨時間變化曲線,從圖中可以看出,隨著時間步的推進,誘導誤差在遠離接觸間斷的同時,其幅值大小也會逐漸降低,但不會消失。結合圖6和圖7可以看出,對定常問題,在計算時間足夠長的情況下,誤差會運動出有限的計算區域,流場結構會恢復原本的狀態。但對非定常問題,波動對流場的影響無法消除,在參考文獻[16-18]的算例中可以看到波動對間斷下游流場參數的污染。

此外,這逐步遠離間斷的非物理波動不應歸于Van Leer格式的數值耗散,而是一種數值誤差,因為數值耗散的作用是“抹平”間斷,而這兩個非物理波動對密度曲線的過渡區并無影響。

(a) 密度參數分布曲線(a) Density parameter distribution curves

(b) 速度參數分布曲線(b) Velocity parameter distribution curves

(c) 壓力參數分布曲線(c) Pressure parameter distribution curves

(a) 密度參數分布曲線(a) Density parameter distribution curves

(b) 速度參數分布曲線(b) Velocity parameter distribution curves

(c) 壓力參數分布曲線(c) Pressure parameter distribution curves

(a) 密度參數分布曲線(a) Density parameter distribution curves

(b) 速度參數分布曲線(b) Velocity parameter distribution curves

(c) 壓力參數分布曲線(c) Pressure parameter distribution curves

圖6 波動位置隨時間變化曲線(初始條件1)Fig.6 Curve of fluctuation position with time(initial condition 1)

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

2.2 對數值計算結果的理論推導及機理分析

滿足Euler方程的初始間斷在計算過程中會由數學間斷變成數值過渡區,按照CFD理論,一維有限差分格式的數值解在過渡區滿足式(7):

(7)

因為一個網格節點無法存儲兩組數據,因此初始間斷必占據兩個網格節點。假設初始間斷的中心在半節點i+1/2處(又稱為界面),那么i點的初始流場參數和i-1點的相等。

對于FVS格式,一階迎風格式離散Euler方程寫成:

(8)

式中,上標“+”和“-”分別代表通量為正、負數值通量。

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

(9)

展開有:

(10)

對于FDS格式,離散Euler方程后得到:

(11)

以Roe格式為例,采用存儲在網格節點上的流動參數重構界面通量,一階精度的界面通量為:

(12)

其中,廣義系數矩陣A*是根據Roe平均公式計算界面處左右變量以后得到。接觸間斷兩側參數滿足條件:ui+1=ui=u*,pi+1=pi=p*,可知界面密度變化量對動量方程和能量方程沒有影響:

(13)

(14)

對于混合格式AUSM+,其界面通量可以寫為:

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

(15)

考察i+1點流場參數更新,當流場為靜止流場時,半點馬赫數Mai+1/2=0,且全場壓力相等,因此有:

(16)

因此,在靜止流場中,使用AUSM+格式計算接觸間斷可保持流場參數不變。

(17)

(18)

(19)

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

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

同時,由以上公式推導可以看出,當初始間斷的中心位于半節點i+1/2處且全場為超聲速條件時,點i位于間斷的上游區域,根據雙曲型方程擾動不向上游傳播的特性,點i處流場參數不應隨時間發生變化,因此一個時間步后,Van Leer格式及AUSM+格式計算的接觸間斷流場,皆為點i+1處的密度參數發生變化,點i處流場參數無變化;而Roe格式計算的接觸間斷流場,點i處的密度參數發生變化,這顯然不符合雙曲型方程擾動不向上游傳播的特性,這或許與“紅玉”現象出現的原因有關,有待進一步驗證。

2.3 二維接觸間斷問題

上文給出了三種不同類型的迎風格式在一維接觸間斷中的數值計算結果,并通過公式推導對不同計算結果出現的原因進行了分析。為了加深對不同迎風格式下接觸間斷的認識,給出以下二維算例。

在二維問題中,將初始接觸間斷放置于x方向為超聲速的均勻流場中,設置計算區域為[0,2]×[-2,2],網格量為200×400。將間斷放置于y=0處,初始流動參數為:

(20)

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

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

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

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

為進一步研究二維接觸間斷對流場結構的影響,設計了同時存在兩個間斷的流場算例。初始流場參數設置為:

(21)

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

圖9給出了雙接觸間斷下,在計算時間t=1.5時的流場壓力誤差分布云圖及渦量分布云圖。從圖中可以看出兩個間斷所產生的誤差在相互干擾過程中會改變每個間斷的誤差分布范圍,且會產生復雜的小尺度結構。

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

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

以上算例可以看成是y方向的一維靜止接觸間斷與x方向超聲速自由流的組合問題,因為y方向速度為0,此時Roe格式和AUSM+格式對接觸間斷數值耗散為0,表現為流場參數保持初始值不變。而Van Leer格式在對接觸間斷有密度數值耗散的同時,其在間斷從初始數學上的間斷變成有厚度的數值剪切層時誘導出的數值誤差相互干擾,會生成復雜的非物理小尺度結構,從而影響流場的結構分布。

3 線性分布流場

第2節給出了使用不同迎風格式數值模擬接觸間斷問題所得到的計算結果,并通過公式推導對不同結果出現的原因進行了理論分析。本節繼續給了以上三種迎風格式在流場參數線性分布的流場中的計算結果。

設置無量綱計算區域為[0,1]×[0,1],網格量為100×100。給定流場密度參數為線性分布,初始參數為:

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

(22)

邊界給初始值,即理論值。計算終止時刻tstop=2.0,時間步長按CFL=0.5(通量無分裂時CFL=0.1)計算。

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

(a) Van Leer

(b) Roe

(c) AUSM+

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

圖11 中心點處密度絕對誤差隨時間變化曲線Fig.11 Time variation curves of density absolute error at central point

考慮到所用空間離散格式為一階迎風格式,而此時流場為二階精度流場,使用一階格式計算二階流場可能會產生誤差,因此進一步使用二階迎風MUSCL格式、五階WCNS格式[26]和五階WENO格式[27]對流場進行了數值模擬,計算結果如圖12所示。

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

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

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

從以上計算結果可以看出,Roe格式和AUSM+格式能夠保持初始流場參數的線性分布不變,而Van Leer格式在空間重構對象為原始變量時可以保持流場參數線性分布不變,在空間重構對象為對流通量時則會有誤差產生,破壞流場的結構分布。

4 結論

對三種迎風格式在接觸間斷中的數值耗散問題進行了數值實驗,并通過公式推導對不同流場參數下數值耗散產生的機理進行了分析。以Roe格式為代表的FDS格式及以AUSM+格式為代表的混合格式,在接觸間斷問題中存在密度數值耗散,其耗散受密度差推動產生,同時受流場速度(馬赫數)影響,但因為此時其質量方程和動量方程、能量方程為解耦關系,所以更新后的密度參數對速度及壓力參數無影響。以Van Leer格式為代表的FVS格式在計算接觸間斷問題時,不僅存在密度數值耗散,在流場靜止或流場內存在亞聲速區域條件下,密度耗散的產生還會誘導出速度擾動誤差及壓力擾動誤差,該擾動誤差對格式的數值耗散大小無影響,但對流場結構的影響無法忽略。特別地,對L-F分裂格式,局部L-F分裂會出現以上非物理振蕩誤差,全局L-F分裂則不會產生非物理振蕩現象。因此,文章中的FVS格式不包括全局L-F分裂格式。

FVS格式在計算一維接觸間斷問題時產生的誘導數值誤差在二維接觸間斷問題中表現為復雜小尺度結構,若二維接觸間斷存在于復雜結構流場中,復雜小尺度結構的產生必會對流場結構分析帶來困難。此外在空間重構對象為對流通量時,使用FVS格式計算密度線性分布的流場,會破壞流場參數原有的梯度,產生較大的數值誤差,在高階格式條件下,誤差大大減小,但并未消除,使得流場的計算精度難以達到二階。在空間重構對象為原始變量時,使用FVS格式可以很好地保持流場梯度。