用隱式高精度間斷伽遼金方法模擬可壓層流和湍流

黨亞斌, 劉凱禮, 孫一峰, 楊小權(quán),2,*

(1. 中國商用飛機(jī)有限責(zé)任公司 上海飛機(jī)設(shè)計(jì)研究院, 上海 201210; 2. 中國空氣動力研究與發(fā)展中心 空氣動力學(xué)國家重點(diǎn)實(shí)驗(yàn)室, 四川 綿陽 621052)

0 引 言

近20年來，隨著高精度數(shù)值方法概念的提出和發(fā)展，眾多學(xué)者提出和構(gòu)造了諸多高精度數(shù)值格式。目前高精度方法主要可以分為三大類：有限差分型(Finite Difference, FD)高精度方法，有限體積型(Finite Volume, FV)高精度方法和有限元型(Finite Element, FE)高精度方法。間斷伽遼金方法(Discontinuous Galerkin method, DG)是有限元高精度方法的一種，最先是由Reed和Hill[1]提出，經(jīng)過一些學(xué)者的研究

已經(jīng)能夠用于RANS的計(jì)算[2-5]，但是RANS求解時(shí)的穩(wěn)定性問題成了目前DG方法研究的難點(diǎn)。

DG方法有諸多優(yōu)點(diǎn)[4-6]：1) 能夠保持單元平均值意義下的守恒性，且具有良好的穩(wěn)定性和收斂性；2) 可以簡單通過提高單元內(nèi)分片多項(xiàng)式次數(shù)來提高精度；3) 可以通過使用任意網(wǎng)格來靈活的處理復(fù)雜的物理邊界和幾何區(qū)域；4) 具有緊致性，易于實(shí)施大規(guī)模并行；5) 容易實(shí)現(xiàn)hp自適應(yīng)。與此同時(shí)，DG方法也有一些不足之處，還需要進(jìn)一步解決的問題有如下幾點(diǎn)：1) DG方法使用TVD或TVB限制器抑制數(shù)值振蕩，然而TVD限制器在光滑的極值點(diǎn)附近會降低格式精度，TVB限制器雖然修正了TVD限制器的這一缺陷，卻增加了與問題相關(guān)的參數(shù)，而且對于實(shí)際問題這些參數(shù)只能通過具體問題和相關(guān)經(jīng)驗(yàn)獲得；2) DG方法在如何有效地推廣到NS方程，特別是RANS方程仍然是一個(gè)值得探索的問題，目前，很多學(xué)者依然在探索間斷有限元方法高效處理黏性項(xiàng)的方法和RANS的計(jì)算方法；3) DG方法相對于傳統(tǒng)有限體積方法計(jì)算量較大，對計(jì)算機(jī)計(jì)算效率和存儲空間都提出了更高的要求。

DG方法在處理雙曲守恒律方程和可壓縮無黏Euler方程組是成功而有效的。然而，對于處理擴(kuò)散問題或者含有黏性項(xiàng)的NS方程時(shí)，則帶來了一定的難度。本質(zhì)問題在于，如何構(gòu)造單元界面的間斷問題，即如何構(gòu)造含有間斷的高階黏性數(shù)值通量求解格式。針對這一問題，學(xué)者們提出了各種有效的數(shù)值格式，諸如型罰函數(shù)格式(Symmetric IP, SIP)[7]、局部間斷有限元方法(Local DG, LDG)[8]、(Second Bassi-Rebay scheme, BR2)[3]、BR2格式、LDG的升級版緊致間斷伽遼金格式(Compact DG, CDG)[9]、直接間斷伽遼金格式(Direct DG, DDG)[10-12]。

