張 帆,張力丹,劉 君,陳飆松

(大連理工大學(xué) 工業(yè)裝備結(jié)構(gòu)分析國(guó)家重點(diǎn)實(shí)驗(yàn)室,大連 116024)

0 引 言

迎風(fēng)格式被廣泛用于可壓縮流計(jì)算流體力學(xué)的有限差分法和有限體積法中守恒型方程的通量計(jì)算。由于迎風(fēng)格式固有的耗散性質(zhì),它能夠有效地進(jìn)行包含激波等流場(chǎng)間斷的超聲速、高超聲速流動(dòng)的計(jì)算,這也就是激波捕捉法;而迎風(fēng)格式對(duì)應(yīng)的解析精確Jacobian矩陣也是提高隱式時(shí)間離散格式收斂性和穩(wěn)定性的關(guān)鍵[1-2]。但是,部分迎風(fēng)格式在含激波的多維流場(chǎng)中出現(xiàn)數(shù)值激波不穩(wěn)定現(xiàn)象,這一問(wèn)題在高超聲速計(jì)算中尤為突出[3],并對(duì)計(jì)算精度與穩(wěn)定性都造成了不利影響。

分析表明,經(jīng)典的通量矢量分裂(Flux Vector Splitting,FVS)型迎風(fēng)格式計(jì)算效率高,穩(wěn)定性好,能夠可靠地計(jì)算激波。但是,經(jīng)典的FVS格式不能夠準(zhǔn)確模擬線性波,導(dǎo)致流場(chǎng)中的接觸間斷[4]被抹平,進(jìn)而影響摩阻或熱流等黏性效應(yīng)主導(dǎo)的問(wèn)題的計(jì)算精度[5]。而通量差分分裂(Flux Difference Splitting,FDS)型格式能夠有效地計(jì)算定常激波與接觸間斷,但是需要熵修正[6-7]或存在數(shù)值激波不穩(wěn)定[8]等問(wèn)題。

激波捕捉法具有使用靈活與編程簡(jiǎn)易的優(yōu)點(diǎn),因此,為了改善迎風(fēng)格式的性能,國(guó)內(nèi)外都進(jìn)行了大量研究工作。其中,AUSM格式[9]將流動(dòng)中的對(duì)流機(jī)制與壓力傳播機(jī)制分別進(jìn)行處理,給出了新的格式構(gòu)造思路。AUSM格式結(jié)合了FVS與FDS的特點(diǎn),具有較好的綜合性能。但是AUSM格式及其后續(xù)發(fā)展的同類格式也沒(méi)有完全解決前述的激波計(jì)算不穩(wěn)定問(wèn)題。考慮到實(shí)際流動(dòng)的多維效應(yīng),多種旋轉(zhuǎn)迎風(fēng)格式被構(gòu)造出來(lái)[10-13],其機(jī)理也得到進(jìn)一步研究[14-15]。旋轉(zhuǎn)格式能夠有效地抑制激波不穩(wěn)定現(xiàn)象,但是由于對(duì)兩個(gè)方向進(jìn)行計(jì)算,導(dǎo)致計(jì)算量較大;并且由于旋轉(zhuǎn)通量函數(shù)的非線性性質(zhì),為保證收斂性需要附加處理[10-12]?；诮Y(jié)合不同格式優(yōu)點(diǎn)的思路,混合格式被提出并用于改善迎風(fēng)格式性能[16-20]。研究表明,捕捉激波的穩(wěn)定性與接觸間斷的分辨率是一對(duì)矛盾,單一的迎風(fēng)格式難以同時(shí)在這兩個(gè)方面具有良好的性能,因而需要采用耗散可調(diào)的旋轉(zhuǎn)格式或混合格式。

