郝世熙，丁秋實，魏 桐，劉正先，劉 偉，李孝檢，趙 明，*

(1.天津大學(xué) 機械工程學(xué)院，天津 300350；2.電子科技大學(xué) 計算機科學(xué)與工程學(xué)院，成都 611731)

0 引言

間斷伽遼金（discontinuous Galerkin,DG）有限元方法[1]是一種求解守恒系統(tǒng)的高精度方法，以其高精度、高緊致性、便于大規(guī)模并行計算、適于處理復(fù)雜邊界、易于實現(xiàn)自適應(yīng)等優(yōu)點，成為CFD 高精度流場仿真的重要方法。

在DG 方法框架內(nèi)，學(xué)者們先后提出LDG[2]、BR2[3]和內(nèi)罰（interior penalty,IP）等方法進行Navier-Stokes（N-S）方程黏性項離散。其中，內(nèi)罰方法將黏性項整體作為輔助變量，通過罰函數(shù)進行相鄰單元間的信息交換，具有緊致性好、易實施等優(yōu)勢[4]。在此基礎(chǔ)上，基于DG 方法開展雷諾平均（Reynolds average Navier-Stokes,RANS）湍流仿真成為重要的研究方向。Hartmann、Bassi 等[5-8]將IPDG 方法應(yīng)用到可壓縮N-S 方程模擬中，與標準k-ω湍流模型相結(jié)合，驗證了該方法在RANS 湍流模擬方面的能力。

國內(nèi)關(guān)于DG 方法的研究應(yīng)用已經(jīng)從Euler 方程拓展到可壓縮N-S 方程的湍流模擬[9]。李喜樂、郝海兵等[10]結(jié)合Baldwin-Lomax 零方程湍流模型，驗證了DG 方法在跨聲速流動模擬中的可靠性。Yu、Yan 等[11]結(jié)合Baldwin-Lomax 湍流模型求解了二維流動問題，在激波解析方面表現(xiàn)出較高精度。賀立新、張來平等[12]構(gòu)造了一類DG/FV 重構(gòu)方法并進行了二維算例的湍流模擬。Cheng、Yang 等[13-14]將直接間斷伽遼金有限元方法（direct DG,DDG）應(yīng)用于RANS方程求解。張濤、呂宏強、秦望龍、陳建偉等[15-17]實現(xiàn)了基于S-A 湍流模型的RANS 和DES 仿真并驗證了精度和收斂性，拓寬了DG 方法的應(yīng)用范圍。剪切應(yīng)力輸運（shear stress transport,SST）k-ω模型普適性強，?；雀?，在航空航天領(lǐng)域應(yīng)用廣泛[18]，但其在內(nèi)罰間斷伽遼金（internal penalty discontinuous Galerkin,IPDG）方法內(nèi)的適用性和精度尚待研究。

與有限體積法相似的是，DG 方法也需要求解單元界面的數(shù)值通量?；诓煌慕扑悸罚嗬^出現(xiàn)了Lax-Friedrichs、van Leer 等[19]通量矢量分裂（flux vector split,FVS）格式，同時誕生了HLL、HLLC、Roe格式[20-21]為代表的通量差分（flux difference split,FDS）格式和綜合兩類格式特點的AUSM 類格式[22]。在基于有限體積法的RANS 仿真中，不同通量格式表現(xiàn)出迥異的數(shù)值特性。例如，HLL 格式計算穩(wěn)定，但是耗散較大會抹平間斷。Roe 和AUSM 等格式具有低耗散特性，但會出現(xiàn)激波不穩(wěn)定等問題[23-24]。在基于IPDG 方法的RANS 仿真中，通常會采用Lax-F、HLLC、Roe 等格式進行對流通量計算[6-17,25-28]，尚缺乏不同工況下格式特性的針對性研究。因此，應(yīng)用幾種代表性的對流通量格式開展IPDG 框架內(nèi)的RANS仿真并結(jié)合不同格式的構(gòu)造原理分析其數(shù)值表現(xiàn)具有一定的理論意義。

為此，本文開展了兩方面工作。首先，在IPDG高精度方法框架內(nèi)基于兩方程SSTk-ω模型的湍流模擬，在亞/跨/超聲速工況下驗證SSTk-ω模型在IPDG 湍流模擬中的適用性和準確性。在此基礎(chǔ)上，對比分析了Lax-F、HLL、Roe、AUSM 通量格式在上述湍流模擬框架中的表現(xiàn)。研究發(fā)現(xiàn)AUSM 格式在超聲速工況下激波面“褶皺”現(xiàn)象明顯，Lax-F 格式和HLL 格式魯棒性較好，但數(shù)值精度和激波解析能力略差，Roe 格式具有寬速域內(nèi)較高的數(shù)值精度和激波解析能力。

