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

(大連理工大學 工業(yè)裝備結構分析國家重點實驗室,大連 116024)

0 引 言

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

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

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

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

1 控制方程與基本數值方法

1.1 控制方程及其離散

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

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

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

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

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

完成空間離散求解后,可以采用顯式或隱式方法對式(2)中的時間導數項進行離散求解。本文采用的方法包括顯式一階前差與隱式LU-SGS方法[22]。

1.2 通量計算

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

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

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

對于SLAU格式,其分裂的對流項與壓力項分別為:

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

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

1.3 空間變量重構與限制器

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

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

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

圖1 測試網格示例Fig.1 Schematics of computational grids

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

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

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

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

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

3 混合格式

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

以TV格式為例,借鑒對流迎風分裂格式的構造思路,將其寫為式(4)的形式后,得到以下各項:

其中a為聲速?？梢园l(fā)現(xiàn)TV格式在進行對流項與壓力項的分裂后,其壓力項主要由壓力而不是速度決定,而壓力場的橢圓型方程性質決定它的作用形式與對流項的雙曲型方程擾動形式不同,因而認為壓力項對激波穩(wěn)定性的影響將是次要的。基于這一考慮,下面的混合格式將不考慮壓力項的修正。

為了簡化處理,將TV格式在式(4)中的對流項進行合并,可以恢復其在文獻[24]中的形式:

由此可見,TV格式中對流項計算的關鍵是得到合理的速度近似,準確的速度近似使其獲得了良好的接觸間斷分辨率[37],同時也導致了激波計算穩(wěn)定性的下降[33]。

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

其中Ma±為單元表面左右兩側法向馬赫數。需要強調的是,不僅壓力項的處理保持TV格式原有形式,上式中的加權混合僅對動量方程進行,而質量方程與能量方程的對流項仍保持TV格式原有形式。最后定義上式中權系數如下[20]:

式(11)中的權系數通過計算待求通量所在的面元左右兩側單元的所有表面上的壓力差得到,如圖4所示,k(i)和k(j)分別為單元i和j的所有表面。從而確?？缭郊げǖ木W格面能夠隨著壓力差的增大切換至耗散格式。反之,在壓力差較小的區(qū)域,如均勻流場或接觸間斷附近,混合格式完全保持TV格式的原有形式,從而保持其低耗散、高分辨率的性質。

圖4 混合格式權系數定義示意Fig.4 Schematic of the definition of the weights in the hybrid scheme

4 算例與分析

4.1 亞聲速層流邊界層

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

圖5 層流平板邊界層流動計算網格Fig.5 Computation grid of flat plate boundary layer flow

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

圖6 層流平板邊界層結果對比Fig.6 Comparison of the results of flat plate boundary layer

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

4.2 亞聲速圓柱擾流

平板邊界層流動較為簡單,為此進一步采用圓柱的黏性繞流算例測試多維黏性流動特性的計算精度[15]。計算網格如圖7所示,圓柱附近采用200×25個四邊形網格離散,外部采用13 895個三角形網格離散。計算來流馬赫數為0.3,雷諾數為40。

圖7 亞聲速圓柱繞流計算網格Fig.7 Computation grid of subsonic flow around a cylinder

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

圖8 混合格式計算結果Fig.8 Numerical result of the hybrid scheme

表1 分離渦長度Table 1 Length of the detached eddies

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

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

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

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

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

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

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

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

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

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

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

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

4.5 高超聲速三維鈍錐繞流

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

圖12 鈍錐繞流計算網格Fig.12 Computation grid of hypersonic flow around a blunt cone

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

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

5 結 論

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