由于非結(jié)構(gòu)網(wǎng)格上無(wú)法保證正交性,同時(shí)網(wǎng)格單元與流動(dòng)間斷交錯(cuò),導(dǎo)致非結(jié)構(gòu)網(wǎng)格上的計(jì)算更易于產(chǎn)生非物理擾動(dòng)。本文對(duì)非結(jié)構(gòu)網(wǎng)格上的數(shù)值激波不穩(wěn)定現(xiàn)象進(jìn)行了研究,分析了不穩(wěn)定現(xiàn)象與網(wǎng)格分布的關(guān)系;在此基礎(chǔ)上,根據(jù)迎風(fēng)型格式的耗散機(jī)制,研究并構(gòu)造新型混合型格式;最后,通過(guò)一系列數(shù)值算例對(duì)比研究新型混合格式的性能。

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

1.1 控制方程及其離散

本文控制方程采用積分形式的二維可壓縮Navier-Stokes(NS)方程組:

式中,Ω為積分域,?Ω為積分域邊界,n為積分域邊界外法向;Q為方程守恒變量,Fc(Q)和Fv(Q)分別為對(duì)流通量與黏性通量。關(guān)于可壓縮N-S方程的詳細(xì)推導(dǎo)與說(shuō)明可以參考文獻(xiàn)[21]。

本文采用格心型有限體積法對(duì)上述方程進(jìn)行離散求解,即以網(wǎng)格單元為空間離散后的積分域/控制體,從而得到方程組半離散形式如下:

其中:下標(biāo)i表示第i個(gè)單元,Nf為該單元的表面的個(gè)數(shù)。Fk、nk和Sk分別為該單元第k個(gè)表面的通量、外法向單位矢量和面積。第i個(gè)單元內(nèi)的變量分布可以描述為:

此處下標(biāo)i表示變量定義于單元i的形心。?i為限制器值,?Qi為單元內(nèi)的變量梯度,Δr為相對(duì)單元形心的矢徑。采用式(3)描述變量分布時(shí),有限體積法的空間離散具有二階精度;當(dāng)變量梯度為零時(shí),空間離散精度降為一階。

完成空間離散求解后,可以采用顯式或隱式方法對(duì)式(2)中的時(shí)間導(dǎo)數(shù)項(xiàng)進(jìn)行離散求解。本文采用的方法包括顯式一階前差與隱式LU-SGS方法[22]。

1.2 通量計(jì)算

式(2)中的單元表面通量包括了對(duì)流通量與黏性通量。其中,黏性通量采用Blazek的方法[21]構(gòu)造面元梯度,而后采用中心差分計(jì)算通量。對(duì)流通量可以采用中心型格式或迎風(fēng)型格式計(jì)算,本文僅討論迎風(fēng)型格式。

采用迎風(fēng)格式計(jì)算單元表面對(duì)流通量,可以寫(xiě)為以下統(tǒng)一形式[9]:

其中上標(biāo)＋、－分別表示面元左右兩側(cè)變量。為簡(jiǎn)便,后續(xù)討論中省略上式中的面元通量下標(biāo)。在上式的形式下,對(duì)流通量被分為對(duì)流項(xiàng)Φ±與壓力項(xiàng)并分別進(jìn)行計(jì)算。本文主要對(duì)比三種對(duì)流迎風(fēng)分裂型格式,其中包括屬于AUSM類型的SLAU格式[23],以及將能量方程中壓力相關(guān)信息全部由壓力項(xiàng)p~處理的TV格式[24]和K-CUSP-X格式[25]。這三個(gè)迎風(fēng)型格式均可以寫(xiě)為式(4)的形式,但是各項(xiàng)的含義不同,將其差異簡(jiǎn)要介紹如下。

對(duì)于SLAU格式,其分裂的對(duì)流項(xiàng)與壓力項(xiàng)分別為:

式中:ρ、u、v、h分別為密度、x方向速度、y方向速度以及總焓。是格式對(duì)壓力項(xiàng)進(jìn)行近似描述的標(biāo)量函數(shù)。由上式可見(jiàn),AUSM系列格式認(rèn)為對(duì)流項(xiàng)決定了能量通量,而壓力項(xiàng)對(duì)能量方程沒(méi)有直接影響。對(duì)于TV格式與K-CUSP-X格式,能量方程中與壓力相關(guān)的部分被移到壓力項(xiàng)中,如式(6)所示:

