HLLE++格式在高馬赫數空腔流動模擬中的應用

張培紅， 羅 磊， 賈洪印， 趙 煒， 張耀冰， 吳曉軍

(中國空氣動力研究與發展中心 計算空氣動力研究所，綿陽 621000)

1 引 言

空腔流動現象廣泛存在于內埋武器艙、起落架艙和表面縫隙等航空飛行器工程實際中，由于其巨大的工程實際意義和流動的復雜性，長期以來受到廣泛關注和研究[1]。特別是數值模擬方法的發展，為空腔流動的研究提供了新的強有力的工具和手段。相關的研究表明[2-5]，采用RANS方程進行數值模擬得到的結果，在空腔靜態流動形態、底板中心線平均壓力分布系數等方面與試驗結果吻合較好。如早期采用MacCormack顯式二步預估差分格式計算亞聲速和超聲速空腔流動[3]，以及采用Brailovskays差分算法[4]、三階迎風TVD格式和隱式格式[5]等，都取得了很好的效果。

數值格式是數值求解的基礎和靈魂，數值格式的發展，大大促進了數值模擬在空腔流動研究中的應用。Jameson等[6]提出了著名的JST格式，通過構造人工粘性很好地解決了激波附近的振蕩問題。隨后，迎風格式的研究極大推進了數值求解技術的發展，Godunov格式的提出[7]為迎風格式的發展提供了新的方向。特別是針對Godunov格式中求解Riemann問題計算量大的缺陷，Roe[8]提出了對激波和接觸間斷都有很好分辨率的近似Riemann求解器，更是把數值模擬的應用推向了巨大成功。Hamed等[9]采用三階Roe格式，研究了馬赫數1.19的空腔流動，得到了與實驗數據較為一致的結果。Arunajatesan等[10]采用數值模擬方法研究了空腔不同長寬比的影響。

侯中喜等[11]采用高精度和高分辨率格式開展了超聲速開式空腔流動特性研究。馬明生等[12]通過求解可壓縮流動，研究了不同L/D對空腔流動特性的影響，分析了空腔流動類型隨L/D增大轉變的機理。李曉東等[13]通過求解N-S方程，對亞音速空腔流動進行研究。肖虹等[14]利用求解雷諾平均N-S方程研究了飛行器腹部空腔繞流情況。文獻[15，16]對亞、跨和超聲速條件下空腔流場特性進行了較為深入系統的研究。

上述研究多是針對馬赫數2以下的空腔流動開展研究，而對于馬赫數大于2的高馬赫數空腔流動研究較少。高馬赫數空腔流動的強激波和強剪切的流動特點，對數值格式的數值精度和穩定性提出了更高的要求。HLLE++格式兼具了Roe格式的數值精度與HLLE+格式的穩定性，具有可在2個網格內捕捉到激波、無熵增、邊界層耗散小、精確捕捉接觸間斷和可大大降低紅玉現象的發生等優點。

本文結合非結構網格和HLLE++格式的特點，通過改進激波探測的求解，把HLLE++格式應用到高馬赫數空腔流動的數值模擬中，有效解決了高馬赫數空腔流動中激波捕捉和流動分離模擬對數值精度和穩定性的要求，取得了較好的效果，為開展高馬赫數空腔流動研究提供了可靠手段。

2 HLLE++格式及其改進

HLL族格式是由Harten等[17]提出的一族Godunov格式，HLL格式用內部狀態的平均值來簡化單元界面處的Riemann問題解。最簡單的HLL格式假設在最大左行波和右行波之間只有一個內部狀態，如圖1所示。Einfeldt[18]基于最簡單的HLL格式假設，提出了HLLE格式。HLLE格式具有很好的魯棒性，但存在擴散性過大和接觸間斷信息丟失的缺點。

圖1 HLL格式假設

針對HLLE格式的缺點，Einfeldt[18]在HLLE格式的基礎上增加了針對線性特征場的反擴散項，發展了HLLEM格式。Park等[19]注意到利用 HLLEM 格式的形式，選擇不同的波速和接觸間斷速度，可以分別得到HLLE，HLLEM和Roe格式，故提出結合HLLEM格式的波速判斷和Roe格式的接觸間斷速度判斷，形成了HLLE+格式。Tramel等[20]在HLLE+的基礎上，采用新的開關函數，同時將特征值替換為HLLE+的特征值和Roe的耗散形式的特征值的混合組合形式，形成HLLE++格式，可表示為

(1)

(2,3)

(4)

(5)

式中β為開關函數，

βHLLE + +=MAX(βHLLE + +, 0.4)

βHLLE + +=MIN(βNEW,C 1,βNEW,C 2)

(6)

Swi,j,k=MAXl =i -1, i +1m =j -1, j +1n =k -1, k +1(kpl,m,n)

(7)

非結構網格不像結構網格那樣具有明確的i,j和k三個方向[21，22]，式(7)的激波探測公式不再適用。本文針對非結構網格的特點，對激波探測的求解方式進行了改進，目的是最大可能地探測網格點i周圍的激波單元，提高格式的數值精度和穩定性，可表示為

