高階氣體動理學格式在湍流數值模擬中的應用

曹貴瑜, 潘 亮, 徐 昆

(1. 南方科技大學 前沿與交叉科學研究院, 深圳 518055; 2. 北京師范大學 數學科學學院, 數學與復雜系統教育部重點實驗室, 北京 100875;3. 香港科技大學 數學系, 香港)

0 引 言

湍流在自然現象和工程應用中無處不在[1]，湍流問題的研究對于航空、航天、天氣預報等領域的國家戰略需求和學科發展具有重要意義。目前，難以采用理論解析方法直接研究湍流問題，因此通常需要借助實驗手段和數值模擬進行研究。

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

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

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

1 高階氣體動理學格式

對于有限體積方法，關鍵過程是通過數值通量更新每個控制體內的守恒變量。GKS的數值通量是基于BGK方程的積分解。無外力場的三維BGK方程為:

式中碰撞不變量矢量

相空間微元為dΞ=dudvdwdξ1…dξN。在連續流區域，對BGK方程做C-E展開，

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

(4)

式中Du=?/?t+ui?/?xi。根據氣體動理論，零階展開f=g對應推導Euler方程；一階展開f=g-τDug對應推導NS方程。這就是BGK方程求解NS方程的基礎。GKS就是通過求解BGK方程式(1)達到求解NS方程的目的。

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

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

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

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

(8)

式中：x′=xi+1/2,jm,kn-u(t-t′)是分子運動軌跡，f0是界面初始速度分布函數，g是界面上沿分子運動軌跡對應的平衡態分布函數?？梢钥吹?，這個積分解包含了時間和空間的所有信息。經過簡單計算，界面分布函數顯式表達為：

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函數，gl和gr對應界面兩側的初始平衡態分布函數，g0是界面上初始平衡態分布。g0通過如下求得。

〈aku+bkv+ckw+Ak〉=0

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

其中Q*為t*=tn+Δt/2時刻的守恒量。可以證明，由式(12)和式(13)給出的Qn+1具有四階時間精度。和四階Runge-Kutta方法相比，僅需要一個中間步就可實現四階精度，提高了數值格式的計算效率[20]。更重要的是，兩步四階格式有著和二階格式相當的數值穩定性。至此，基于二階GKS通量，五階WENO空間重構和兩步時間離散的兩步四階HGKS簡述完畢。

2 數值算例

本節中，選取經典的無黏算例、層流算例、低雷諾數湍流直接數值模擬算例和高雷諾數工程湍流RANS模擬算例，來驗證HGKS計算復雜流動問題的能力。

2.1 Titarev-Toro問題

第一個算例是Titarev-Toro 問題[31]，初值條件如下：

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

圖1 Titarev-Toro問題，t=5的密度數值解和精確解[20]

2.2 雙馬赫反射問題

雙馬赫反射問題是可壓縮無黏流動的經典算例[32]。該算例的計算域為[0,4]×[0,1]。初始時刻x=1/6處有一道與x軸呈60°的馬赫數為10的激波，激波前后的初值條件如下：

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

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

圖2 雙馬赫反射問題，t=0.2的數值解[19]

2.3 平板邊界層問題

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

圖3 平板邊界層問題，無量綱的U和V，數值結果和Blasius 參考解[19]

2.4 低速槽道湍流

槽道湍流是研究壁湍流的典型算例。低速槽道湍流計算條件：馬赫數Ma=0.1，壁面摩擦雷諾數Reτ=180。計算域為[0,2π]×[-1,1]×[0,π]，G1和G2兩套計算網格分別為963和1283。初始速度場采用泊肅葉理論解疊加幅值為當地速度10%的白噪聲隨機場。通過外加質量力驅動流動, 保證流量為常數。流向和展向給定為周期邊界條件，法向壁面為等溫壁。

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

(a) 平均流向速度

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

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

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

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

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

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

(a) 速度場脹量θ概率密度分布函數

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

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

2.6 亞聲速NACA0012翼型湍流

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

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

圖8中顯示了NACA0012翼型周圍的壓力系數分布和翼型上壁面摩擦系數分布。粗網格G3上，相比于二階隱式GKS的結果，隱式HGKS的壓力系數和摩擦系數與實驗數據[38]和CFL3D在密網格G4上參考解更加接近。阻力系數CD對翼型周圍的壓力分布和壁面摩擦力分布非常敏感。 如表1所示，隱式HGKS的升力系數CL和阻力系數CD在粗網格G3上求解結果非常接近密網格G4上CFL3D的參考解。 特別是粗網格G3上隱式HGKS的CD和參考解的差異在3個阻力數(0.0001)之內，達到了工程湍流模擬的要求。 對于二階隱式GKS，升力系數是可以接受的，而它過高地預測了阻力系數達到40個阻力數。 該算例證實了隱式HGKS高精度RANS模擬湍流場的能力。

(a)

表1 NACA0012升力系數CL和阻力系數CD

2.7 跨聲速ARA M100翼身湍流

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

圖9 ARA M100翼身構型[34]

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

(a)

3 結 論

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

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

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