本節(jié)僅簡(jiǎn)要介紹三種格式間的主要差異,格式實(shí)現(xiàn)細(xì)節(jié)限于篇幅不再介紹。數(shù)值計(jì)算表明,上述三種格式均具有低耗散性,能夠準(zhǔn)確地模擬接觸間斷以及黏性流動(dòng)。此外,SLAU格式在低馬赫數(shù)流動(dòng)中耗散較小;K-CUSP-X格式具有速度擾動(dòng)衰減特性,相對(duì)TV格式提高了激波計(jì)算的穩(wěn)定性。

1.3 空間變量重構(gòu)與限制器

如前所述,采用式(3)描述變量分布時(shí),空間離散具有二階精度。本文采用Gauss-Green積分方法[26]計(jì)算單元內(nèi)的梯度,其中積分所需格點(diǎn)變量可以通過(guò)多種方法計(jì)算[27],在此僅采用加權(quán)平均法。需要注意的是,在存在間斷的流場(chǎng)中,二階變量重構(gòu)可能引起振蕩,需要采用限制器確保重構(gòu)變量的單調(diào)性。尤其在高超聲速流動(dòng)計(jì)算中,限制器對(duì)計(jì)算精度、穩(wěn)定性等都有重要意義。本文在含間斷流場(chǎng)計(jì)算的二階格式中采用與Venkatakrishnan限制器[28]相結(jié)合的非結(jié)構(gòu)網(wǎng)格MLP(Multi-dimensional Limiting Process)限制器[29]及其最新改進(jìn)型MLP-pw限制器[30]。非結(jié)構(gòu)網(wǎng)格的MLP限制器應(yīng)用了Venkatakrishnan限制器的限制函數(shù),通過(guò)協(xié)調(diào)考慮單元格心與格點(diǎn)變量分布,提出MLP條件,從而利用更多的單元模板信息以便更好地近似多維流動(dòng)變量,并抑制非物理振蕩。而MLP-pw限制器構(gòu)造了強(qiáng)化和放寬的兩個(gè)MLP型限制條件,分別用于提升激波捕捉的穩(wěn)定性與降低光滑流場(chǎng)的耗散。

2 數(shù)值激波不穩(wěn)定現(xiàn)象

本節(jié)采用1.2節(jié)中的三種迎風(fēng)格式計(jì)算高超聲速鈍體無(wú)黏繞流,測(cè)試它們?cè)诜墙Y(jié)構(gòu)網(wǎng)格上模擬強(qiáng)激波的穩(wěn)定性,并且給出了應(yīng)用廣泛的Roe格式[31]和HLLC格式[32]的結(jié)果作為參考。為此生成了兩種類型的三角形網(wǎng)格。其中,生成規(guī)則三角形網(wǎng)格時(shí),首先將計(jì)算區(qū)域剖分為120×80的四邊形網(wǎng)格,而后對(duì)每個(gè)四邊形單元沿對(duì)角線剖分;不規(guī)則三角形網(wǎng)格采用Delaunay三角化方法將流場(chǎng)剖分為14 946個(gè)不規(guī)則分布的三角形單元,在激波附近與波后區(qū)域,最大和最小單元間的面積比小于5,單元長(zhǎng)寬比均小于2。圖1中給出網(wǎng)格的示例。為了顯示的清晰,圖中網(wǎng)格邊界的離散點(diǎn)只有實(shí)際計(jì)算網(wǎng)格的一半,內(nèi)部網(wǎng)格采用與實(shí)際計(jì)算網(wǎng)格相同的方法生成。

圖1 測(cè)試網(wǎng)格示例Fig.1 Schematics of computational grids

