HLL－HLLC格式的構(gòu)造與應(yīng)用研究

傅 林，高正紅，張曉東

（西北工業(yè)大學(xué)航空學(xué)院 翼型、葉柵空氣動力學(xué)國防科技重點(diǎn)實驗室，陜西 西安 710072）

0 引 言

Godunov首次將Riemann問題及其相關(guān)解引入到計算流體力學(xué)領(lǐng)域，提出了Godunov方法，將問題轉(zhuǎn)化為求解網(wǎng)格單元界面的Godunov型通量。目前比較成熟和應(yīng)用廣泛的為一階通量數(shù)值計算方法，具有較強(qiáng)的數(shù)值穩(wěn)定性并且可以正確的收斂到熵解。主要有基于通量差分思想和基于矢通量分裂思想的兩類格式，包括著名的ROE，Van Leer等格式。ROE格式在跨聲速和超聲速計算中經(jīng)常由于不滿足熵條件，而要引入熵修正，增加了數(shù)值耗散；Van Leer格式本身耗散較大，適合超聲速流場計算。Harten，Lax，Van Leer等人首先提出了一種求解近似黎曼解的方法－HLL格式［1］。采用HLL格式可以直接求解近似黎曼通量，計算效率較高且穩(wěn)定，適應(yīng)范圍廣。然而由于HLL格式在間斷面具有較大的數(shù)值耗散，計算結(jié)果容易將接觸間斷抹平，在跨聲速復(fù)雜流場計算中誤差較大。Toro等人提出了一種改進(jìn)的方法，稱為HLLC格式［2］。HLLC近似黎曼求解器能夠準(zhǔn)確的捕捉激波、接觸間斷和稀疏波。與精確求解黎曼問題的方法相比具有計算量小、低耗散、高分辨率的特點(diǎn)，且嚴(yán)格滿足熵條件。在低速及跨聲速可壓縮流場中，得到了成功且廣泛的應(yīng)用。

經(jīng)數(shù)值研究發(fā)現(xiàn)在可壓縮流中，當(dāng)馬赫數(shù)較低時，HLLC格式可以獲得比ROE相對于實驗更好的結(jié)果［1］。然而隨著馬赫數(shù)繼續(xù)增加，HLLC格式在強(qiáng)激波區(qū)域會出現(xiàn)激波不穩(wěn)定的現(xiàn)象［3］。相反，HLL格式不會出現(xiàn)激波不穩(wěn)定現(xiàn)象，但是由于HLL格式的耗散過大，得到的解略差于ROE格式。本文研究了HLL，HLLC格式的數(shù)學(xué)原理，構(gòu)建了混合的HLL－HLLC格式。其在低馬赫數(shù)下仍然能獲得較低的耗散，在強(qiáng)激波流場中，可以克服HLLC格式的激波不穩(wěn)定現(xiàn)象。為此，優(yōu)先設(shè)計了激波探測器，用來準(zhǔn)確捕捉流場中的強(qiáng)激波區(qū)域。在強(qiáng)激波區(qū)域附近，可實現(xiàn)耗散較大的具有HLL性質(zhì)的格式，以克服激波不穩(wěn)定現(xiàn)象；而在遠(yuǎn)離強(qiáng)激波區(qū)域采用低耗散的具有HLLC性質(zhì)的格式。從而使HLL－HLLC格式在各種馬赫數(shù)條件下都能得到較好的解。

文中計算采用的求解器為作者自行研究開發(fā)的具有定常、非定常數(shù)值模擬能力的三維有限體積代碼－LMNS3D LMNS3D擁有基于CPU計算的MPI并行技術(shù)和基于GPU計算的CUDA并行技術(shù)，能夠進(jìn)行大規(guī)模并行計算。通量計算可選擇不同格式，如ROE，HLL，HLLC［4－6］，采用 MUSCL插值獲得二階空間離散精度，同時具有四種通量限制器可供選擇。湍流模型為在航空界應(yīng)用廣泛的SA一方程模型和K－Omega SST兩方程模型。粘性項采用二階中心差分離散，擁有薄層近似及全NS計算能力。時間推進(jìn)技術(shù)采用LU－SGS隱式方法，同時采用低速預(yù)處理技術(shù)進(jìn)行低速流場數(shù)值模擬。

