高精度IPDG 湍流模擬及通量格式數(shù)值特性

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

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

0 引言

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

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

國內(nèi)關(guān)于DG 方法的研究應(yīng)用已經(jīng)從Euler 方程拓展到可壓縮N-S 方程的湍流模擬[9]。李喜樂、郝海兵等[10]結(jié)合Baldwin-Lomax 零方程湍流模型，驗(yàn)證了DG 方法在跨聲速流動(dòng)模擬中的可靠性。Yu、Yan 等[11]結(jié)合Baldwin-Lomax 湍流模型求解了二維流動(dòng)問題，在激波解析方面表現(xiàn)出較高精度。賀立新、張來平等[12]構(gòu)造了一類DG/FV 重構(gòu)方法并進(jìn)行了二維算例的湍流模擬。Cheng、Yang 等[13-14]將直接間斷伽遼金有限元方法（direct DG,DDG）應(yīng)用于RANS方程求解。張濤、呂宏強(qiáng)、秦望龍、陳建偉等[15-17]實(shí)現(xiàn)了基于S-A 湍流模型的RANS 和DES 仿真并驗(yàn)證了精度和收斂性，拓寬了DG 方法的應(yīng)用范圍。剪切應(yīng)力輸運(yùn)（shear stress transport,SST）k-ω模型普適性強(qiáng)，模化精度高，在航空航天領(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）格式，同時(shí)誕生了HLL、HLLC、Roe格式[20-21]為代表的通量差分（flux difference split,FDS）格式和綜合兩類格式特點(diǎn)的AUSM 類格式[22]。在基于有限體積法的RANS 仿真中，不同通量格式表現(xiàn)出迥異的數(shù)值特性。例如，HLL 格式計(jì)算穩(wěn)定，但是耗散較大會(huì)抹平間斷。Roe 和AUSM 等格式具有低耗散特性，但會(huì)出現(xiàn)激波不穩(wěn)定等問題[23-24]。在基于IPDG 方法的RANS 仿真中，通常會(huì)采用Lax-F、HLLC、Roe 等格式進(jìn)行對(duì)流通量計(jì)算[6-17,25-28]，尚缺乏不同工況下格式特性的針對(duì)性研究。因此，應(yīng)用幾種代表性的對(duì)流通量格式開展IPDG 框架內(nèi)的RANS仿真并結(jié)合不同格式的構(gòu)造原理分析其數(shù)值表現(xiàn)具有一定的理論意義。

為此，本文開展了兩方面工作。首先，在IPDG高精度方法框架內(nèi)基于兩方程SSTk-ω模型的湍流模擬，在亞/跨/超聲速工況下驗(yàn)證SSTk-ω模型在IPDG 湍流模擬中的適用性和準(zhǔn)確性。在此基礎(chǔ)上，對(duì)比分析了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）方程

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

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

黏性通量分量表示為：

式中熱流率矢量為：

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

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

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

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

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

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

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

式中a1、?的具體計(jì)算見文獻(xiàn)[29]。

1.3 松耦合求解方法

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

2 空間離散方法

2.1 DG 方法

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

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

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

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

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

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

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

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

2.3 對(duì)流通量離散格式

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

2.3.1 Lax-F 格式

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

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

2.3.2 HLL 格式

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

HLL 通量表達(dá)式為：

2.3.3 Roe 格式

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

2.3.4 AUSM 格式

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

AUSM 格式的具體表達(dá)式為：

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

3 算例分析

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

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

3.1 亞聲速工況

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

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

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

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

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

3.2 跨聲速工況

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

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

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

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

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

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

圖6 Lax-F 格式翼型中段激波腳后動(dòng)能(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 種通量格式計(jì)算的翼型表面壓力系數(shù)分布曲線如圖7 所示，4 種格式的壓力系數(shù)均與參考值吻合[35]，進(jìn)一步證明本文建立的湍流仿真方法在超聲速流動(dòng)中的適用性和準(zhǔn)確性。

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

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

圖8 超聲速工況4 種通量格式的翼型當(dāng)?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-ω湍流模型的湍流仿真方法，魯棒性強(qiáng)。在此基礎(chǔ)上進(jìn)行亞/跨/超聲速工況下的NACA 0012翼型繞流湍流模擬，結(jié)果表明SSTk-ω湍流模型在IPDG 方法中具有寬速域適用性。

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

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