該問(wèn)題的計(jì)算中,設(shè)定來(lái)流馬赫數(shù)為8的單一組分均勻來(lái)流,流動(dòng)受半圓形鈍體阻礙形成弓形激波。本節(jié)中均采用一階空間精度進(jìn)行計(jì)算。為了突出問(wèn)題,首先討論不規(guī)則網(wǎng)格的計(jì)算結(jié)果。這里采用Roe格式、HLLC格式和TV格式計(jì)算,都出現(xiàn)了典型的“Carbuncle”現(xiàn)象。而后,采用SLAU格式和KCUSP-X格式分別在TV格式得到的流場(chǎng)基礎(chǔ)上繼續(xù)計(jì)算,仍然存在振蕩變化的異常現(xiàn)象。圖2(f)中給出K-CUSP-X格式異常區(qū)域附近的流線,出現(xiàn)了明顯的分離渦。需要說(shuō)明,K-CUSP-X在均勻初場(chǎng)條件下能夠得到僅有微弱波動(dòng)的流場(chǎng),而SLAU格式在均勻初場(chǎng)或TV格式計(jì)算結(jié)果基礎(chǔ)上繼續(xù)計(jì)算都表現(xiàn)出類似于圖2(d)的結(jié)果。

由于網(wǎng)格分布不規(guī)則,難以尋找不穩(wěn)定的非物理解與網(wǎng)格或流動(dòng)之間的關(guān)系。因此在規(guī)則網(wǎng)格上進(jìn)行計(jì)算,其它計(jì)算條件保持不變,采用均勻來(lái)流初始化流場(chǎng)。圖3給出SLAU格式在規(guī)則三角形網(wǎng)格上的計(jì)算結(jié)果。合理的流場(chǎng)分布應(yīng)該上下對(duì)稱,但是由于網(wǎng)格分布不同,流場(chǎng)下半?yún)^(qū)出現(xiàn)了非物理解;反之,在上半?yún)^(qū)域則獲得了合理的計(jì)算結(jié)果。限于篇幅,其它格式的計(jì)算結(jié)果并未給出,但均出現(xiàn)不同程度的不穩(wěn)定現(xiàn)象。TV格式的異?，F(xiàn)象更為嚴(yán)重,對(duì)流場(chǎng)整體分布的影響更加顯著;而K-CUSP-X格式的異?，F(xiàn)象影響區(qū)域則較小。

圖2 激波不穩(wěn)定現(xiàn)象(密度云圖)Fig.2 Shock instabilities(density contour)

從圖3(c)中可以看出,出現(xiàn)激波計(jì)算異常結(jié)果區(qū)域中,流線與網(wǎng)格線幾乎相切。關(guān)于數(shù)值激波不穩(wěn)定現(xiàn)象,文獻(xiàn)[33]認(rèn)為具有良好接觸間斷分辨率的格式將出現(xiàn)不穩(wěn)定現(xiàn)象,文獻(xiàn)[34]中指出質(zhì)量通量與網(wǎng)格切向動(dòng)量通量的作用,文獻(xiàn)[35]進(jìn)一步強(qiáng)調(diào)穿過(guò)激波的動(dòng)量擾動(dòng)將導(dǎo)致激波不穩(wěn)定的發(fā)生。在本算例中,網(wǎng)格穿過(guò)激波,兩側(cè)單元的動(dòng)量變化劇烈;網(wǎng)格與流線相切又使擾動(dòng)量沿網(wǎng)格方向傳播。而SLAU格式對(duì)于接觸間斷的低耗散性質(zhì)使動(dòng)量擾動(dòng)難以被耗散,導(dǎo)致出現(xiàn)激波計(jì)算的異常。反之,若采用van Leer格式[36]等無(wú)法準(zhǔn)確捕捉接觸間斷的迎風(fēng)格式,由于包括面元切向速度在內(nèi)的接觸間斷被逐漸抹平,動(dòng)量擾動(dòng)也被抑制,因此不會(huì)在本算例中出現(xiàn)異常結(jié)果。值得注意的是,K-CUSP-X格式在設(shè)計(jì)時(shí)考慮了速度擾動(dòng)的耗散特性,因此其穩(wěn)定性相對(duì)SLAU格式和TV格式更好;但是它對(duì)密度擾動(dòng)并不具備耗散能力,因此并不能完全抑制動(dòng)量擾動(dòng)并避免異?，F(xiàn)象。本節(jié)結(jié)論進(jìn)一步佐證了文獻(xiàn)[33]的結(jié)論,即接觸間斷分辨率良好的格式傾向于出現(xiàn)激波不穩(wěn)定現(xiàn)象。因此,盡管上述格式在多種問(wèn)題中體現(xiàn)出良好的性能,其穩(wěn)定性仍需要進(jìn)一步提高,以改善其對(duì)復(fù)雜網(wǎng)格條件的適應(yīng)能力。