1 控制方程

1.1 無量綱雷諾平均Navier-Stokes（RANS）方程

本文基于笛卡爾坐標系下的RANS 方程組開展研究：

式中U代表守恒變量 [ρ,ρu,ρv,ρw,ρE]T，對流通量具體表示為：

黏性通量分量表示為：

式中熱流率矢量為：

剪切應(yīng)力張量為：

式中：黏性系數(shù) μ=μL+μt，層流黏性系數(shù) μL通過蘇士南（Sutherland）公式計算得到，湍流黏性系數(shù)μt通過 SSTk-ω湍流模型方程計算。

1.2 無量綱SST k-ω 湍流模型方程

SSTk-ω湍流模型方程[29]綜合了k-ε和k-ω湍流模型的特點，具有較好的近壁模擬能力和分離預(yù)測能力，較S-A 湍流模型適用于逆壓梯度流動的模擬[30]。

笛卡爾坐標系下的無量綱化SSTk-ω湍流模型：

式中，σk、σω、β、γ 等參數(shù) Φ通過式（8）計算得到；相應(yīng)的 σk1、σω1、β1、γ1等模式參數(shù) Φ1和 σk2、σω2、β2、γ2等模式參數(shù) Φ2，以及 β?、F1、F2等參數(shù)的具體計算見文獻[29]。

渦黏系數(shù)計算公式為：

式中a1、?的具體計算見文獻[29]。

1.3 松耦合求解方法

本文采取松耦合數(shù)據(jù)傳遞實現(xiàn)湍流模擬方程的求解，即獨立求解RANS 方程組和SSTk-ω方程組。在上述框架內(nèi)，RANS 方程組和湍流模型方程組求解采取不同精度：RANS 方程組采用高階DG 方法求解，湍流方程組采用有限體積法求解。高階精度的流場信息通過加權(quán)投影到湍流模型中。上述方法可以有效避免矩陣剛性帶來的收斂性問題，同時規(guī)避湍流模型高精度求解的穩(wěn)定性問題。

2 空間離散方法

2.1 DG 方法

記計算域 ?的邊界為 Γ，劃分單元 ?k，n為單元邊界外法矢量。方程（1）通過與測試函數(shù) φh相乘后在計算域內(nèi)積分得到其弱形式：

式中：Uh表示理論解在有限元空間中的近似；為多項式階數(shù)，uj為單元解的自由度，φj為基函數(shù)，在DG 方法內(nèi)測試函數(shù) φh取為基函數(shù) φj。通過提高單元內(nèi)多項式階數(shù)可以實現(xiàn)DG方法框架內(nèi)的高精度數(shù)值仿真。

2.2 內(nèi)罰黏性項離散方法

由于方程含有變量梯度，采用改進的內(nèi)罰函數(shù)方法進行降階處理，構(gòu)建IPDG 方法[25-26]。將黏性項作為中間變量：

約定 [·] 表示兩側(cè)值的階躍，{·}表示兩側(cè)值的平均，定義提升算子為：

式（16）第二項中單元邊界上黏性通量取為：

懲罰因子Cip建議取值范圍為4 至20[5]。通過數(shù)值算例試驗表明，Cip小于4 時，計算難以收斂，Cip大于20 時計算結(jié)果誤差相對較大，在4 到20 之間其取值對結(jié)果影響不明顯，故本文在保證收斂性的基礎(chǔ)上將懲罰因子取為4。p為多項式階數(shù)，hf表示面的尺度。

對IP 格式空間半離散的式（16）采取BDF2 格式進行時間離散，并采用Newton -GMRES 隱式迭代方法求解，其完整離散表達形式為：

2.3 對流通量離散格式

IPDG 方法需要求解單元界面對流項的數(shù)值通量。本文采用如下4 種格式進行計算：

2.3.1 Lax-F 格式

Lax-F 格式因構(gòu)造簡單而在DG 方法中被廣泛采用，其構(gòu)造表達式為：

