曹貴瑜, 潘 亮, 徐 昆

(1. 南方科技大學(xué) 前沿與交叉科學(xué)研究院, 深圳 518055; 2. 北京師范大學(xué) 數(shù)學(xué)科學(xué)學(xué)院, 數(shù)學(xué)與復(fù)雜系統(tǒng)教育部重點(diǎn)實(shí)驗(yàn)室, 北京 100875;3. 香港科技大學(xué) 數(shù)學(xué)系, 香港)

0 引 言

湍流在自然現(xiàn)象和工程應(yīng)用中無(wú)處不在[1]，湍流問(wèn)題的研究對(duì)于航空、航天、天氣預(yù)報(bào)等領(lǐng)域的國(guó)家戰(zhàn)略需求和學(xué)科發(fā)展具有重要意義。目前，難以采用理論解析方法直接研究湍流問(wèn)題，因此通常需要借助實(shí)驗(yàn)手段和數(shù)值模擬進(jìn)行研究。

隨著數(shù)值計(jì)算方法和超級(jí)計(jì)算機(jī)的發(fā)展，直接數(shù)值模擬[2-4](Direct Numerical Simulation，DNS)逐漸成為研究湍流機(jī)理的重要手段。對(duì)于不可壓縮湍流的直接數(shù)值模擬，譜方法[3](Spectral Method)、偽譜法[5](Pesudo-Spectral Method)以及格子玻爾茲曼方法[6](Lattice Boltzmann Method，LBM)得到了充分的發(fā)展。但是，以上方法都不能應(yīng)用于帶間斷的可壓縮湍流模擬。目前，高階有限差分法[7-8]被廣泛應(yīng)用于可壓縮湍流的直接數(shù)值模擬，但由于其數(shù)值不穩(wěn)定性，仍難以適用于超聲速湍流馬赫數(shù)流動(dòng)。為了在間斷區(qū)域避免非物理振蕩、光滑流動(dòng)區(qū)域保持高精度，Wang等結(jié)合WENO重構(gòu)和緊致有限差分格式發(fā)展了高精度混合格式[9]，并成功應(yīng)用于超聲速湍流馬赫數(shù)的各向同性湍流直接數(shù)值模擬中。鑒于有限體積格式的物理守恒性，且能適用于復(fù)雜網(wǎng)格和復(fù)雜幾何外形，二階有限體積格式被廣泛應(yīng)用于可壓縮工程湍流的計(jì)算中[10]。但是，對(duì)可壓縮湍流的直接數(shù)值模擬卻鮮有報(bào)道。在可壓縮湍流場(chǎng)中通常包含激波、邊界層、渦以及它們之間的相互作用等復(fù)雜的流動(dòng)結(jié)構(gòu)。二階精度的數(shù)值方法具有較大的數(shù)值耗散和數(shù)值色散，難以捕捉精細(xì)的流場(chǎng)結(jié)構(gòu)。近年來(lái)，高階有限體積格式逐漸興起，克服了二階格式的精度不足問(wèn)題，逐漸成為可壓縮湍流直接數(shù)值模擬的重要研究工具。