值得注意的是，很多學(xué)者提出并發(fā)展了一類重構(gòu)型和混合型的間斷有限元方法。這一類方法的基本思想是在DG方法的框架下，借鑒DG方法的優(yōu)勢，通過重構(gòu)方式在原有DG解函數(shù)自由度的基礎(chǔ)上，重構(gòu)高階自由度從而達(dá)到高階精度的目的。Dumbser等人將DG方法和高精度有限體積方法統(tǒng)一到同一PnPm框架下，提出了一類PnPm方法[13]。Luo提出并發(fā)展了一類重構(gòu)型DG方法(reconstruction DG, rDG)[6]，并將這種方法應(yīng)用到求解Euler方程組、NS方程組和RANS方程組，并取得了很好的效果。Van Leer 等人提出并發(fā)展了另一類重構(gòu)DG方法(Recovery DG, RDG)[14]。Zhang等人提出并發(fā)展了一類混合DG/FV方法[15]，通過DG方法和有限體積方法結(jié)合在一起，取得了很好的效果。

采用DG方法求解RANS方程是目前DG方法研究領(lǐng)域的挑戰(zhàn)性問題。雖然國內(nèi)外諸多學(xué)者嘗試將RANS方程組和Spalart-Allmaras (SA)湍流模型耦合在一起[3-4,16]，采用統(tǒng)一的數(shù)值方法進(jìn)行求解，取得了一些研究進(jìn)展，但是采用DG進(jìn)行RANS計(jì)算依然不夠成熟。影響DG進(jìn)行RANS穩(wěn)定求解的因素主要包括以下幾個(gè)方面：(1) 湍流模型的穩(wěn)定性。一方程和兩方程湍流模型在數(shù)值求解中難免會出現(xiàn)求解物理量負(fù)值的情況，這對于高精度算法的穩(wěn)定性影響非常大。因此，Allmaras[17]和Oliver[18]等人分別采用不同的方法對Spalart-Allmaras模型進(jìn)行了修正，使得SA模型的穩(wěn)定性大大提高，能夠用于高精度算法的求解。(2) 隱式時(shí)間格式的穩(wěn)定性。在DG求解過程中，由于顯式格式時(shí)間步長受限，易導(dǎo)致計(jì)算效率低下。因而隱式格式自然變成了一種可以依賴的高效計(jì)算方法。傳統(tǒng)的隱式格式一般采用LU-SGS求解。近年來，有些學(xué)者采用GMRES求解線化后的系統(tǒng)方程組，并利用LU-SGS進(jìn)行預(yù)處理，這樣做的好處是使線性方程組的Jacobian矩陣剛性有所降低，有利于收斂。但是對于黏性計(jì)算，如果采用不精確的無黏通量Jacobian矩陣進(jìn)行預(yù)處理，由于無黏Jacobian矩陣和黏性控制方程數(shù)值格式所對應(yīng)的Jacobian矩陣相差太遠(yuǎn)，必然導(dǎo)致隱式格式的穩(wěn)定性下降，計(jì)算CFL數(shù)降低。在隱式格式進(jìn)行計(jì)算時(shí)，Jacobian矩陣的處理方式有三種：1) 無黏Jacobian矩陣，不考慮對流通量數(shù)值格式的數(shù)值耗散、邊界通量和黏性通量的影響；2) 通過解析數(shù)值格式獲得近似精確的無黏Jacobian矩陣，通常用黏性譜半徑代替黏性通量的貢獻(xiàn)；3) 通過解析數(shù)值通量格式(包括黏性通量)獲得與數(shù)值格式相對應(yīng)的解析精確Jacobian矩陣，對于諸如DG方法的緊致格式，可以通過鏈?zhǔn)椒▌t獲得與數(shù)值離散格式相對應(yīng)的解析精確Jacobian矩陣。前兩種處理方法對于常規(guī)二階有限體積方法是穩(wěn)定的，但是對于高精度數(shù)值格式而言，Jacobian矩陣的精確性直接決定了隱式數(shù)值格式的穩(wěn)定性和收斂性，因而采用與數(shù)值格式相對應(yīng)的精確Jacobian矩陣進(jìn)行隱式GMRES計(jì)算意義重大。