(8)

式中j為與單元i共點的所有單元，包括i自身。

(9)

式中j和k為與單元i共面的單元。

δ為輸入參數，取值范圍為0～1，取值越小，格式耗散越大，穩定性越好，一般取為0.1。

3 HLLE++格式驗證

為了驗證HLLE++格式對高馬赫數流動的模擬能力，以超聲速無粘斜坡為算例，采用Roe，HLLC和HLLE++三種格式進行對比計算分析。無粘斜坡的計算條件為M=4.0，H=0 km，斜坡坡度為30°。圖2給出了無粘斜坡計算網格。圖3～圖5分別給出了不同格式計算得到的密度云圖。

圖2 無粘斜坡計算網格

圖3 Roe格式計算結果

圖4 HLLC格式計算結果

圖5 HLLE++格式計算結果

可以看出，對于不是沿激波排列的網格，Roe格式和HLLC格式在斜坡超聲速激波附近出現了明顯的紅玉現象，沿網格分布激波出現明顯的臺階狀，產生的誤差傳播到激波的下游。出現這種現象的原因是Roe和HLLC等近似Riemann解格式都是從一維流動推導而來，當推廣到二維和三維時，會帶來數值誤差源。近似Riemann解格式對網格沿激波排列具有高度的敏感性，當網格沒有沿激波排列時，在不同方向上的Riemann問題就會發生交叉耦合現象。穿過激波的守恒變量的強跳躍受Riemann求解器在非排列網格的各個方向上錯誤的表征為其他特征變量的強跳躍，對于具有數值耗散的格式，這將會帶來格式的數值誤差。由于其直接放大了單元界面的法向速度，如果單元面法向速度方向剛好是激波的切線方向，會導致本來很小的速度分量錯誤放大。而在單元界面的切向方向由于耗散的缺乏，造成交叉耦合擾動或者放大，或者保持中立穩定，即擾動無阻尼地傳播。在這兩種情況下，解的質量嚴重退化，就會產生非物理特性，通常稱為紅玉現象。HLLE++格式大大改善了這種缺陷，在捕捉高超聲速激波時不會發生紅玉現象，同時保持在光滑區域的低數值耗散特性。

為進一步驗證HLLE++在邊界層區域的低耗散特性和在激波捕捉時能夠避免紅玉現象發生的能力，采用本文數值方法模擬了高馬赫數經典標準算例雙橢球的流場和熱流分布。雙橢球模型為兩個具有不同軸長的同心半橢球垂直相貫形成，后段接一等截面橢圓柱構型，模型全長215 mm，橫向最大尺寸131.6 mm，最大高度105.3 mm。計算馬赫數7.8，攻角0°，來流靜溫74.42 K，壁溫296.0 K。計算時網格量為1322184，其中四面體539466，三棱柱777832，金字塔4886。物面網格為267582個三角形，附面層第一層為1.6×10-6m，對應網格雷諾數為45，增長率為1.2。圖6給出了計算得到的雙橢球對稱面流場?？梢钥闯鲎撉暗姆蛛x、再附和二次分離以及頭部弓形激波、激波/激波干擾等流動特征得到很好的模擬，在頭部激波附近未出現紅玉現象。圖7給出了計算得到的雙橢球對稱面中心線熱流分布與試驗比較。熱流計算對數值方法在邊界層區域的低耗散特性提出了更高要求，可以看出，計算值與試驗值吻合較好，說明本文數值方法在邊界層區域保持了低耗散特性。

圖6 雙橢球對稱面流場

圖7 雙橢球對稱面中心線熱流分布

4 高馬赫數空腔流動數值模擬

4.1 網格影響研究

為驗證HLLE++格式模擬高馬赫數空腔流動的可靠性和有效性，以文獻[23]的空腔試驗模型為算例，開展了數值模擬和網格影響研究。

該試驗模型的空腔后壁板和側壁板均可活動，用以調節空腔的長深比和長寬比，空腔前后壁板和底板上沿對稱面布置了多個測壓孔，如圖8所示。本次計算采用的是長深比為24的空腔模型，空腔模型的主要幾何參數和來流參數列入表1。

表1 空腔模型參數和來流參數

圖8 空腔試驗模型

計算網格是影響數值模擬精度的主要因素之一，只有網格分布合理和網格量足夠時，模擬結果才會受網格影響較小或者與網格無關，因此數值模擬中通常會開展網格的無關性或網格的收斂性研究。本文采用4套網格由疏到密，網格量分別約為400萬、1000萬、2400萬和4200萬，并按此順序分別命名為Grid1,Grid2,Grid3和Grid4。網格劃分采用非結構的四面體+三棱柱+六面體單元的混合網格策略，在空腔內全部使用六面體(矩形)單元以保證較好的正交性，空腔外為混合網格單元。圖9 給出了Grid3的對稱面網格。計算時，無粘通量采用HLLE++格式離散，湍流模型采用k-ωSST兩方程模型。

圖9 空腔對稱面網格