1 數(shù)值方法

1.1 控制方程

在物理坐標(biāo)系下NS方程為：

經(jīng)過雅克比轉(zhuǎn)換，NS方程轉(zhuǎn)化為計算系下的守恒形式，為本文求解的模型方程：Q′代表有限體體心值，等代表界面上的無粘數(shù)值通量，由一階單調(diào)近似黎曼求解器給出；R′、S′、T′代表粘性數(shù)值通量。要得到高精度格式，要求提高界面上守恒變量值的插值精度。

1.2 HLL及HLLC近似黎曼求解器

HLLC近似黎曼求解器基于HLL發(fā)展而來，HLLC格式將黎曼問題求解區(qū)域由向左傳最快波系、接觸間斷、向右傳最快波系分成四個部分（如圖1所示）。最左邊和最右邊的波系可能為激波或稀疏波，用符號表示。經(jīng)過這些波系，流場變量發(fā)生變化，特別是激波使流場解變得不連續(xù)。中間的波系為接觸間斷，左右兩邊的流場變量為。求解時假設(shè)在各自的扇形區(qū)域里是常量。Harten等人通過對每一個扇形域中守恒變量做積分平均，提出了一種方法來簡化黎曼問題解的結(jié)構(gòu)。Harten、Lax、Van Leer等人首先提出了一種近似黎曼求解器HLL［1］，忽略了中間接觸間斷的影響，即認(rèn)為接觸間斷兩邊的是相同的。理論研究表明［1］：如果解收斂，那么可以收斂到滿足熵條件的物理解。HLL格式是一種非常有效和魯棒的格式，可以直接應(yīng)用于Euler／NS方程。

圖1 黎曼問題示意圖Fig.1 Riemann problem schematic

其解向量的表達(dá)式為：

然而HLL求解器具有容易抹平接觸間斷的特點(diǎn)。為了克服這一弱點(diǎn)，Toro等人提出了HLLC（C表示contact）格式來精確捕捉接觸間斷［2］。其解的結(jié)構(gòu)為：

HLLC格式需要預(yù)先估計最慢波速、中間波系波速和最快波速。本文采用基于ROE平均的直接波速估計方法，具有很好的魯棒性和精度［8］。

則HLLC格式的通量為：

HLLC格式具有以下優(yōu)點(diǎn)［9］：

（1）精確捕捉激波和接觸間斷。

（2）保證計算過程中的標(biāo)量衡正。

（3）增強(qiáng)的熵條件，不需要熵修正，適合跨聲速和低速計算。

（4）不需要知道通量雅克比矩陣的細(xì)節(jié)。

1.3 HLL－HLLC格式的構(gòu)造

1.3.1 激波不穩(wěn)定現(xiàn)象的產(chǎn)生機(jī)理

Liou［10］分析了不同數(shù)值通量計算格式的耗散特性。在分析中，把每種格式的耗散項展開，用主變量ρ、u和p的差分形式表達(dá)，定義密度差和壓力差的系數(shù)為算子D（ρ）、D（p）。研究不同格式（AUSM ＋，AUSMDV，HLLE和ROE）耗散項中的。研究結(jié)果表明：激波不穩(wěn)定現(xiàn)象出現(xiàn)的根源是在求解質(zhì)量通量的算法中包含了D（p）的影響；當(dāng)一個格式求解質(zhì)量通量算法中，具有D（p）＝0的特性，則是一個激波穩(wěn)定的格式。Pandolfi和D′Ambrosio通過進(jìn)一步研究發(fā)現(xiàn)，當(dāng)一種格式顯式的處理了接觸間斷，就會顯示出明顯的激波不穩(wěn)定現(xiàn)象［11］。從已有的研究中可以看出，HLLC格式之所以在強(qiáng)激波區(qū)域出現(xiàn)激波不穩(wěn)定現(xiàn)象，是因為顯式的表達(dá)了接觸間斷。