為了提高DG方法在采用隱式時(shí)間推進(jìn)格式進(jìn)行黏性計(jì)算的穩(wěn)定性和收斂性，本文提出了一種基于鏈?zhǔn)椒▌t解析求解與數(shù)值格式相對應(yīng)的精確Jacobian矩陣的GMRES隱式方法，該方法不僅有利于提高DG隱式求解的穩(wěn)定性，而且還能夠大幅度提高計(jì)算效率。與此同時(shí)，在RANS方程求解中，針對SA湍流模型源項(xiàng)隱式化時(shí)出現(xiàn)的非物理解，對模型的生成源項(xiàng)進(jìn)行修正，并給出合理的隱式化方式。文中通過層流和湍流算例驗(yàn)證了本文發(fā)展方法的有效性。

1 控制方程

笛卡爾坐標(biāo)系下耦合修正的負(fù)SA湍流模型的RANS方程組為[17]：

(1)

其中，U是守恒變量、Fi是無黏通量、Fvi是黏性通量和S是源項(xiàng)。

2 間斷伽遼金有限元方法

式(1)在弱積分下的DG方法可表示為：

(2)

黏性通量離散采用BR2格式[2]，該格式的形式如下：

(3)

為了保證格式的穩(wěn)定性，h取1～6。r為局部算子，R是全局算子，定義如下：

(4)

(5)

這兩種算子的離散形式為：

(6)

(7)

(8)

(9)

(10)

基函數(shù)Bi(x)采用Luo提出來的Taylor基函數(shù)[6]。需要說明的是文中的BR2P1是二階DG格式，BR2P2是三階DG格式。

3 隱式時(shí)間推進(jìn)格式

式(4)寫成半離散形式為：

(11)

其中，M是質(zhì)量矩陣,U是解矢量，R是右端殘值。R包括了無黏通量、黏性通量和它們對應(yīng)的域積分項(xiàng)。這里無黏通量的求解采用HLLC格式[19]，黏性通量的求解采用BR2格式[2]。

(12)

式(11)線化后可以表示為：

(13)

其中?Rn/?U是Jacobian矩陣，Δt是當(dāng)?shù)貢r(shí)間步長。當(dāng)Δt趨于無窮大時(shí)，系統(tǒng)方程的求解就變成了典型的牛頓迭代。

3.1 GMRES方法

GMRES在求解線性方程組時(shí)具有較高計(jì)算效率。由于本文計(jì)算所求解的控制方程非線性非常強(qiáng)，需要在GMRES求解時(shí)進(jìn)行預(yù)處理，本文采用LU-SGS進(jìn)行預(yù)處理[20]，具體如下：

P-1AΔU=P-1Rn

(14)

P是預(yù)處理矩陣，

P=(L+D)D-1(U+D)

(15)

其中，L、U和D分別是下三角、上三角和對角矩陣。

在GMRES的求解過程中，和Jacobian矩陣密切相關(guān)的兩個(gè)步驟是預(yù)處理和矩陣矢量生成。預(yù)處理可以采用近似的Jacobian進(jìn)行計(jì)算，而矩陣矢量生成則需要精確計(jì)算AΔU。傳統(tǒng)的處理方式是考慮到黏性通量和其數(shù)值格式的復(fù)雜性，通常采用無黏Jacobian矩陣，或者考慮到黏性通量的影響采用譜半徑代替黏性Jacobian矩陣進(jìn)行預(yù)處理。在矩陣矢量生成時(shí)必須采用有限差分進(jìn)行AΔU計(jì)算，具體如下：

(16)

其中ε是個(gè)小量，決定了式(16)的計(jì)算精度。如果能夠獲得與數(shù)值通量離散格式相對應(yīng)的精確Jacobian矩陣L、U和D，不但可以采用精確Jacobian矩陣進(jìn)行預(yù)處理，而且還可以采用精確Jacobian矩陣獲得更為精確的AΔU，即有:

(17)

對比近似Jacobian矩陣預(yù)處理加有限差分矩陣矢量生成與精確Jacobian矩陣進(jìn)行預(yù)處理和矩陣矢量生成，如果Jacobian矩陣不精確，雖然可以進(jìn)行預(yù)處理，但是無法直接進(jìn)行矩陣矢量生成，必須采用有限差分進(jìn)行矩陣矢量生成，這種方式一方面預(yù)處理的不精確影響GMRES的穩(wěn)定性和系統(tǒng)方程的收斂性，另一方面有限差分的計(jì)算精度取決于ε，容易影響AΔU的計(jì)算精度，與此同時(shí)采用有限差分增加了內(nèi)存存儲R(U+εΔU)。而采用精確Jacobian矩陣進(jìn)行計(jì)算只需要在求解GMRES之前計(jì)算一次L、U和D矩陣，且采用精確Jacobian矩陣預(yù)處理降低了系統(tǒng)剛性，提高了AΔU的計(jì)算精度，從而大幅度提高了GMRES的穩(wěn)定性和系統(tǒng)方程的收斂性。

3.2 解析精確Jacobian矩陣

像DG這類緊致格式任意高階數(shù)值通量所對應(yīng)的Jacobian矩陣都可以精確求出。對于采用HLLC格式求解無黏通量所對應(yīng)的解析精確Jacobian矩陣可以寫成：

(18)

對于采用BR2格式求解黏性通量所對應(yīng)的解析精確Jacobian矩陣可以寫成為:

(19)

其中face為總交界面,I和J代表界面左右的單元標(biāo)號。式(18)和式(19)中L和U矩陣以edge形式存儲信息，D矩陣則以element的形式存儲。式(18)和(19)的求解可以采用鏈?zhǔn)椒▌t求解：

(20)

(21)

其中?Uh/?U和?Ux/?U分別是B和Bx的對角陣。有關(guān)式(18)和(19)的詳細(xì)求解可見參考文獻(xiàn)[11]。需要說明的是，解析精確Jacobian的求解完全按照數(shù)值格式進(jìn)行鏈?zhǔn)椒▌t離散求解，中間沒有任何假設(shè)，從數(shù)值離散角度講，解析精確Jacobian矩陣精確地對應(yīng)了相應(yīng)的?R/?U，因而在每一步迭代計(jì)算中Jacobian矩陣是精確的。

3.3 SA湍流模型

對于SA模型而言，在源項(xiàng)中由于渦量容易為0，從而導(dǎo)致式(9)中的分母為0，因此對其進(jìn)行了修正，即：

(22)

因此，對應(yīng)的導(dǎo)數(shù)為：

(23)

修正后既保證了湍流生成源項(xiàng)P中渦量為0時(shí)不產(chǎn)生非物理解，又限制了湍流生成項(xiàng)在求解Jacobian矩陣時(shí)因分母為0而導(dǎo)致的計(jì)算失敗。

4 算例分析

4.1 層流圓柱繞流

層流圓柱繞流的計(jì)算條件為Ma=0.2、雷諾數(shù)Re=40。計(jì)算網(wǎng)格采用二階曲線三角形和四邊形混合網(wǎng)格，總共有8046個(gè)三角形單元和4568個(gè)四邊形單元，網(wǎng)格如圖1所示。網(wǎng)格離壁面最近距離為0.0338倍直徑，其中圓柱直徑為1。

圖2給出了CFL在30步之內(nèi)從1增加到1×1010時(shí)BR2P1和BR2P2的殘值收斂曲線。從圖可知，采用精確Jacobian矩陣進(jìn)行計(jì)算，在有限的計(jì)算步數(shù)內(nèi)CFL數(shù)可以取到很大的數(shù)，且BR2P1在迭代25步時(shí)殘值下降10個(gè)量級，相比而言BR2P2需要迭代32步。就計(jì)算時(shí)間而言，BR2P2的CPU耗時(shí)約是BR2P1的3倍以上。