圖3 SLAU格式在規(guī)則三角形網(wǎng)格上的不穩(wěn)定現(xiàn)象[20]Fig.3 Shock instabilities of SLAU scheme on the regular triangular grid

3 混合格式

鑒于迎風(fēng)格式的接觸間斷分辨率與激波穩(wěn)定性存在矛盾,因而組合具有不同特點(diǎn)的迎風(fēng)格式,即設(shè)計(jì)混合迎風(fēng)格式以改進(jìn)迎風(fēng)格式的綜合性能。同時(shí),考慮到動(dòng)量分布的擾動(dòng)被認(rèn)為是激波不穩(wěn)定現(xiàn)象的原因,混合格式的構(gòu)造應(yīng)當(dāng)修正與動(dòng)量相關(guān)的項(xiàng)。

以TV格式為例,借鑒對(duì)流迎風(fēng)分裂格式的構(gòu)造思路,將其寫(xiě)為式(4)的形式后,得到以下各項(xiàng):

其中a為聲速。可以發(fā)現(xiàn)TV格式在進(jìn)行對(duì)流項(xiàng)與壓力項(xiàng)的分裂后,其壓力項(xiàng)主要由壓力而不是速度決定,而壓力場(chǎng)的橢圓型方程性質(zhì)決定它的作用形式與對(duì)流項(xiàng)的雙曲型方程擾動(dòng)形式不同,因而認(rèn)為壓力項(xiàng)對(duì)激波穩(wěn)定性的影響將是次要的?；谶@一考慮,下面的混合格式將不考慮壓力項(xiàng)的修正。

為了簡(jiǎn)化處理,將TV格式在式(4)中的對(duì)流項(xiàng)進(jìn)行合并,可以恢復(fù)其在文獻(xiàn)[24]中的形式:

由此可見(jiàn),TV格式中對(duì)流項(xiàng)計(jì)算的關(guān)鍵是得到合理的速度近似,準(zhǔn)確的速度近似使其獲得了良好的接觸間斷分辨率[37],同時(shí)也導(dǎo)致了激波計(jì)算穩(wěn)定性的下降[33]。

Van Leer格式具有良好的激波計(jì)算穩(wěn)定性,文獻(xiàn)[37]分析了van Leer格式對(duì)接觸間斷的耗散機(jī)制?；旌细袷浇梃bvan Leer格式的馬赫數(shù)分裂函數(shù)[36],并通過(guò)壓力權(quán)函數(shù)得到混合格式的對(duì)流項(xiàng)如下:

其中Ma±為單元表面左右兩側(cè)法向馬赫數(shù)。需要強(qiáng)調(diào)的是,不僅壓力項(xiàng)的處理保持TV格式原有形式,上式中的加權(quán)混合僅對(duì)動(dòng)量方程進(jìn)行,而質(zhì)量方程與能量方程的對(duì)流項(xiàng)仍保持TV格式原有形式。最后定義上式中權(quán)系數(shù)如下[20]:

式(11)中的權(quán)系數(shù)通過(guò)計(jì)算待求通量所在的面元左右兩側(cè)單元的所有表面上的壓力差得到,如圖4所示,k(i)和k(j)分別為單元i和j的所有表面。從而確保跨越激波的網(wǎng)格面能夠隨著壓力差的增大切換至耗散格式。反之,在壓力差較小的區(qū)域,如均勻流場(chǎng)或接觸間斷附近,混合格式完全保持TV格式的原有形式,從而保持其低耗散、高分辨率的性質(zhì)。

圖4 混合格式權(quán)系數(shù)定義示意Fig.4 Schematic of the definition of the weights in the hybrid scheme

4 算例與分析

4.1 亞聲速層流邊界層