具體分析HLLC格式發(fā)現(xiàn)：連續(xù)性方程并不顯式的包含壓力差的貢獻(xiàn)，但是由公式（19）預(yù)估波速的時候，求解中間波系波速的公式中包含壓力差項的影響。

所以HLLC格式無法避免出現(xiàn)激波不穩(wěn)定現(xiàn)象。而同樣的分析可以發(fā)現(xiàn)，由于HLL格式并沒有顯式的表達(dá)中間波系，所以不會出現(xiàn)激波不穩(wěn)定現(xiàn)象。

1.3.2 克服HLLC格式激波不穩(wěn)定現(xiàn)象的方法

本文試圖在HLLC的框架下調(diào)整通過黎曼解中間波系（接觸間斷）的通量差分來消除HLLC的激波不穩(wěn)定現(xiàn)象。減小通過中間波系的解向量的變化，三波系的黎曼問題會退化為兩波系黎曼問題，因此新格式將具有HLL格式的性質(zhì)（如圖2所示）。

圖2 狀態(tài)解構(gòu)造示意圖Fig.2 Construction schematic of state solution

f為本文定義的激波捕捉函數(shù)［12］：

其中pL和pR是作用在網(wǎng)格面上左右兩邊的壓力。在激波以外的區(qū)域，壓力連續(xù)的變化，所以混合函數(shù)f的值為1.0，此時格式具體表現(xiàn)為HLLC的性質(zhì)；在強(qiáng)激波附近，壓力差pR－pL很大，混合函數(shù)f為一個小于1的值，格式表現(xiàn)為HLLC和HLL的混合型；當(dāng)f很小時，格式表現(xiàn)為完全的HLL性質(zhì)，從而抑制激波不穩(wěn)定現(xiàn)象的產(chǎn)生。在三維應(yīng)用中，三個方向的混合函數(shù)f的定義為：

2 HLL－HLLC格式在可壓縮速流動中的應(yīng)用

2.1 高超聲速雙楔流動數(shù)值模擬

Neuenhahm等人為了研究前緣半徑Rn大小和壁面溫度對高超聲速雙楔流動的影響做了一系列實驗［13］。超聲速雙楔流動包含強(qiáng)激波附面層干擾和小分離現(xiàn)象，是考察數(shù)值計算格式和求解器的典型算例。本算例計算模型如圖3。

圖3 高超聲速雙楔流動模型Fig.3 Hypersonic flow on the model of double wedge

計算條件為：來流馬赫數(shù)為8.1，靜壓為520Pa，靜溫為106K，單位長度雷諾數(shù)為3.8×106。壁面為等溫?zé)o滑移邊界條件，Tw為300K。如圖4所示，計算網(wǎng)格為三塊點(diǎn)對接網(wǎng)格，y＋的最大值為0.9。無粘通量分別采用四種數(shù)值格式－ROE（不加熵修正）、HLL、HLLC、HLL－HLLC。計算采用層流模型。

圖4 超聲速雙楔模型計算網(wǎng)格Fig.4 Computational grid of double wedge

圖5為四種格式計算壓力云圖的比較，可以看出HLLC格式在這種強(qiáng)激波流動中出現(xiàn)了激波不穩(wěn)定現(xiàn)象，激波后出現(xiàn)非物理解。ROE、HLL和本文構(gòu)造的HLL－HLLC格式計算結(jié)果較好。圖6為不同格式計算壓力系數(shù)分布與實驗對比，由于實驗洞壁干擾等因素，實驗中取了三次實驗值分別為EXP＿1、EXP＿2、EXP＿3。HLLC格式出現(xiàn)激波不穩(wěn)定現(xiàn)象，壓力分布出現(xiàn)非物理震蕩；同時可以看出，HLL格式耗散過大，分離點(diǎn)靠后與實驗結(jié)果差距較大。本文構(gòu)造的格式HLL－HLLC，計算穩(wěn)定，結(jié)果優(yōu)于HLL格式，接近于ROE格式。從圖7可以看出，文中設(shè)計的激波感受因子能精確捕捉強(qiáng)激波。圖8給出了HLL－HLLC格式計算的雙楔交界處的小分離區(qū)域流線。