圖10給出了4套網格的殘值收斂曲線?？梢钥闯?，四套網格計算收斂都很好，殘差下降5個量級以上。圖11給出了不同網格計算得到的對稱面上空腔底部、前緣平板和后緣平板上的壓力系數分布與試驗結果比較。圖12給出了不同網格計算得到的對稱面上空腔前壁板和后壁板上的壓力系數分布與試驗結果比較。其中，FP,RP,FL,FF和RF分別指空腔前緣平板、空腔后緣平板、空腔底板、空腔前壁板和空腔后壁板上的試驗測量結果。

圖10 不同網格殘差收斂曲線

從圖11和圖12可以看出，對于上述4套網格，空腔上游流動、前壁板分離流動和兩激波間附著平行流動處的模擬結果與試驗符合較好。數值計算結果與試驗結果的差異主要體現在后壁板分離流動和空腔下游流動的預測上。從圖11可以看出，數值計算預測的后壁板分離流動在底板上的分離點位置相比于試驗靠前，且網格越粗預測的分離點位置越靠前，與試驗值偏離越遠；隨網格不斷加密，分離點位置逐漸向后移動，靠向試驗值；可以看出，對于較密的兩套網格Grid3和Grid4，其計算結果幾乎沒有差別，但分離點位置與試驗值相比仍有差異，這也導致了后壁板分離渦處空腔底板以及空腔后緣平板上的壓力系數分布與試驗存在一定差異。圖12的結果同樣表明，空腔后壁板上壓力分布的預測結果隨著網格加密逐漸接近試驗結果，但網格加密到一定程度后，計算結果不再變化，且與試驗結果存在一定差異。上述不同疏密網格的計算結果表明，本文采用的數值計算方法具有網格收斂性，空間離散誤差隨著網格加密迅速降低。同時，說明本文建立的HLLE++方法可以較好地模擬高馬赫數空腔流動。

圖11 不同網格對空腔對稱面壓力分布影響

圖12 不同網格對空腔前后壁壓力分布影響

4.2 湍流模型影響研究

對于激波、剪切層和大分離渦之間具有顯著相互作用的高馬赫數空腔流動而言，湍流模型會對數值模擬結果產生較為明顯的影響?；?.1節的空腔算例，考察了工程中常用的S -A一方程湍流模型和Menter的k-ωSST兩方程湍流模型對模擬結果的影響。計算時網格采用Grid3網格。

圖13和圖14給出了S -A一方程模型和k-ωSST兩方程模型的計算結果對比?？梢钥闯?，兩種湍流模型對于空腔后臺階(前壁板)處分離渦尺度的模擬結果和對空腔內兩分離渦之間的壓力平臺模擬結果基本一致，均與試驗值符合較好。對于空腔前臺階(后壁板)處的分離流動，兩種模型對分離點位置的預測差異不大；分離點后兩種湍流模型計算結果，特別是后壁板附近底部壓力分布差異較大，SST模型的計算結果與試驗值更接近，而S -A模型的計算結果明顯高于試驗值。這說明兩種湍流模型對空腔前臺階處分離區流動的計算結果存在明顯差異。

圖13 不同湍流模型對空腔對稱面壓力分布影響

圖14 不同湍流模型對空腔前后壁壓力分布影響

圖15～圖17分別給出了兩種湍流模型計算得到的空腔對稱面壓力系數云圖、馬赫數云圖和總壓恢復系數云圖。其中，總壓恢復系數定義為當地總壓與來流總壓之比。通過對比發現S -A模型計算的空腔前臺階(后壁板)處的分離渦尺度要明顯大于SST模型的計算結果，主要體現在方向z上分離渦高度，S -A模型計算結果明顯大于SST模型計算結果。同時可以注意到，S -A模型計算的剪切層厚度要明顯大于SST模型，這對空腔外流動的影響表現為激波后總壓恢復系數存在一定差異。由此可以看出，湍流模型主要影響的是剪切流動的計算結果，具體體現在剪切層/邊界層厚度、分離區尺度和剪切輸運效果的差異上，這些差異又同時導致了空腔外的流動狀態，如激波強度和壓縮/膨脹波傳播角度等存在差異。結合圖13和圖14壓力系數對比結果，可以看出，k-ωSST兩方程湍流模型在對高馬赫數流動問題的模擬上，相比于S -A模型更加接近真實流動情況。

圖15 不同湍流模型得到的壓力云圖比較

圖16 不同湍流模型得到的馬赫數云圖比較

圖17 不同湍流模型得到的總壓恢復系數云圖比較

5 結 論

針對非結構混合網格，通過改進激波探測的求解方法，建立了基于非結構混合網格的HLLE++計算方法，并應用于高馬赫數空腔流動的模擬，可以得到以下結論。

(1) HLLE++格式在邊界層區域可以較好地保持低數值耗散特性，在捕捉高超聲速激波時可以有效避免紅玉現象的發生。

(2) HLLE++格式可以較好地模擬高馬赫數空腔流動特性，數值計算方法具有網格收斂性。

(3)k-ωSST兩方程湍流模型在對高馬赫數流動問題的模擬上，相比于S -A模型更加接近真實流動情況。