HLLE++格式在高馬赫數(shù)空腔流動(dòng)模擬中的應(yīng)用

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

(中國空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，綿陽 621000)

1 引 言

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

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

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

上述研究多是針對(duì)馬赫數(shù)2以下的空腔流動(dòng)開展研究，而對(duì)于馬赫數(shù)大于2的高馬赫數(shù)空腔流動(dòng)研究較少。高馬赫數(shù)空腔流動(dòng)的強(qiáng)激波和強(qiáng)剪切的流動(dòng)特點(diǎn)，對(duì)數(shù)值格式的數(shù)值精度和穩(wěn)定性提出了更高的要求。HLLE++格式兼具了Roe格式的數(shù)值精度與HLLE+格式的穩(wěn)定性，具有可在2個(gè)網(wǎng)格內(nèi)捕捉到激波、無熵增、邊界層耗散小、精確捕捉接觸間斷和可大大降低紅玉現(xiàn)象的發(fā)生等優(yōu)點(diǎn)。

本文結(jié)合非結(jié)構(gòu)網(wǎng)格和HLLE++格式的特點(diǎn)，通過改進(jìn)激波探測的求解，把HLLE++格式應(yīng)用到高馬赫數(shù)空腔流動(dòng)的數(shù)值模擬中，有效解決了高馬赫數(shù)空腔流動(dòng)中激波捕捉和流動(dòng)分離模擬對(duì)數(shù)值精度和穩(wěn)定性的要求，取得了較好的效果，為開展高馬赫數(shù)空腔流動(dòng)研究提供了可靠手段。

2 HLLE++格式及其改進(jìn)

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

圖1 HLL格式假設(shè)

針對(duì)HLLE格式的缺點(diǎn)，Einfeldt[18]在HLLE格式的基礎(chǔ)上增加了針對(duì)線性特征場的反擴(kuò)散項(xiàng)，發(fā)展了HLLEM格式。Park等[19]注意到利用 HLLEM 格式的形式，選擇不同的波速和接觸間斷速度，可以分別得到HLLE，HLLEM和Roe格式，故提出結(jié)合HLLEM格式的波速判斷和Roe格式的接觸間斷速度判斷，形成了HLLE+格式。Tramel等[20]在HLLE+的基礎(chǔ)上，采用新的開關(guān)函數(shù)，同時(shí)將特征值替換為HLLE+的特征值和Roe的耗散形式的特征值的混合組合形式，形成HLLE++格式，可表示為

(1)

(2,3)

(4)

(5)

式中β為開關(guān)函數(shù)，

β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)

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

(8)

式中j為與單元i共點(diǎn)的所有單元，包括i自身。

(9)

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

δ為輸入?yún)?shù)，取值范圍為0～1，取值越小，格式耗散越大，穩(wěn)定性越好，一般取為0.1。

3 HLLE++格式驗(yàn)證

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

圖2 無粘斜坡計(jì)算網(wǎng)格

圖3 Roe格式計(jì)算結(jié)果

圖4 HLLC格式計(jì)算結(jié)果

圖5 HLLE++格式計(jì)算結(jié)果

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

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

圖6 雙橢球?qū)ΨQ面流場

圖7 雙橢球?qū)ΨQ面中心線熱流分布

4 高馬赫數(shù)空腔流動(dòng)數(shù)值模擬

4.1 網(wǎng)格影響研究

為驗(yàn)證HLLE++格式模擬高馬赫數(shù)空腔流動(dòng)的可靠性和有效性，以文獻(xiàn)[23]的空腔試驗(yàn)?zāi)Ｐ蜑樗憷_展了數(shù)值模擬和網(wǎng)格影響研究。

該試驗(yàn)?zāi)Ｐ偷目涨缓蟊诎搴蛡?cè)壁板均可活動(dòng)，用以調(diào)節(jié)空腔的長深比和長寬比，空腔前后壁板和底板上沿對(duì)稱面布置了多個(gè)測壓孔，如圖8所示。本次計(jì)算采用的是長深比為24的空腔模型，空腔模型的主要幾何參數(shù)和來流參數(shù)列入表1。

表1 空腔模型參數(shù)和來流參數(shù)

圖8 空腔試驗(yàn)?zāi)Ｐ?/p>

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

圖9 空腔對(duì)稱面網(wǎng)格

圖10給出了4套網(wǎng)格的殘值收斂曲線。可以看出，四套網(wǎng)格計(jì)算收斂都很好，殘差下降5個(gè)量級(jí)以上。圖11給出了不同網(wǎng)格計(jì)算得到的對(duì)稱面上空腔底部、前緣平板和后緣平板上的壓力系數(shù)分布與試驗(yàn)結(jié)果比較。圖12給出了不同網(wǎng)格計(jì)算得到的對(duì)稱面上空腔前壁板和后壁板上的壓力系數(shù)分布與試驗(yàn)結(jié)果比較。其中，F(xiàn)P,RP,FL,FF和RF分別指空腔前緣平板、空腔后緣平板、空腔底板、空腔前壁板和空腔后壁板上的試驗(yàn)測量結(jié)果。

圖10 不同網(wǎng)格殘差收斂曲線

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

圖11 不同網(wǎng)格對(duì)空腔對(duì)稱面壓力分布影響

圖12 不同網(wǎng)格對(duì)空腔前后壁壓力分布影響

4.2 湍流模型影響研究

對(duì)于激波、剪切層和大分離渦之間具有顯著相互作用的高馬赫數(shù)空腔流動(dòng)而言，湍流模型會(huì)對(duì)數(shù)值模擬結(jié)果產(chǎn)生較為明顯的影響。基于4.1節(jié)的空腔算例，考察了工程中常用的S -A一方程湍流模型和Menter的k-ωSST兩方程湍流模型對(duì)模擬結(jié)果的影響。計(jì)算時(shí)網(wǎng)格采用Grid3網(wǎng)格。

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

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

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

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

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

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

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

5 結(jié) 論

針對(duì)非結(jié)構(gòu)混合網(wǎng)格，通過改進(jìn)激波探測的求解方法，建立了基于非結(jié)構(gòu)混合網(wǎng)格的HLLE++計(jì)算方法，并應(yīng)用于高馬赫數(shù)空腔流動(dòng)的模擬，可以得到以下結(jié)論。

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

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

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