圖5 壓力系數(shù)分布云圖Fig.5 Comparison of Cp distribution

圖6 壓力系數(shù)分布與實驗比較Fig.6 Comparison of aerodynamic coefficient

圖7 激波感受因子捕捉到的激波示意圖Fig.7 Shock captured by shock feeling factor

圖8 小分離區(qū)域流線示意圖Fig.8 Streamline schematic of small separation region

2.2 超聲速后臺階流動數(shù)值模擬

超聲速后臺階流動是典型的復(fù)雜流動問題（如圖9），流動現(xiàn)象包括：湍流邊界層、膨脹波、流動分離和再附著、再附激波等，從而成為考察CFD求解器和計算格式的典型算例之一。Smith［14］等人對此問題進(jìn)行了一系列實驗，分別研究了不同臺階高度（0in、0.433in、0.75in），馬赫數(shù)從2.5到5.0，雷諾數(shù)從2.5×105到1.8×106等溫壁面下的流動情況。

圖9 超聲速后臺階流動模型Fig.9 Supersonic flow over the model of backward－facing step

本文計算條件為：來流馬赫數(shù)為2.5，靜壓為1.5×104Pa，靜溫為170K，單位長度雷諾數(shù)為1.8×107。如圖10，臺階高度為0.433in，計算網(wǎng)格為三塊對接網(wǎng)格，y 的最大值為1.25 算例假設(shè)壁面條件為無滑移絕熱壁面，計算區(qū)域的上邊界為滑動壁面。分別采用ROE格式和本文發(fā)展的HLL－HLLC格式計算無粘通量，湍流模型采用SA一方程湍流模式。

圖10 超聲速后臺階流動計算網(wǎng)格Fig.10 Computational grid of backward－facing step

圖10中，x軸和y軸均采用臺階高度做歸一化處理，按照實驗中的方法，壓力采用x／H＝－3.39處的壓力值進(jìn)行歸一化處理。圖11顯示，兩種格式計算的空間壓力分布云圖一致。從圖12馬赫數(shù)云圖可以看出，Prandtl－Meyer膨脹波、再附著激波和分離線非常明顯。圖13表明，本文兩種格式計算的再附點(diǎn)位置均為x／H＝2.5512附近，文獻(xiàn)［15］采用SA湍流模型計算的分離點(diǎn)位置為x／H＝2.53。從圖14可以看出，壁面壓力分布和實驗值吻合較好。其中HLL－HLLC格式和不加熵修正的ROE格式計算結(jié)果一致，均表現(xiàn)出較高的分辨率；在分離區(qū)的壓力分布好于加熵修正的ROE格式。從本算例可以看出，在中等超聲速問題中本文發(fā)展的HLL－HLLC格式，和ROE格式具有相似的耗散性和精度。

圖11 壓力分布云圖比較Fig.11 Comparison of Cp distribution

圖12 馬赫數(shù)云圖比較Fig.12 Comparison of Ma distribution

圖13 流線圖比較Fig.13 Comparison of streamline

圖14 壓力分布與實驗比較Fig.14 Comparison of aerodynamic coefficient

2.3 高超聲速鈍頭體數(shù)值模擬

計算條件為：來流馬赫數(shù)10，單位長度雷諾數(shù)1×105，壁面為等溫壁面274K。模型為半徑1m的半球。分別采用 HLLC，HLL－HLLC、加熵修正的ROE和AUSM＋四種格式計算，計算模式為層流。

從圖15可以看出，HLLC格式計算結(jié)果出現(xiàn)了激波不穩(wěn)定現(xiàn)象，捕捉到的激波為非物理解；另外三種格式HLL－HLLC、加熵修正的ROE和AUSM＋均正確反映了鈍頭體前的弓形激波。其中加熵修正的ROE格式捕獲的激波厚度最厚，HLL－HLLC和AUSM＋格式計算得到的激波厚度基本相當(dāng)；HLLHLLC相對于AUSM＋格式得到的激波后壓力分布更加光滑。另，作者嘗試在本算例中使用未加熵修正的ROE格式，但計算失敗。由此可以看出，本文發(fā)展的HLL－HLLC近似黎曼求解器具有較低的數(shù)值耗散性和較高的精度。