在過(guò)去的三十年里，氣體動(dòng)理學(xué)格式[11-12](Gas-Kinetic Scheme, GKS)得到了系統(tǒng)性的發(fā)展，目前已成功應(yīng)用于從低速到高超聲速的全速域流動(dòng)問(wèn)題的數(shù)值模擬。在有限體積格式的框架下，基于Bhatnagar-Gross-Krook(BGK)方程[13]的積分解直接構(gòu)造網(wǎng)格界面上的分布函數(shù)，求矩得到相應(yīng)的數(shù)值通量進(jìn)而更新下一時(shí)刻網(wǎng)格單元內(nèi)的宏觀守恒量。相比于傳統(tǒng)的Riemann求解器，可以同時(shí)考慮從自由分子流到連續(xù)流的多尺度演化過(guò)程，不僅能準(zhǔn)確給出流場(chǎng)光滑區(qū)域的解，而且能在間斷區(qū)域很好地捕捉激波?；贑hanman-Enskog(C-E)展開(kāi)，可以同時(shí)處理無(wú)黏和黏性問(wèn)題，避免額外的黏性項(xiàng)離散。此外，氣體動(dòng)理學(xué)格式是真正的多維格式，可以將切向信息包含在數(shù)值通量中，這對(duì)于三維流動(dòng)和間斷流動(dòng)十分重要。近年來(lái)，高階精度氣體動(dòng)理學(xué)格式(High-Order Gas-Kinetic Scheme, HGKS)發(fā)展也極為迅速?；赪ENO空間重構(gòu)[14-15]和速度分布函數(shù)高階展開(kāi)，發(fā)展了單步三階的氣體動(dòng)理學(xué)格式[16-17]。該格式可避免使用傳統(tǒng)高精度格式中界面上通量的Gauss積分和Runge-Kutta時(shí)間離散。但是，高階分布函數(shù)的形式十分復(fù)雜，給發(fā)展更高階數(shù)值格式和三維實(shí)戰(zhàn)編程帶來(lái)巨大的困難。相比于時(shí)間空間解耦的Riemann求解器，時(shí)空耦合的氣體動(dòng)理學(xué)求解器可以提供通量的時(shí)間導(dǎo)數(shù)。基于兩步四階時(shí)間離散和高階精度空間重構(gòu)[18]，我們發(fā)展了具有四階時(shí)間精度的氣體動(dòng)理學(xué)格式[19]。與原有的單步高階氣體動(dòng)理學(xué)方法相比，兩步四階時(shí)間離散方法僅需要二階GKS通量，通量形式得到大大簡(jiǎn)化，大大降低計(jì)算量。目前，結(jié)構(gòu)網(wǎng)格和非結(jié)構(gòu)網(wǎng)格上的兩步四階氣體動(dòng)理學(xué)格式都得到了充分發(fā)展，并對(duì)無(wú)黏流和層流的標(biāo)準(zhǔn)算例都進(jìn)行了驗(yàn)證[20-21]。該格式不僅具有更高的數(shù)值精度和更好的穩(wěn)定性，而且有處理復(fù)雜流動(dòng)問(wèn)題的能力。

近幾年，GKS也不斷被國(guó)內(nèi)外研究者應(yīng)用于湍流的數(shù)值模擬。對(duì)于低雷諾數(shù)湍流問(wèn)題，二階GKS和基于WENO重構(gòu)的空間高分辨率GKS已被用作直接數(shù)值模擬工具[22-23]。對(duì)于高雷諾數(shù)工程湍流問(wèn)題，通過(guò)適當(dāng)?shù)亟o定BGK方程湍流弛豫時(shí)間[24]能引入湍流黏性，可以將氣體動(dòng)理學(xué)格式與渦黏湍流模式耦合，實(shí)現(xiàn)不可解網(wǎng)格上RANS模擬。目前，耦合各類(lèi)RANS模型、LES模型和混合湍流模型[25-30]的二階GKS和高分辨率GKS已被國(guó)內(nèi)外多個(gè)課題組成功應(yīng)用于高雷諾數(shù)工程湍流問(wèn)題的計(jì)算中?？紤]到湍流數(shù)值模擬對(duì)空間和時(shí)間都要求高分辨率，例如直接數(shù)值模擬需要解析Kolmogorov時(shí)間尺度和空間尺度上的所有湍流結(jié)構(gòu)，將高精度高數(shù)值穩(wěn)定性的兩步四階HGKS推廣到湍流直接數(shù)值模擬和RANS計(jì)算中十分必要。鑒于此，本文回顧HGKS在湍流直接數(shù)值模擬和RANS模擬中的應(yīng)用。

1 高階氣體動(dòng)理學(xué)格式

對(duì)于有限體積方法，關(guān)鍵過(guò)程是通過(guò)數(shù)值通量更新每個(gè)控制體內(nèi)的守恒變量。GKS的數(shù)值通量是基于BGK方程的積分解。無(wú)外力場(chǎng)的三維BGK方程為:

式中碰撞不變量矢量

相空間微元為dΞ=dudvdwdξ1…dξN。在連續(xù)流區(qū)域，對(duì)BGK方程做C-E展開(kāi)，