圖3給出了層流圓柱的流線圖，從圖可知本文計(jì)算方法能夠準(zhǔn)確地捕捉到層流圓柱的一對尾渦，這和相關(guān)文獻(xiàn)結(jié)果[21]吻合。表1給出了計(jì)算的層流圓柱阻力系數(shù)和回流區(qū)長度與文獻(xiàn)值的對比。從表1可知本文計(jì)算結(jié)果和文獻(xiàn)值吻合較好。

SchemeForm dragCDpSkin friction dragCDvTotal dragCDRecirculation lengthBR2P11.01440.52471.53912.32BR2P21.01280.52381.53652.30Ref.

4.2 零壓力梯度平板湍流繞流

為了驗(yàn)證本文發(fā)展方法計(jì)算精度的網(wǎng)格收斂性，算例選擇零壓力梯度平板湍流繞流。計(jì)算條件為：馬赫數(shù)Ma=0.2，雷諾數(shù)Re=6×106。這一算例是NASA TMR網(wǎng)站[22]上的湍流模型驗(yàn)證算例，網(wǎng)格取自該網(wǎng)站，計(jì)算采用的三套網(wǎng)格大小分別為35×25，69×49和137×97，如圖4所示。

圖5給出了CFL在200步之內(nèi)從1增加到1×1010時(shí)BR2P1和BR2P2格式的殘值收斂曲線。由圖5可知，殘值下降10個(gè)量級，就迭代步數(shù)而言，BR2P2的迭代步數(shù)是BR2P1的1/3；就計(jì)算效率而言，BR2P2格式的CPU耗時(shí)是BR2P1格式的近3倍。圖6給出了兩種不同精度DG方法計(jì)算的阻力系數(shù)隨網(wǎng)格量的收斂狀況。由圖6可知本文發(fā)展的DG方法有較高的計(jì)算精度，在較粗的網(wǎng)上就能取得收斂解，特別是BR2P2格式能夠在最粗網(wǎng)格上獲得收斂解。與傳統(tǒng)的成熟CFD軟件CFL3D和FUN3D的計(jì)算結(jié)果相比較，本文發(fā)展DG方法的精算精度要明顯高于CFL3D和FUN3D，能夠在較稀的網(wǎng)格上得到收斂解，雖然最終收斂解的值有所差別，但是這種差別不超1%。

4.3 30P30N多段翼型湍流繞流

30P30N多段翼型湍流繞流是具有挑戰(zhàn)性的算例，這是由于在前緣縫翼縫道內(nèi)和襟翼艙內(nèi)有分離渦，加之在縫翼前緣存在跨聲速區(qū)，這些復(fù)雜的流動導(dǎo)致了計(jì)算的不穩(wěn)定。計(jì)算條件為：馬赫數(shù)Ma=0.2，迎角α=16°，雷諾數(shù)Re=9.0×106。計(jì)算網(wǎng)格采用二階曲線三角形和四邊形混合網(wǎng)格，總共有20 461個(gè)三角形單元和10 160個(gè)四邊形單元，第一層離壁面距離為5.0×10-6，遠(yuǎn)場距離為500倍弦長。網(wǎng)格示意圖如圖7所示。