式中，λ=max[(V·n+c)L,(V·n+c)R]，V代表當?shù)貑卧乃俣仁噶?，c代表當?shù)芈曀?，下標L 和R 分別表示單元界面兩側(cè)，F(xiàn)代表單側(cè)通量。

2.3.2 HLL 格式

HLL 格式是一種近似黎曼格式，該格式采用雙波假設(shè)，波速分別為sL、sR，計算較為簡單，魯棒性強。

HLL 通量表達式為：

2.3.3 Roe 格式

Roe 格式是一種經(jīng)典的矢通量差分（FDS）格式。通過求解線性近似的Riemann 問題解析間斷結(jié)構(gòu)。Roe 格式表達式為：

2.3.4 AUSM 格式

AUSM 格式綜合了矢通量分裂格式與矢通量差分格式的特點，具有數(shù)值耗散小和計算效率高的優(yōu)點。AUSM 格式將無黏數(shù)值通量分解成由對流速度控制的對流項和由聲速控制的壓力項兩部分，分別進行處理。其中，對流項部分根據(jù)逆變速度采取迎風(fēng)格式，壓力項部分在亞聲速時包含兩側(cè)信息，在超聲速時采取一側(cè)迎風(fēng)信息。

AUSM 格式的具體表達式為：

式中，總焓H=E+p/ρ。逆變馬赫數(shù)和壓力都分解為單元兩側(cè)之和，具體計算式為：

3 算例分析

本課題組在前期工作中基于IPDG 方法的層流仿真精度驗證及數(shù)值測試體現(xiàn)了IPDG 在數(shù)值仿真中的高精度特點[25-26]。本文采取松耦合方式，分別求解RANS 方程及湍流模型方程，通過數(shù)據(jù)傳輸實現(xiàn)流場迭代，基于二階精度IPDG 方法和SSTk-ω湍流模型對NACA 0012翼型進行亞聲速、跨聲速和超聲速工況下的翼型數(shù)值模擬。具體工況參數(shù)見表1。

對半徑20 倍弦長的圓形計算域劃分網(wǎng)格，第一層網(wǎng)格y+=1，法向增長率1.1，邊界層外網(wǎng)格增長率1.6（如圖1 所示）。計算采用遠場壓力邊界條件和絕熱壁面邊界條件。黏性項采用IP 方法處理，對流項分別采用Lax-F、HLL、Roe、AUSM 通量格式進行計算。在跨/超聲速工況下，為保證激波附近計算穩(wěn)定性，采用基于壓力階躍指標的人工黏性方法處理。

3.1 亞聲速工況

亞聲速大迎角工況下，基于SSTk-ω湍流模型的高精度IPDG 流場仿真收斂至定常解。圖2 表明4 種格式在IPDG 湍流模擬中的翼面壓力數(shù)值仿真結(jié)果均和實驗數(shù)據(jù)吻合[32]，充分驗證本文所建立的湍流仿真方法在亞聲速工況下的準確性。

圖2 亞聲速工況4 種通量格式的翼型表面壓力系數(shù)曲線Fig.2 Surface pressure coefficient curves of the airfoil under subsonic condition for four flux schemes

該工況下的翼型實驗參考值升力系數(shù)約為1.45、阻力系數(shù)約為0.015。4 種通量格式計算結(jié)果基本一致，均能準確模擬流場，其中Roe 格式和AUSM 格式精度較高，如表2 所示。當?shù)伛R赫數(shù)分布如圖3 所示，Roe 格式計算的上翼面前緣高速區(qū)馬赫數(shù)略大（圖3 黑框區(qū)域），故上下翼面壓差較大。

圖3 亞聲速工況4 種通量格式的翼型當?shù)伛R赫數(shù)分布圖Fig.3 Local Mach number of the airfoil under subsonic condition for four flux schemes

表2 4 種通量格式計算的升阻力系數(shù)Table 2 Lift/Drag coefficients calculated with four flux schemes

3.2 跨聲速工況

跨聲速工況中翼型兩側(cè)產(chǎn)生激波。本文建立的湍流仿真方法能夠有效分辨激波結(jié)構(gòu)，具備跨聲速流場模擬能力。如圖4 所示[33]，4 種通量格式的激波捕捉能力不同，激波分辨率和位置略有差異。

圖4 跨聲速工況4 種通量格式的上翼面壓力系數(shù)曲線Fig.4 Upper surface pressure coefficient of the airfoil under transonic condition for four flux schemes