f=g-τDug+τDu(τDug)+…

(4)

式中Du=?/?t+ui?/?xi。根據(jù)氣體動(dòng)理論，零階展開(kāi)f=g對(duì)應(yīng)推導(dǎo)Euler方程；一階展開(kāi)f=g-τDug對(duì)應(yīng)推導(dǎo)NS方程。這就是BGK方程求解NS方程的基礎(chǔ)。GKS就是通過(guò)求解BGK方程式(1)達(dá)到求解NS方程的目的。

BGK方程對(duì)碰撞不變量ψ求矩，并在有限控制體上離散：

式中，|Ωijk|為有限控制體體積。以x方向上數(shù)值通量為例，

式中，xi+1/2,jm,kn=(xi+1/2,yjm,zkn)T，(yjm,zkn)是在界面yj×zk上的高斯點(diǎn)坐標(biāo)，ωmn是高斯點(diǎn)權(quán)重。對(duì)于HGKS，核心任務(wù)就是獲得界面高斯點(diǎn)上的分布函數(shù)f(xi+1/2,jm,kn,t,u,ξ)，進(jìn)而通過(guò)取矩求得界面通量以更新下一時(shí)刻單元內(nèi)的宏觀守恒量。BGK方程局部形式解提供了界面上的分布函數(shù)：

f(xi+1/2,jm,kn,t,u,ξ)=

(8)

式中：x′=xi+1/2,jm,kn-u(t-t′)是分子運(yùn)動(dòng)軌跡，f0是界面初始速度分布函數(shù)，g是界面上沿分子運(yùn)動(dòng)軌跡對(duì)應(yīng)的平衡態(tài)分布函數(shù)?？梢钥吹剑@個(gè)積分解包含了時(shí)間和空間的所有信息。經(jīng)過(guò)簡(jiǎn)單計(jì)算，界面分布函數(shù)顯式表達(dá)為：

f(xi+1/2,jm,kn,t,u,ξ)=(1-e-t/τ)g0+

e-t/τgl[1-(τ+t)(alu+blv+clw)-τAl)]H(u)+

e-t/τgr[1-(τ+t)(aru+brv+crw)-

τAr)](1-H(u))

(9)

其中H(…)為Heaviside函數(shù)，gl和gr對(duì)應(yīng)界面兩側(cè)的初始平衡態(tài)分布函數(shù)，g0是界面上初始平衡態(tài)分布。g0通過(guò)如下求得。

〈aku+bkv+ckw+Ak〉=0

為了獲得時(shí)空上的高階精度，高階有限體積格式需要高階空間重構(gòu)和多步時(shí)間離散?；诟呔瓤臻g重構(gòu)，可以完全確定界面上的分布函數(shù)如式(11)，對(duì)其求矩可以得到界面上時(shí)空耦合的數(shù)值通量。為了獲得高階空間精度，在結(jié)構(gòu)網(wǎng)格的計(jì)算中，采用經(jīng)典的五階WENO重構(gòu)，具體細(xì)節(jié)見(jiàn)參見(jiàn)文獻(xiàn)[14-15]。對(duì)于二維和三維問(wèn)題，采用一維一維重構(gòu)方式，來(lái)獲得界面高斯點(diǎn)值和多維空間梯度值，進(jìn)而獲得時(shí)空耦合的數(shù)值通量。為獲得時(shí)間高階精度，基于守恒律方程的Lax-Wendroff形求解器，Li等發(fā)展了兩步四階的時(shí)間離散方法[18-19]。對(duì)于時(shí)空耦合的演化過(guò)程，可以利用數(shù)值通量和數(shù)值通量的時(shí)間導(dǎo)數(shù)值實(shí)現(xiàn)高階時(shí)間精度。對(duì)于半離散有限體積格式式(5)，可以采用以下的兩步格式進(jìn)行離散：