由于30P30N多段翼型構(gòu)型和流動均非常復(fù)雜。為了提高計(jì)算的穩(wěn)定性和計(jì)算效率，先采用較小CFL數(shù)進(jìn)行計(jì)算使流場達(dá)到穩(wěn)定后再采用快速增加CFL數(shù)到1×1010。所以計(jì)算CFL數(shù)設(shè)定遵從如下方式：首先在1000步以內(nèi)CFL數(shù)從1增大到1000，當(dāng)殘值下降2.5個(gè)量級時(shí)CFL數(shù)在接下來的100步以內(nèi)從當(dāng)前值增大到1×1010。圖8給出了計(jì)算的殘值收斂曲線。由圖8可以看出，兩種不同精度的BR2格式殘值收斂歷程大體相同，BR2P2的CPU耗時(shí)是BR2P1的3倍。需要指出的是BR2P2計(jì)算時(shí)的初始值是以自由來流且沒有采用限制器和人工耗散，計(jì)算依然能夠在大CFL數(shù)下穩(wěn)定并且殘值最終下降了10個(gè)量級。這進(jìn)一步驗(yàn)證了本文發(fā)展的高精度方法具有較高的穩(wěn)定性，能夠用于解決諸如30P30N這樣復(fù)雜構(gòu)型的氣動計(jì)算。

圖9給出了BR2P1和BR2P2計(jì)算的30P30N多段翼型表面壓力系數(shù)分布和DLR結(jié)果[23]的對比。由圖9可知，本文計(jì)算結(jié)果和DLR的計(jì)算值吻合較好。表2給出了BR2P1和BR2P2計(jì)算的升阻力系數(shù)和文獻(xiàn)值[24]的對比(其中Ndof為自由度個(gè)數(shù))，可以看出，然本文的網(wǎng)格單元更少，但是本文的計(jì)算結(jié)果和文獻(xiàn)結(jié)果吻合度較高。

SchemeNdofCLCDBR2P19,18634.098760.054366BR2P2182,7264.109420.052248Ref. DGP1[24]195,8994.122890.052933Ref. DGP2[24]4198054.102820.052202

5 結(jié) 論

為了解決DG在RANS計(jì)算中的穩(wěn)定性問題，本文發(fā)展了一種基于精確Jacobian矩陣的高效隱式間斷伽遼金方法，用于求解可壓縮流動問題。通過層流和湍流算例驗(yàn)證了發(fā)展的方法的有效性。典型算例結(jié)果表明：

1) 發(fā)展的基于鏈?zhǔn)椒▌t求解精確Jacobian矩陣的隱式GMRES方法，不僅能夠提高計(jì)算的穩(wěn)定性，而且還能夠提高計(jì)算效率，可以實(shí)現(xiàn)超大CFL數(shù)的計(jì)算。

2) 針對SA湍流模型生成源項(xiàng)的修正較為有效，能夠改善RANS-SA系統(tǒng)方程在采用隱式方法計(jì)算時(shí)的穩(wěn)定性。

3) 發(fā)展的高精度間斷伽遼金方法不僅能夠用于層流流場的計(jì)算，而且還能夠用于復(fù)雜湍流流場的數(shù)值模擬，結(jié)果和國際上相關(guān)成熟軟件吻合較好。這表明所發(fā)展的高精度間斷伽遼金方法能夠用于求解復(fù)雜二維流動問題，具有較為廣闊應(yīng)用前景。

參 考 文 獻(xiàn)：

[1]Reed W H, Hill T R. Triangular mesh methods for the neutron transport[R]. Los Alamos Scientific Laboratory Report, LAUR-73-479, 1973.

[2]Cockburn B, Hou S, Shu C W. TVD Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws IV: the multidimensional case[J]. Mathematics of Compution, 1990, 54(190): 545-581.

[3]Bassi F, Rebay S. Discontinuous Galerkin solution of the Reynolds-averaged Navier-Stokes andk-εturbulence model equations[J]. Computers & Fluids, 2005, 34: 507-540.

[4]Yang X Q, Cheng J, Liu X D, et al. A Reconstructed direct discontinuous Galerkin method for compressible turbulent flows on Hybrid Grids[R]. AIAA 2016-3332, 2016.

[5]Cheng J, Yang X Q, Liu T G, et al. The direct discontinuous Galerkin method for 2D compressible Navier-Stokes equations on arbitrary grid[R]. AIAA 2016-1334, 2016.