Roe 格式和AUSM 格式激波空間位置相對靠后，激波腳后展向速度梯度相對較大，激波分辨率較高。AUSM 格式和Roe 格式均蘊含迎風(fēng)思想。

Roe 格式為避免通量雅可比矩陣特征值較小時違背熵條件所帶來的非物理解，引入熵修正。但該方法帶來了冗余耗散，影響邊界層內(nèi)計算準確性[34]。AUSM 格式對于對流項采用迎風(fēng)思路而壓力項則考慮了兩側(cè)貢獻，較Roe 格式其更為穩(wěn)定，激波分辨率更高（如圖5 所示）。

圖5 跨聲速工況4 種通量格式的翼型當?shù)伛R赫數(shù)分布圖Fig.5 Local Mach number contours of the airfoil under transonic condition for four flux schemes

Lax-F 格式和HLL 格式激波空間位置相對靠前，激波與邊界層相互干擾，激波腳后展向的速度梯度小，激波分辨率相對較低。其中Lax-F 格式的數(shù)值黏性使動能耗散顯著，誘導(dǎo)產(chǎn)生局部分離渦（如圖6所示）。

圖6 Lax-F 格式翼型中段激波腳后動能(Ek)分布云圖及流線圖Fig.6 Kinetic energy (Ek) contour and streamlines behind the shock foot in the middle section of the airfoil computed by the Lax-F scheme

3.3 超聲速工況

超聲速工況下基于4 種通量格式計算的翼型表面壓力系數(shù)分布曲線如圖7 所示，4 種格式的壓力系數(shù)均與參考值吻合[35]，進一步證明本文建立的湍流仿真方法在超聲速流動中的適用性和準確性。

圖7 超聲速工況4 種通量格式的翼型表面壓力系數(shù)曲線Fig.7 Surface pressure coefficient curves of the airfoil under supersonic condition for four flux schemes

圖8 是超聲速流動馬赫數(shù)等值線圖，4 種格式均能準確模擬翼型頭部脫體激波。Lax-F 格式、HLL 格式數(shù)值耗散大，激波面相對較光滑；AUSM 格式數(shù)值耗散小，激波面表現(xiàn)出“褶皺”的Carbuncle 現(xiàn)象，該現(xiàn)象通常認為是由于激波面平行方向缺少足夠的數(shù)值耗散來抑制數(shù)值擾動[34]；Roe 格式激波面最為光滑。

圖8 超聲速工況4 種通量格式的翼型當?shù)伛R赫數(shù)等值線圖Fig.8 Local Mach number iso-lines of the airfoil under supersonic condition for four flux schemes

4 結(jié)論

通過松耦合-數(shù)據(jù)傳遞，在高精度IPDG 框架內(nèi)建立了基于SSTk-ω湍流模型的湍流仿真方法，魯棒性強。在此基礎(chǔ)上進行亞/跨/超聲速工況下的NACA 0012翼型繞流湍流模擬，結(jié)果表明SSTk-ω湍流模型在IPDG 方法中具有寬速域適用性。

比較了AUSM、Lax-F、HLL、Roe 4 種通量格式在IPDG 湍流模擬中的特性。結(jié)果表明：在低速工況下4 種格式的數(shù)值結(jié)果均與實驗數(shù)據(jù)吻合，低耗散的Roe 格式和AUSM 格式精度略高。跨聲速工況下Lax-F 格式和HLL 格式耗散大，激波解析能力較差，Lax-F 格式激波腳后產(chǎn)生局部流動分離。Roe 格式和AUSM 格式的數(shù)值耗散小，激波分辨率較高。超聲速工況下AUSM 格式誘導(dǎo)產(chǎn)生明顯的Carbuncle 現(xiàn)象。

綜合來看，在IPDG 框架內(nèi)基于SSTk-ω模型的RANS 湍流模擬中，Lax-F 格式和HLL 格式數(shù)值耗散大，計算穩(wěn)定易收斂，但激波解析精度不高；AUSM格式在亞聲速和跨聲速均體現(xiàn)較高精度，但在超聲速流動中激波面“褶皺”較明顯；Roe 格式適用于本文涉及到的寬速域湍流模擬。本文采用原始AUSM格式，其修正發(fā)展的AUSM+類格式未在討論范圍，未來將進一步探究AUSM+類格式在IPDG 湍流模擬中的表現(xiàn)，并發(fā)展基于SSTk-ω湍流模型的轉(zhuǎn)捩模式。