首先,為了證明混合格式引入van Leer格式的馬赫數(shù)分裂函數(shù)并未降低黏性流動(dòng)的計(jì)算精度,采用亞聲速層流平板邊界層算例考核混合格式[12]。此時(shí),混合格式的接觸間斷分辨率直接影響邊界層流動(dòng)的計(jì)算精度。均勻來(lái)流馬赫數(shù)為0.3,流場(chǎng)采用圖5所示的四邊形網(wǎng)格離散,除平板前緣外,保證速度邊界層內(nèi)分布10~20層網(wǎng)格。

圖5 層流平板邊界層流動(dòng)計(jì)算網(wǎng)格Fig.5 Computation grid of flat plate boundary layer flow

計(jì)算中采用二階格式,由于流場(chǎng)中沒(méi)有間斷,因此不采用限制器;采用LU-SGS隱式方法計(jì)算至殘差下降8個(gè)量級(jí)后認(rèn)為計(jì)算收斂。給出van Leer格式、TV格式以及混合格式的計(jì)算結(jié)果與Blasius解的對(duì)比,如圖6所示。

圖6 層流平板邊界層結(jié)果對(duì)比Fig.6 Comparison of the results of flat plate boundary layer

可以看到,混合格式和TV格式都獲得了準(zhǔn)確的邊界層速度型分布。這是由于邊界層中壓力梯度較小,引入的修正耗散項(xiàng)的影響可以忽略不計(jì)。該結(jié)果表明混合格式并未影響接觸間斷的分辨率。此外,SLAU格式與K-CUSP-X格式在此問(wèn)題中也能夠得到準(zhǔn)確結(jié)果。

4.2 亞聲速圓柱擾流

平板邊界層流動(dòng)較為簡(jiǎn)單,為此進(jìn)一步采用圓柱的黏性繞流算例測(cè)試多維黏性流動(dòng)特性的計(jì)算精度[15]。計(jì)算網(wǎng)格如圖7所示,圓柱附近采用200×25個(gè)四邊形網(wǎng)格離散,外部采用13 895個(gè)三角形網(wǎng)格離散。計(jì)算來(lái)流馬赫數(shù)為0.3,雷諾數(shù)為40。

圖7 亞聲速圓柱繞流計(jì)算網(wǎng)格Fig.7 Computation grid of subsonic flow around a cylinder

混合格式的計(jì)算結(jié)果如圖8所示,可見(jiàn)圓柱后形成一對(duì)穩(wěn)定對(duì)稱的分離渦。為進(jìn)一步比較各個(gè)格式,將計(jì)算得到的分離渦長(zhǎng)度列于表1(圓柱直徑為1)。

圖8 混合格式計(jì)算結(jié)果Fig.8 Numerical result of the hybrid scheme

表1 分離渦長(zhǎng)度Table 1 Length of the detached eddies

計(jì)算得到的分離渦長(zhǎng)度同時(shí)受到網(wǎng)格剖分與格式耗散的影響,因此此處并不直接比較相對(duì)實(shí)驗(yàn)的計(jì)算結(jié)果差異。主要的結(jié)論是,在相同的網(wǎng)格條件下,基于TV格式的混合格式計(jì)算結(jié)果與TV格式相差不足0.5%,表明在光滑流場(chǎng)中增加的耗散可以忽略不計(jì)。另外Roe格式、SLAU格式、K-CUSP-X格式耗散也較小,而van Leer格式耗散較大,計(jì)算得到的分離渦最小。

4.3 高超聲速鈍體繞流：一階格式

采用混合格式計(jì)算第二節(jié)中的高超聲速鈍體繞流問(wèn)題,考核其激波穩(wěn)定性。分別采用規(guī)則與不規(guī)則三角形網(wǎng)格,計(jì)算條件不變。計(jì)算結(jié)果如圖9所示。

圖9 混合格式在三角形網(wǎng)格上的密度云圖與流線(一階)Fig.9 Density contour and streamlines of hybrid scheme on triangular grids(first-order)