其中Q*為t*=tn+Δt/2時(shí)刻的守恒量?？梢宰C明，由式(12)和式(13)給出的Qn+1具有四階時(shí)間精度。和四階Runge-Kutta方法相比，僅需要一個(gè)中間步就可實(shí)現(xiàn)四階精度，提高了數(shù)值格式的計(jì)算效率[20]。更重要的是，兩步四階格式有著和二階格式相當(dāng)?shù)臄?shù)值穩(wěn)定性。至此，基于二階GKS通量，五階WENO空間重構(gòu)和兩步時(shí)間離散的兩步四階HGKS簡(jiǎn)述完畢。

2 數(shù)值算例

本節(jié)中，選取經(jīng)典的無(wú)黏算例、層流算例、低雷諾數(shù)湍流直接數(shù)值模擬算例和高雷諾數(shù)工程湍流RANS模擬算例，來(lái)驗(yàn)證HGKS計(jì)算復(fù)雜流動(dòng)問(wèn)題的能力。

2.1 Titarev-Toro問(wèn)題

第一個(gè)算例是Titarev-Toro 問(wèn)題[31]，初值條件如下：

這個(gè)算例描述激波和高頻振蕩波的相互作用。左邊界采用無(wú)反射邊界條件，右邊界給定固定的初值波形。采用1000個(gè)均勻網(wǎng)格進(jìn)行計(jì)算，并基于兩步四階、兩步五階、三步五階氣體動(dòng)理學(xué)格式和WENO重構(gòu)測(cè)試該算例[20]。在圖1中給出Titarev-Toro問(wèn)題在t=1.8時(shí)刻的數(shù)值解和精確解。相比與傳統(tǒng)的基于Runge-Kutta時(shí)間離散和Riemann求解器的數(shù)值格式相比，HGKS能很好地捕捉流場(chǎng)中的高頻波。這也證明了時(shí)間空間耦合的數(shù)值通量對(duì)于高階格式的重要性。

圖1 Titarev-Toro問(wèn)題，t=5的密度數(shù)值解和精確解[20]

2.2 雙馬赫反射問(wèn)題

雙馬赫反射問(wèn)題是可壓縮無(wú)黏流動(dòng)的經(jīng)典算例[32]。該算例的計(jì)算域?yàn)閇0,4]×[0,1]。初始時(shí)刻x=1/6處有一道與x軸呈60°的馬赫數(shù)為10的激波，激波前后的初值條件如下：

(ρ,U,V,p)=(1.4,0,0,1)

從x=1/6開(kāi)始壁面上采用反射邊條件，其余部分的下邊界和左側(cè)邊界一樣采用來(lái)流條件，右側(cè)邊界采用出口條件，并根據(jù)激波的移動(dòng)位置給出上邊界的邊界條件。圖2中給出WENO-Z格式采用1440×480均勻網(wǎng)格時(shí)的局部密度分布。數(shù)值結(jié)果表明HGKS有很好的健壯性，并且能夠很好地分辨剪切層的不穩(wěn)定性。

圖2 雙馬赫反射問(wèn)題，t=0.2的數(shù)值解[19]

2.3 平板邊界層問(wèn)題

二維平板邊界層問(wèn)題是低速黏性流的經(jīng)典算例。算例中，來(lái)流馬赫數(shù)Ma=0.1，雷諾數(shù)Re=1×105。壁面上采用無(wú)滑移條件，平板之前采用對(duì)稱(chēng)條件，其它邊界根據(jù)Riemann不變量給定無(wú)反射條件。網(wǎng)格數(shù)為80×40的非均勻網(wǎng)格。圖3中給出無(wú)量綱的U和V，數(shù)值結(jié)果和Blasius 參考解吻合得很好。HGKS僅用五個(gè)網(wǎng)格就能解析邊界層速度型。

圖3 平板邊界層問(wèn)題，無(wú)量綱的U和V，數(shù)值結(jié)果和Blasius 參考解[19]

2.4 低速槽道湍流

槽道湍流是研究壁湍流的典型算例。低速槽道湍流計(jì)算條件：馬赫數(shù)Ma=0.1，壁面摩擦雷諾數(shù)Reτ=180。計(jì)算域?yàn)閇0,2π]×[-1,1]×[0,π]，G1和G2兩套計(jì)算網(wǎng)格分別為963和1283。初始速度場(chǎng)采用泊肅葉理論解疊加幅值為當(dāng)?shù)厮俣?0%的白噪聲隨機(jī)場(chǎng)。通過(guò)外加質(zhì)量力驅(qū)動(dòng)流動(dòng), 保證流量為常數(shù)。流向和展向給定為周期邊界條件，法向壁面為等溫壁。