圖15 壓力分布云圖比較Fig.15 Comparison of Cp distribution

3 結(jié) 論

（1）近似黎曼求解器是否出現(xiàn)激波不穩(wěn)定現(xiàn)象強(qiáng)烈依賴于其是否顯式捕捉接觸間斷。HLLC格式精確的顯式表達(dá)接觸間斷，在高超聲速流場中，強(qiáng)激波附近容易出現(xiàn)激波不穩(wěn)定現(xiàn)象；HLL格式由于抹平接觸間斷不會出現(xiàn)這種現(xiàn)象。本文發(fā)展的HLLHLLC格式有效的解決了HLLC格式的激波不穩(wěn)定現(xiàn)象，同時相對于HLL格式具有更低的耗散性和更高的精度。

（2）本文設(shè)計的激波捕捉函數(shù)能夠準(zhǔn)確捕捉強(qiáng)激波。在強(qiáng)激波附近，HLL－HLLC表現(xiàn)為 HLL性質(zhì)的格式，以克服激波不穩(wěn)定現(xiàn)象；而在遠(yuǎn)離強(qiáng)激波區(qū)域過渡到HLLC性質(zhì)的格式。通過三個算例證明了本文所構(gòu)造格式的有效性和魯棒性。

［1］TORO E F.Riemann solvers and numerical methods for fluid dynamics［M］.Springer：Berlin，1999.

［2］TORO E F，SPRUCE M，SPEARES W.Restoration of the contact surface in the HLL－Riemann solver［J］.Shock Waves，1994，4（1）：25－34.

［3］PANDOLFI M，D’AMBROSIO D.Numerical instabilities in upwind methods：analysis and cures for the‘Carbuncle’phenomenon［J］.Journal of Computational Physics，2001，166（2）：271－301.

［4］NIU Y Y，CHANG C H，TSENG W Y，et al.Numerical simulation of an aortic flow based on a HLLC type incompressible flow solver［J］.Communications In Computational Physics，2009，5（1）：142－162.

［5］HARTEN A，LAX P D，LEER B V.On upstream difference and Godunov－type schemes for hyperbolic conservation laws［J］.SIAM Review，1983，25（1）：35－61.

［6］KITAMURA K，ROE P，ISMAIL F.Evaluation of Euler fluxes for hypersonic heating computations［J］.AIAA Journal，2009，47（1）：44－52.

［7］REMAKI L，XIE Z Q，HASSAN O，et al.A high order finite volume－HLLC solver and anisotropic Delaunay mesh adaptation［R］.AIAA 2009－1498，2009.

［8］BATTEN P，CLARKE N，LAMBERT C，et al.On the choice of wave speeds for the HLLC Riemann solver［J］.SIAM Journal on Scientific Computing，1997，18（6）：1553－1570.

［9］NICHOLS R H，TRAMEL R W，BUNING P G.Solver and turbulence model upgrades to OVERFLOW 2 for unsteady and high－speed applications［R］.AIAA－2006－2824，2006.

［10］LIOU M S.Mass flux schemes and connection to shock instability［J］.Journal of Computational Physics，2000，160（2）：623－648.

［11］KIM S S，KIM C G，RHO O H，et al.Cure for the shock instability：development of a shock－stable Roe scheme［J］.Journal of Computational Physics，2003，185（2）：342－374.

［12］QUIRK J J.A contribution to the great Riemann solver debate［J］.International Journal for Numerical Methods in Fluids，1994，18（6）：555－574.

［13］REINARTZ B U，BALLMANN J，BOYCE R.Numerical investigation of wall temperature and entropy layer effects on double wedge shock／boundary layer interactions［R］.AIAA 2006－8137，2006.

［14］SMITH H E.The flow field and heat transfer downstream of rearward facing step in supersonic flow［R］.ARL－67－0056，1967.

［15］DIPPOLD V，MOHLER S R.Validation of the wind－US unstructured flow solver for wall－bounded flows［R］.AIAA 2006－637，2006.