可以看到,混合格式?jīng)]有出現(xiàn)第二節(jié)中幾種迎風(fēng)格式表現(xiàn)出的激波不穩(wěn)定現(xiàn)象。盡管網(wǎng)格非對(duì)稱,并且不規(guī)則三角形網(wǎng)格中鈍體前緣亞聲速區(qū)的密度分布仍與實(shí)際物理解存在差異,但是混合格式得到了基本對(duì)稱的流場(chǎng)分布,觀察流線也可以看到?jīng)]有出現(xiàn)第二節(jié)中的非物理分離渦。因此可以認(rèn)為,混合格式不僅保持了TV格式黏性流動(dòng)計(jì)算的精度,也改善了強(qiáng)激波計(jì)算的穩(wěn)定性。

4.4 高超聲速鈍體繞流：二階格式

空間二階精度條件下,式(4)中的單元表面兩側(cè)變量差值減小,從而使迎風(fēng)格式耗散減小,提高了流場(chǎng)模擬的精度。另一方面,盡管限制器的作用避免了流場(chǎng)變量重構(gòu)的振蕩,但是由于耗散減小,激波計(jì)算的不穩(wěn)定現(xiàn)象更容易產(chǎn)生。因此本節(jié)在上一小節(jié)算例的基礎(chǔ)上,采用二階格式結(jié)合MLP限制器進(jìn)行計(jì)算,檢驗(yàn)混合格式的性能。

本節(jié)算例的計(jì)算條件和網(wǎng)格與上一小節(jié)相同。鑒于對(duì)高超聲速流場(chǎng)進(jìn)行計(jì)算,激波前后物理量變化劇烈,理論最大壓力比超過(guò)70倍,因此MLP限制器中的Venkatakrishnan限制函數(shù)的參數(shù)K設(shè)為0.1。計(jì)算結(jié)果如圖10所示。

相對(duì)于一階格式的計(jì)算結(jié)果,二階格式中流動(dòng)的對(duì)稱性下降,表明隨著耗散降低,對(duì)擾動(dòng)的抑制能力下降。但是,計(jì)算中并未出現(xiàn)明顯的激波不穩(wěn)定現(xiàn)象,特別是弓形激波的位置與形態(tài)都較為準(zhǔn)確。作為對(duì)比,給出第2節(jié)中穩(wěn)定性最好的K-CUSP-X格式的計(jì)算結(jié)果,如圖11所示,可以看到在規(guī)則三角形網(wǎng)格上波后存在較為明顯的波動(dòng),而在不規(guī)則三角形網(wǎng)格上計(jì)算得到的激波形態(tài)存在異常,波后亞聲速區(qū)存在非定常、非物理的渦結(jié)構(gòu)。上述結(jié)果表明,在二階精度條件下采用混合格式能夠有效改善計(jì)算穩(wěn)定性。

圖10 混合格式在三角形網(wǎng)格上的密度云圖(二階)Fig.10 Density contour of hybrid scheme on triangular grids(second-order)

圖11 K-CUSP-X格式在三角形網(wǎng)格上的密度云圖(二階)Fig.11 Density contour of K-CUSP-X scheme on triangular grids(second-order)

若進(jìn)一步對(duì)比混合格式在一階、二階精度空間離散中的差異,可以發(fā)現(xiàn)采用一階格式得到的流場(chǎng)更為合理。而二階格式在規(guī)則三角形網(wǎng)格上出現(xiàn)了類似于圖3、但更為為微弱的不合理密度分布。這主要是在二階精度條件下耗散更低導(dǎo)致。文獻(xiàn)[39]討論了限制器對(duì)二階格式計(jì)算穩(wěn)定性的影響,并發(fā)現(xiàn)即使采用較為穩(wěn)定的van Leer格式,也可能出現(xiàn)流場(chǎng)的異常波動(dòng)。因此,在二階空間精度格式的計(jì)算中,不僅需要考慮通量計(jì)算的影響,還需要考慮限制器的作用。特別是在非結(jié)構(gòu)網(wǎng)格的條件下,網(wǎng)格的影響因素更為復(fù)雜,因此需要對(duì)限制器和迎風(fēng)型格式的作用進(jìn)行綜合的研究。而在實(shí)際工程應(yīng)用中,通過(guò)生成高質(zhì)量網(wǎng)格或者采用網(wǎng)格自適應(yīng)技術(shù),避免或降低網(wǎng)格條件可能帶來(lái)的不利影響,提高計(jì)算精度與穩(wěn)定性。顯然,隨著網(wǎng)格的加密,二階格式將更快地在光滑流場(chǎng)得到準(zhǔn)確的計(jì)算結(jié)果,而激波附近的異常也會(huì)減弱。