圖4(a)中給出了歸一化的流向平均速度分布。其中，作為參考解的譜方法計(jì)算網(wǎng)格為129×192×160。LBM作為模擬不可壓縮流動(dòng)流行的介觀方法，其計(jì)算網(wǎng)格為200×400×200，DUGKS則采用1283計(jì)算網(wǎng)格[33]?？梢钥闯鯤GKS、LBM、DUGKS和參考解譜方法[2]計(jì)算結(jié)果都符合得很好。圖4(b)顯示了平均雷諾剪切應(yīng)力曲線。與譜方法結(jié)果相比，G2網(wǎng)格上HGKS和DUGKS結(jié)果明顯優(yōu)于LBM結(jié)果，尤其是在近壁區(qū)域。

(a) 平均流向速度

圖5給出了用壁面摩擦速度無(wú)量綱化的三個(gè)方向的脈動(dòng)速度均方根值。由圖可見(jiàn)，流向速度脈動(dòng)峰值明顯大于另外兩個(gè)方向的脈動(dòng)，且在槽道中心區(qū)域湍流脈動(dòng)趨于各向同性。從圖5還可以看出，LBM在近中心線區(qū)域表現(xiàn)較好而在近壁區(qū)欠佳，DUGKS在近壁區(qū)域表現(xiàn)較好而在近中心線區(qū)域欠佳。HGKS在近壁區(qū)域和近中心線區(qū)域都與譜方法結(jié)果更接近，優(yōu)于二階LBM和DUGKS。此外，與LBM相比，HGKS結(jié)果是在較大計(jì)算域的粗網(wǎng)格分辨率獲得的。由于LBM需要等距網(wǎng)格，因此網(wǎng)格被限制為一個(gè)極小的值以求解黏性底層，即在相同尺寸的計(jì)算域中，LBM所需的網(wǎng)格數(shù)將遠(yuǎn)超HGKS。

(a) 流向脈動(dòng)速度均方根

2.5 超聲速均勻各向同性衰減湍流

可壓縮均勻各向同性衰減湍流常用來(lái)檢驗(yàn)湍流直接數(shù)值模擬中數(shù)值格式的魯棒性。算例R1和R2的初始泰勒微尺度雷諾數(shù)分別為Reλ=72和Reλ=120，初始湍流馬赫數(shù)均為Mat=2.0，動(dòng)力黏性系數(shù)按冪律指數(shù)為0.76給出。計(jì)算域?yàn)閇0,2π]×[0,2π]×[0,2π]，R1計(jì)算網(wǎng)格為3843，R2計(jì)算網(wǎng)格為5123，HGKS網(wǎng)格的分辨率和時(shí)間步長(zhǎng)以參考準(zhǔn)則[35]設(shè)定。初始速度脈動(dòng)場(chǎng)按無(wú)散條件給定，能譜為:

式中，A0=0.00013確定初始湍動(dòng)能，κ是波數(shù)，κ0=8確定能譜峰值。密度和壓力場(chǎng)為均勻分布。三個(gè)方向均采用周期邊界條件。