[6]Luo H, Luo L Q, Nourgaliev R. A reconstructed discontinuous Galerkin method for the compressible Euler equations on arbitrary grids[J]. Communication in Computational Physcis, 2012, 12: 1495-1519.

[7]Arnold D N. An interior penalty finite element method with discontinuous elements[J]. SIAM Journal on Numerical Analysis, 1982, 19(4):742-760.

[8]Cockburn B, Shu C W. The local discontinuous Galerkin method for time dependent convection diffusion systems[J]. SIAM Journal on Numerical Analysis, 1998, 35(6): 2440-2463.

[9]Peraire J, Persson P O. The compact discontinuous Galerkin (CDG) method for elliptic problems[J]. SIAM J Sci Comput 2008, 30: 1806-1824.

[10]Liu H L. Optimal error estimates of the direct discontinuous Galerkin method for convection-diffusion equations[J]. Math Comput, 2015, 84: 2263-2295.

[11]Yang X Q, Cheng J, Wang C J, et al. A fast, implicit discontinuous Galerkin method based on analytical jacobians for the compressible Navier-Stokes equations[R]. AIAA 2016-1326, 2016.

[12]Cheng J, Yang X Q, Liu T G, et al. A direct discontinuous Galerkin method for 2D compressible Navier-Stokes equations on arbitrary grid[J]. Journal of Computational Physics, 2016, 327: 484-502.

[13]Dumbser M, Valsara D S, Toro E F, et al. A unified framework for the construction of one-step finite volume and discontinuous Galerkin schemes on unstructured meshes[J]. Journal of Computational Physics, 2008, 227: 8209-8253.

[14]Van L B, Nomura S. Discontinuous Galerkin method for diffusion[R]. AIAA 2005-5108, 2005.

[15]Li M, Liu W, Zhang L P, et al. Applications of high order hybrid DG/FV methods for 2D viscous flows[J]. Acta Aerodynamica Sinica, 2015, 33(1): 17-24. (in Chinese)李明, 劉偉, 張來平, 等. 高階精度DG/FV混合方法在二維黏性流動模擬中的推廣[J]. 空氣動力學(xué)學(xué)報(bào), 2015, 33(1): 17-24.

[16]Yang X Q, Yang A, Sun G. An efficient numerical method for coupling the RANS equations with Spalart-Allmaras turbulence model equation[J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(9): 2007-2018. (in Chinese)楊小權(quán), 楊愛明, 孫剛. 一種強(qiáng)耦合Spalart-Allmaras 湍流模型的RANS 方程的高效數(shù)值計(jì)算方法[J]. 航空學(xué)報(bào), 2013, 34(9): 2007-2018.

[17]Allmaras S R, Johnson F T, Spalart P R. Modifications and clarifications for the implementation of the Spalart-Allmaras turbulence model[C]//Seventh International Conference on Computational Fluid Dynamics, 2012.

[18]Todd A O. A high-order, adaptive, discontinuous Galerkin——nite element method for the Reynolds averaged Navier-Stokes equations[D]. Massachusetts Institute of Technology, 2008.

[19]Toro E F. Riemann solvers and numerical methods for fluid dynamics[M]. Verlag, Beilin: Springer Press, 1999.

[20]Luo H, Baum J D, L?hner R. A fast, matrix-free implicit method for compressible flows on unstructured grids[J]. J Comput Phys, 1998, 146: 664-690.

[21]Tseng Y H, Ferziger J H. A ghost-cell immersed boundary method for flow in complex geometry[J]. Journal of Computational Physics, 2003, 192: 593-623.

[22]https://turbmodels. larc.nasa.gov/[23]http://www.ipacs-benchmark. org/[24]Peraire J, Persson P O. The compact discontinuous Galerkin (CDG) method for elliptic problems[J]. SIAM Journal on Scientific Computing, 2008, 30(4): 1806-1824.