4.5 高超聲速三維鈍錐繞流

本節(jié)對(duì)比混合格式與TV格式在三維高超聲速外形含激波流場(chǎng)計(jì)算中的性能,采用MLP-pw限制器[30]實(shí)現(xiàn)二階空間精度的激波捕捉計(jì)算。數(shù)值模擬的鈍錐外形[40]半錐角為4.7°。來(lái)流馬赫數(shù)為6,迎角0°。計(jì)算網(wǎng)格的詳細(xì)參數(shù)可參考文獻(xiàn)[41],局部網(wǎng)格如圖12所示。

圖12 鈍錐繞流計(jì)算網(wǎng)格Fig.12 Computation grid of hypersonic flow around a blunt cone

對(duì)比圖13中TV格式與混合格式計(jì)算得到的溫度分布,可見(jiàn)混合格式有效地消除了激波不穩(wěn)定現(xiàn)象,從而得到了更為合理的溫度分布。這對(duì)于氣動(dòng)加熱問(wèn)題中的熱流計(jì)算是必要的。但是必須說(shuō)明的是,準(zhǔn)確的熱流計(jì)算還需要考慮黏性通量計(jì)算、流場(chǎng)(溫度)梯度計(jì)算、網(wǎng)格分布乃至湍流等因素,需要進(jìn)行系統(tǒng)的分析。本文結(jié)果僅用于說(shuō)明混合格式對(duì)數(shù)值激波不穩(wěn)定這一影響因素的抑制作用。

圖13 鈍錐繞流溫度云圖Fig.13 Pressure contours of the computation of hypersonic flow around a blunt cone

5 結(jié) 論

本文對(duì)迎風(fēng)型通量格式的激波捕捉能力進(jìn)行了分析對(duì)比,得到的數(shù)值結(jié)果表明多種性能優(yōu)秀的迎風(fēng)格式在非結(jié)構(gòu)網(wǎng)格上表現(xiàn)出數(shù)值激波不穩(wěn)定現(xiàn)象,并且再次驗(yàn)證了激波橫斷方向的不均勻動(dòng)量擾動(dòng)是激波計(jì)算異常的原因。準(zhǔn)確分辨接觸間斷的格式難以耗散網(wǎng)格切向的動(dòng)量擾動(dòng),而此類擾動(dòng)在非結(jié)構(gòu)網(wǎng)格上由于網(wǎng)格與激波交錯(cuò)分布而易于產(chǎn)生。因此本文得出結(jié)論:修正動(dòng)量方程通量中的對(duì)流項(xiàng)是抑制動(dòng)量擾動(dòng)、提高穩(wěn)定性的關(guān)鍵。通過(guò)結(jié)合van Leer格式的馬赫數(shù)分裂函數(shù),本文構(gòu)造了一種TV-van Leer混合格式,引入了van Leer格式的耗散項(xiàng)形式。該格式僅對(duì)TV格式對(duì)流項(xiàng)的動(dòng)量方程分量進(jìn)行修正,而壓力項(xiàng)以及質(zhì)量、能量方程保持TV格式原有形式。數(shù)值計(jì)算結(jié)果表明該混合格式保持了TV格式的黏性流動(dòng)計(jì)算精度,并且改善了在非結(jié)構(gòu)網(wǎng)格上捕捉強(qiáng)激波的穩(wěn)定性,提高了計(jì)算高超聲速流動(dòng)的能力。而對(duì)于工程CFD分析中常用的二階空間精度非結(jié)構(gòu)有限體積法,其數(shù)值穩(wěn)定性與計(jì)算精度的改進(jìn)工作必須綜合考慮限制器與計(jì)算網(wǎng)格的影響,展開(kāi)系統(tǒng)的分析和研究。