圖6(a)給出了無(wú)量綱時(shí)間t/τto=0.5和t/τto=1.0(τto為大渦翻轉(zhuǎn)時(shí)間)速度場(chǎng)脹量θ概率密度分布函數(shù)。所有脹量PDF都顯示出負(fù)尾巴，這是由激波結(jié)構(gòu)導(dǎo)致的可壓縮各向同性湍流中最顯著的流動(dòng)結(jié)構(gòu)。圖6(b)給出了無(wú)量綱時(shí)間t/τto=0.5和t/τto=1.0的脹量和x方向速度分量在x=0和z=0截線分布。截線預(yù)示，較強(qiáng)的可壓縮區(qū)域和高膨脹區(qū)域頻繁且隨機(jī)出現(xiàn)在湍流場(chǎng)中。為了更直觀地顯示流場(chǎng)結(jié)構(gòu)，圖7給出了無(wú)量綱時(shí)刻t/τto=0.5時(shí)θ/〈θ〉*云圖。其中〈θ〉*是流場(chǎng)脹量均方根。通常，θ/〈θ〉*≤-3的強(qiáng)壓縮區(qū)域被識(shí)別為小激波[36]。這些隨機(jī)分布的小激波和高膨脹區(qū)域?qū)е铝鲌?chǎng)中的劇烈變化的空間梯度和時(shí)間梯度，給高階數(shù)值格式穩(wěn)定性帶來(lái)巨大的挑戰(zhàn)。很少有高階數(shù)值格式能模擬超聲速湍流馬赫數(shù)的均勻各向同性衰減湍流。HGKS在高達(dá)Mat=2.0的高湍流馬赫數(shù)的DNS成功應(yīng)用證實(shí)了HGKS優(yōu)秀的魯棒性，為高可壓縮湍流提供了有力的計(jì)算工具。

(a) 速度場(chǎng)脹量θ概率密度分布函數(shù)

(b) 脹量θ和x方向速度分量在x=0和z=0截線分布

圖7 超聲速各向同性衰減湍流算例R1在無(wú)量綱t/τto=0.5時(shí)間θ/〈θ〉*云圖[34]

2.6 亞聲速NACA0012翼型湍流

亞聲速NACA0012翼型湍流是NASA湍流模擬網(wǎng)站[37]提供的標(biāo)準(zhǔn)湍流模型驗(yàn)證算例。自由來(lái)流條件如下：馬赫數(shù)Ma=0.15，雷諾數(shù)Re=3.0×106，翼型攻角15°，弦長(zhǎng)c=1.0，升力系數(shù)和阻力系數(shù)計(jì)算中參考面積為1.0。 計(jì)算域和邊界條件與NASA網(wǎng)站保持一致，其中G3和G4兩套計(jì)算網(wǎng)格設(shè)置分別為225×65和897×257。粗網(wǎng)格G3被應(yīng)用于隱式HGKS(Implicit HGKS, IHGKS)和二階隱式GKS(Implicit GKS, IGKS)模擬，參考解為CFL3D在密網(wǎng)格G4計(jì)算結(jié)果。對(duì)于RANS計(jì)算，式(1)中分子弛豫時(shí)間τ調(diào)整為:

其中μt為湍流模型[30]提供的湍流黏性系數(shù)。本文耦合k-ωSST 湍流模型求得湍流黏性系數(shù)，進(jìn)而調(diào)整BGK方程中的湍流弛豫時(shí)間，時(shí)間離散采用LUSGS隱式算法。

圖8中顯示了NACA0012翼型周?chē)膲毫ο禂?shù)分布和翼型上壁面摩擦系數(shù)分布。粗網(wǎng)格G3上，相比于二階隱式GKS的結(jié)果，隱式HGKS的壓力系數(shù)和摩擦系數(shù)與實(shí)驗(yàn)數(shù)據(jù)[38]和CFL3D在密網(wǎng)格G4上參考解更加接近。阻力系數(shù)CD對(duì)翼型周?chē)膲毫Ψ植己捅诿婺Σ亮Ψ植挤浅Ｃ舾小?如表1所示，隱式HGKS的升力系數(shù)CL和阻力系數(shù)CD在粗網(wǎng)格G3上求解結(jié)果非常接近密網(wǎng)格G4上CFL3D的參考解。 特別是粗網(wǎng)格G3上隱式HGKS的CD和參考解的差異在3個(gè)阻力數(shù)(0.0001)之內(nèi)，達(dá)到了工程湍流模擬的要求。 對(duì)于二階隱式GKS，升力系數(shù)是可以接受的，而它過(guò)高地預(yù)測(cè)了阻力系數(shù)達(dá)到40個(gè)阻力數(shù)。 該算例證實(shí)了隱式HGKS高精度RANS模擬湍流場(chǎng)的能力。

(a)

表1 NACA0012升力系數(shù)CL和阻力系數(shù)CD

2.7 跨聲速ARA M100翼身湍流

三維跨聲速ARA M100翼身湍流來(lái)流條件：馬赫數(shù)Ma=0.8027，基于弦長(zhǎng)雷諾數(shù)Relc=1.31×107(弦長(zhǎng)lc=0.245 )，翼身攻角2.873°。采用CFL3D第6版網(wǎng)站[39]提供的計(jì)算域和網(wǎng)格，其中C-O型網(wǎng)格設(shè)置為321×57×49。 圖9顯示了ARA M100機(jī)翼機(jī)身，其中黑色部分為機(jī)翼，綠色部分為機(jī)身。

圖9 ARA M100翼身構(gòu)型[34]

圖10(a)給出了翼尖附近截面Z/b=0.935的馬赫數(shù)的云圖，該云圖顯示激波-邊界層相互作用，證實(shí)了當(dāng)前隱式HGKS在激波捕捉方面的魯棒性。CFL3D網(wǎng)站的實(shí)驗(yàn)數(shù)據(jù)、隱式HGKS、二階隱式GKS和基于SA模型的CFL3D的在翼根附近典型截面Z/b=0.123壓力系數(shù)Cp的結(jié)果比較如圖9(b)所示。所有方法在該截面壓力系數(shù)Cp都與實(shí)驗(yàn)數(shù)據(jù)吻合較好。以上結(jié)果表明相比于數(shù)值離散的誤差，湍流模型的誤差在跨聲速三維復(fù)雜RANS模擬中占主導(dǎo)地位。該算例驗(yàn)證了隱式HGKS的魯棒性和模擬三維工程湍流的能力。

(a)

3 結(jié) 論

本文回顧了高階氣體動(dòng)理學(xué)格式及其在湍流數(shù)值模擬中的應(yīng)用研究。時(shí)空耦合的氣體動(dòng)理學(xué)求解器可以提供通量的時(shí)間導(dǎo)數(shù)?；趦刹剿碾A時(shí)間離散和高精度空間重構(gòu)，我們發(fā)展了具有四階時(shí)間精度的氣體動(dòng)理學(xué)格式。該格式不僅具有更高的數(shù)值精度和更好的穩(wěn)定性，而且有處理復(fù)雜流動(dòng)問(wèn)題的能力。數(shù)值結(jié)果表明，高階氣體動(dòng)理學(xué)格式可以為可壓縮湍流的數(shù)值模擬提供有力的計(jì)算工具。未來(lái)，將使用高階氣體動(dòng)理學(xué)格式研究更具有挑戰(zhàn)性的可壓縮湍流問(wèn)題，例如超聲速湍流邊界層和激波邊界層相互作用等。希望為可壓縮強(qiáng)間斷湍流的流動(dòng)機(jī)理研究探索新方向。

同時(shí)值得指出的是，氣體動(dòng)理學(xué)格式的主要優(yōu)勢(shì)是給出了在網(wǎng)格邊界上氣體分布函數(shù)的演化解。除了給出對(duì)應(yīng)的數(shù)值通量外，在網(wǎng)格邊界上的其它物理量也能夠顯示地給出來(lái)，比如網(wǎng)格邊界上在這一時(shí)間步長(zhǎng)結(jié)束時(shí)的守恒量。所以格式除了更新網(wǎng)格里面的守恒量的值，每個(gè)網(wǎng)格內(nèi)的守恒量的梯度也可以根據(jù)高斯定律由邊界上的值積分直接確定。這一性質(zhì)對(duì)構(gòu)造緊致高階氣體動(dòng)理學(xué)格式提供了可靠的網(wǎng)格局部的動(dòng)力學(xué)基礎(chǔ)，從而避免了用很遠(yuǎn)處和本網(wǎng)格沒(méi)有直接動(dòng)力學(xué)關(guān)系的量來(lái)做高階重構(gòu)。這方面的研究在過(guò)去幾年得到了飛速發(fā)展[21]。這一基于局部高階動(dòng)力學(xué)演化解的指導(dǎo)方向是發(fā)展下一代高階計(jì)算流體力學(xué)格式的可行之道。

致謝:感謝廣州天河超算中心提供計(jì)算資源。