李啟兵, 張 潮

(清華大學(xué) 航天航空學(xué)院, 北京 100084)

0 引 言

計(jì)算效率是CFD格式在復(fù)雜流動(dòng)模擬中面臨的核心問(wèn)題。非結(jié)構(gòu)網(wǎng)格能更好地適應(yīng)復(fù)雜外形，網(wǎng)格生成和自適應(yīng)動(dòng)態(tài)調(diào)整自動(dòng)化程度高，能極大地減少網(wǎng)格生成工作量。發(fā)展非結(jié)構(gòu)網(wǎng)格上的高效數(shù)值格式，具有重要的理論研究和工程應(yīng)用價(jià)值，也具有很大的挑戰(zhàn)性。緊致高階精度重構(gòu)和高效通量演化方法是其中的關(guān)鍵。此外，還需要考慮格式的并行效率以及穩(wěn)健性等因素。

傳統(tǒng)有限體積法由于每個(gè)單元只包含單元平均值，難以實(shí)現(xiàn)緊致高階重構(gòu)。這里的緊致是指重構(gòu)模板最多只包含面相鄰單元。緊致的計(jì)算模板不僅能更好地與可壓縮流動(dòng)的局域性相容，而且能有效避免非緊致帶來(lái)的一些弊端。例如，由于非結(jié)構(gòu)網(wǎng)格分布的任意性，距離較遠(yuǎn)的模板單元的選取較為復(fù)雜。此外，模板單元過(guò)多也會(huì)導(dǎo)致緩存利用率和并行效率降低、計(jì)算域邊界處理更加困難，而且網(wǎng)格質(zhì)量對(duì)重構(gòu)的影響較大。

內(nèi)自由度方法通過(guò)在網(wǎng)格單元內(nèi)設(shè)置足夠的自由度來(lái)實(shí)現(xiàn)緊致重構(gòu)，并利用控制方程對(duì)每個(gè)內(nèi)自由度加以記錄和更新。這里的自由度表示的是每個(gè)單元內(nèi)記錄的離散信息，如內(nèi)點(diǎn)值、子單元平均值或多項(xiàng)式展開(kāi)系數(shù)，與氣體動(dòng)理論中涉及的氣體分子的自由度是不同的概念。內(nèi)自由度方法中最具有代表性的是間斷伽遼金(DG)方法[1-2]，也是目前最為流行的非結(jié)構(gòu)網(wǎng)格高精度方法。Huynh 針對(duì)一維雙曲守恒律提出了一種簡(jiǎn)單高效的內(nèi)自由度方法，稱(chēng)為通量重構(gòu)(Flux Reconstruction, FR)格式[3-4]。FR首先計(jì)算求解點(diǎn)(內(nèi)點(diǎn))上的通量值，并插值得到單元內(nèi)的通量分布，再結(jié)合單元界面通量對(duì)其進(jìn)行修正，從而可以直接計(jì)算通量散度值來(lái)更新守恒量。相比于傳統(tǒng)的DG方法，F(xiàn)R方法具有很高的效率。LCP 方法則通過(guò)引入提升算子的方式進(jìn)行通量修正，從而使之適用于三角形和四面體網(wǎng)格[5]。FR和LCP后來(lái)統(tǒng)稱(chēng)CPR方法，為多種內(nèi)自由度方法提供了一個(gè)統(tǒng)一框架[6-7]。

相比于傳統(tǒng)有限體積法，內(nèi)自由度方法不僅易于實(shí)現(xiàn)緊致重構(gòu)，而且由于單元內(nèi)多個(gè)自由度共享同一個(gè)連續(xù)的多項(xiàng)式分布，重構(gòu)效率更高。同時(shí)，單元內(nèi)的連續(xù)分布使得計(jì)算通量時(shí)不必求解Riemann問(wèn)題，從而更為簡(jiǎn)單高效。但是，在捕捉流場(chǎng)間斷時(shí)，單元內(nèi)的連續(xù)分布難以描述間斷結(jié)構(gòu)，重構(gòu)時(shí)需要針對(duì)單元整體進(jìn)行限制，從而導(dǎo)致分辨率的降低，相當(dāng)于損失了自由度。在這方面，譜體積方法(SV)具有獨(dú)特的優(yōu)勢(shì)[8-9]。它結(jié)合了內(nèi)自由度方法和有限體積法的特點(diǎn)，在每個(gè)單元內(nèi)劃分若干個(gè)子單元，將子單元平均值作為內(nèi)自由度，按照有限體積法更新，保證了高階重構(gòu)的緊致性。由于可以分別對(duì)每個(gè)子單元進(jìn)行限制，因而能夠避免內(nèi)自由度損失，實(shí)現(xiàn)在子單元尺度上捕捉間斷[10]。已有研究對(duì)DG和SV進(jìn)行了混合[11-12]。但SV 格式對(duì)于子單元?jiǎng)澐值囊笫挚量?，?duì)格式穩(wěn)定性影響較大，不易推廣到四邊形或六面體網(wǎng)格。

利用面相鄰子單元參與重構(gòu)，Pan 等發(fā)展了針對(duì)一維及二維四邊形網(wǎng)格上的雙曲守恒方程的子單元有限體積法(Subcell Finite Volume, SCFV)。相比于SV，SCFV劃分的子單元更少，也更簡(jiǎn)單靈活[13-14]。這種通過(guò)引入面相鄰單元自由度參與重構(gòu)的方法，集合了內(nèi)自由度方法和有限體積方法的優(yōu)勢(shì)，與DG結(jié)合后也能有效降低計(jì)算量，增大時(shí)間步長(zhǎng)[15-16]。為獲得穩(wěn)定的數(shù)值格式，重構(gòu)時(shí)子單元的守恒性非常重要。SCFV通過(guò)一種簡(jiǎn)單的平移修正來(lái)間接保證守恒性。但是，這會(huì)在單元內(nèi)部的子單元界面上額外引入間斷，從而導(dǎo)致子單元界面的通量計(jì)算效率低于SV 方法。

除了高精度重構(gòu)外，通量求解器也對(duì)CFD方法的精度和效率起著重要的影響作用。傳統(tǒng)CFD方法多是基于宏觀Euler方程的Riemann解來(lái)計(jì)算數(shù)值通量。由于傳統(tǒng)的Riemann解只有一階時(shí)間精度，通常需要結(jié)合多步Runge-Kutta(R-K)方法進(jìn)行時(shí)間推進(jìn)。多步法中每個(gè)子步都需要進(jìn)行重構(gòu)、通量計(jì)算以及并行通信，影響格式的計(jì)算效率。Li和Du利用通量及其一階時(shí)間導(dǎo)數(shù)，提出了兩步四階方法，通過(guò)兩個(gè)子步達(dá)到四階時(shí)間精度[17]。基于廣義Riemann問(wèn)題(GRP)通量求解器，兩步四階方法成功用于Euler方程的高效求解[18-19]。此外，傳統(tǒng)高精度格式需要分別計(jì)算無(wú)黏和黏性通量。由于計(jì)算黏性通量需要用到單元界面上的物理量的一階導(dǎo)數(shù)，而重構(gòu)得到的兩側(cè)導(dǎo)數(shù)值通常是不連續(xù)的，從而使得計(jì)算往往較為復(fù)雜。

徐昆等提出的氣體動(dòng)理學(xué)格式(GKS)為發(fā)展高精度格式提供了一條新的途徑[20-21]。GKS利用介觀BGK方程的局部解析解來(lái)計(jì)算網(wǎng)格界面上的數(shù)值通量。該局部解在連續(xù)流區(qū)可以用平衡態(tài)分布及其時(shí)空導(dǎo)數(shù)，即宏觀量及其時(shí)空導(dǎo)數(shù)來(lái)表示，因而GKS可以直接記錄和更新宏觀守恒量，避免了在微觀速度空間離散導(dǎo)致的巨大計(jì)算資源開(kāi)銷(xiāo)。與N-S方程不同，BGK方程能描述任意重構(gòu)初始間斷的演化，因而GKS能夠引入更為合理的數(shù)值耗散，具有很強(qiáng)的穩(wěn)健性，同時(shí)也能更好地考慮多維特性[22]。由于缺乏多維Riemann解，這在傳統(tǒng)格式中是難以做到的。此外，GKS通量自動(dòng)耦合了無(wú)黏和黏性影響，且包含時(shí)間的演化，只需要單步即可達(dá)到二階時(shí)間精度[23]。二階GKS在高馬赫數(shù)流動(dòng)、多介質(zhì)流動(dòng)以及湍流等多種流動(dòng)中顯示了優(yōu)異的性能[24-27]。

為進(jìn)一步提高計(jì)算效率，可以基于時(shí)空二階Taylor展開(kāi)，利用BGK方程的局部解析解，建立時(shí)空三階精度的動(dòng)理學(xué)格式(HGKS)[28-29]。由于通量中包含了三階時(shí)空演化，具有真正的多維特性，HGKS可以通過(guò)直接積分來(lái)計(jì)算界面通量，避免了采用數(shù)值積分時(shí)多個(gè)積分點(diǎn)上的通量計(jì)算[30-31]。同時(shí)，時(shí)間方向只需要單步推進(jìn)即可達(dá)到三階精度，避免了傳統(tǒng)格式中采用多步法帶來(lái)的不足[32-33]。

將高階GKS拓展到非結(jié)構(gòu)網(wǎng)格上時(shí)，同樣面臨緊致重構(gòu)的問(wèn)題?？紤]到GKS構(gòu)造了BGK方程的局部時(shí)空演化解來(lái)計(jì)算數(shù)值通量，該演化解可以用來(lái)直接計(jì)算界面上的宏觀量，從而為重構(gòu)提供額外的信息，提高緊致性[34-36]。按照這一思路，已建立了三角形網(wǎng)格上的緊致三階HWENO-GKS[37-38]。如何更好地利用GKS提供的額外信息，進(jìn)一步提高格式精度以及穩(wěn)定性，是一個(gè)很有意思的研究方向。結(jié)合緊致最小二乘重構(gòu)，也可以發(fā)展非結(jié)構(gòu)網(wǎng)格上的高階GKS[39]，不過(guò)需要求解大型代數(shù)方程組，在顯式格式中計(jì)算量相對(duì)較大。Ren等基于DG框架發(fā)展了非結(jié)構(gòu)網(wǎng)格上的高階GKS[40]，能夠有效保證緊致性?；诟痈咝У膬?nèi)自由度方法，如CPR，有望發(fā)展更為高效的高階GKS。

在時(shí)間推進(jìn)格式方面，由于GKS通量包含時(shí)間演化，可以直接積分得到每個(gè)時(shí)間步內(nèi)的數(shù)值通量，從而建立單步高精度格式[29]。也可以借鑒兩步四階方法，利用相對(duì)更為簡(jiǎn)單二階GKS通量及其一階時(shí)間導(dǎo)數(shù)，通過(guò)兩個(gè)子步獲得四階時(shí)間精度[41]?？紤]到重構(gòu)過(guò)程通常需要較大的計(jì)算量，更少的子步意味著更小的計(jì)算量，以及更高的并行效率。同時(shí)，GKS無(wú)需額外計(jì)算黏性通量，且內(nèi)涵自適應(yīng)耗散機(jī)制，將其與高效緊致重構(gòu)方法相結(jié)合，充分挖掘其多維及時(shí)空演化特性，有望構(gòu)造出高精度、高效率和強(qiáng)穩(wěn)健性的數(shù)值格式。

本文旨在探討GKS通量求解器與高效緊致重構(gòu)方法的有機(jī)結(jié)合，發(fā)展非結(jié)構(gòu)網(wǎng)格上的高效高精度格式。首先將簡(jiǎn)要介紹課題組近年來(lái)發(fā)展的幾種基于內(nèi)自由度緊致重構(gòu)的高精度動(dòng)理學(xué)格式。將三階GKS通量求解器與SCFV和CPR相結(jié)合，構(gòu)造了兩種單步時(shí)空三階格式SCFV-GKS和CPR-GKS。進(jìn)一步將二階GKS通量與CPR結(jié)合，利用兩步四階方法構(gòu)造了時(shí)空四階格式CPR-GKS，并結(jié)合SCFV及曲邊網(wǎng)格技術(shù)增強(qiáng)對(duì)流場(chǎng)間斷和復(fù)雜外形的模擬能力。在此基礎(chǔ)上，通過(guò)多種典型流動(dòng)對(duì)這幾種格式進(jìn)行驗(yàn)證，對(duì)比分析不同格式的精度、效率以及激波捕捉性能。研究表明，結(jié)合高效的重構(gòu)方法和時(shí)空演化的GKS通量，充分利用各自?xún)?yōu)勢(shì)，能夠建立更加高效、穩(wěn)健的新型數(shù)值方法。

1 基于內(nèi)自由度的緊致高精度GKS

1.1 GKS的基本原理

GKS基于介觀氣體動(dòng)理論中的Boltzmann-BGK方程，

g=ρ(2πRT)-(K+3)/2e-[(u-U)2+ξ2]/(2RT)

(1)

其中，f=f(x,u,ξ,t)為氣體分子速度分布函數(shù)，u為分子平動(dòng)速度，ξ為內(nèi)自由度，內(nèi)自由度數(shù)為K。g為平衡態(tài)分布函數(shù)，(ρ,U,T)分別為氣體的密度、速度和溫度，τ為分子平均碰撞時(shí)間，R為氣體常數(shù)。對(duì)BGK方程取矩，可以得到宏觀守恒量Q及守恒方程：

Fα=〈uαf〉,α=1,2,3

Ψ=(1,u,(u2+ξ2)/2)T

(2)

在連續(xù)流區(qū)，通過(guò)一階Chapmann-Enskog展開(kāi)，

將分布函數(shù)近似用平衡態(tài)及其時(shí)空導(dǎo)數(shù)表示，可以導(dǎo)出宏觀的N-S方程。由此可以建立逼近N-S方程的動(dòng)理學(xué)格式，且只需要記錄和更新宏觀量(對(duì)應(yīng)平衡態(tài)分布)，避免在速度空間和內(nèi)自由度空間的離散，使得計(jì)算量和傳統(tǒng)直接基于宏觀方程的CFD方法相當(dāng)。

由方程(2)可以建立有限體積形式的GKS，

其中，Si為離散物理空間上網(wǎng)格單元i的界面面積，Ωi為單元體積。為計(jì)算網(wǎng)格單元界面數(shù)值通量，GKS利用BGK方程(1)的形式解：

e-t/τf0(x-ut,u,ξ)，

x′=x-u(t-t′)

(5)

來(lái)描述界面上任意重構(gòu)初值的演化，其中，x′=x-u(t-t′)。該演化解描述了從初始分布函數(shù)f0向平衡態(tài)分布的時(shí)空演化，自動(dòng)耦合了分子的自由運(yùn)動(dòng)和相互碰撞，這是確?；谠撗莼獾腉KS成為真正多尺度方法的關(guān)鍵[23]。從宏觀上看，不僅耦合處理了無(wú)黏與黏性項(xiàng)，而且引入了隨時(shí)間變化的自適應(yīng)數(shù)值耗散，保證了其在可壓縮黏性流動(dòng)模擬中的優(yōu)異性能。

在求解連續(xù)流區(qū)流動(dòng)時(shí)，GKS基于式(3)，結(jié)合時(shí)空Taylor展開(kāi)來(lái)構(gòu)造初始分布函數(shù)f0，并且相應(yīng)地用Taylor展開(kāi)構(gòu)造平衡態(tài)分布g，從而可以借助式(5)獲得界面上的演化解，進(jìn)而直接計(jì)算數(shù)值通量。其中Taylor展開(kāi)的階數(shù)對(duì)應(yīng)于重構(gòu)多項(xiàng)式階數(shù)，直接取決于格式精度的要求。通過(guò)二階Taylor展開(kāi)，可以得到單元界面上具有時(shí)間和空間三階截?cái)嗾`差的分布函數(shù)(以二維為例，界面坐標(biāo)x=0):

f(0,y,t,u,ξ)=flH(u)+fr[1-H(u)]+fc,

fc={C4+C5akuk+C6A+C8(A2+B)/2+

C7(akam+dkm)ukum/2+C9(Aak+bk)uk+

[C4a2+C6(Aa2+b2)+C5(aka2+dk2)uk]y+

C4(a2)2+d22)y2/2}g0

(6)

其中H為Heaviside 函數(shù)；上標(biāo)s=l,r；下標(biāo)k=1,2；m=1,2，系數(shù):

C1=τ+t,C2=τt+t2/2,C3=τt,

C4=1-e-t/τ,C5=-τ+(τ+t)e-t/τ,

C6=t-τ+τe-t/τ,C7=-(t2+2τt)e-t/τ,

C8=t2-2τt,C9=-τt(1+e-t/τ)

(7)

分布函數(shù)中的各項(xiàng)展開(kāi)系數(shù)可以由重構(gòu)出的宏觀物理量各階導(dǎo)數(shù)及相容條件給出：

〈ak〉=?Q/?xk,

〈akam+dkm〉=?2Q/?xk?xm,

〈akum+A〉=0,

〈(akam+dkm)uk+Aak+bk〉=0,

〈(Aak+bk)uk+A2+B〉=0

(8)

單元界面上的守恒量及其一、二階導(dǎo)數(shù)可由高階重構(gòu)得到。對(duì)氣體分布函數(shù)求矩即可得到二階時(shí)空分布的守恒量通量F(0,y,t)。由于其中包含了高階時(shí)空導(dǎo)數(shù)，只需要直接對(duì)時(shí)間和界面切向坐標(biāo)積分就可以得到界面數(shù)值通量，從而建立單步時(shí)空三階精度GKS，且不需要布置高斯積分點(diǎn)，計(jì)算量得到顯著減少。

需要說(shuō)明的是，式(6)中包含了界面兩側(cè)的初始間斷分布，從而有利于捕捉流場(chǎng)間斷。實(shí)際上，對(duì)于光滑區(qū)流動(dòng)，比如CPR內(nèi)點(diǎn)，式(6)可以簡(jiǎn)化為：

f(x,t,u,ξ)={1-τakuk+(t-τ)A-

τt(Aak+bk)uk+(t2/2-τt)(A2+B)+

[ak-τ(akam+dkm)um+

(t-τ)(Aak+bk)]xk+

(akam+dkm)xkxm/2}g0

(9)

因此可以獲得整個(gè)單元內(nèi)通量的時(shí)空分布F(x,t)，從而一次計(jì)算出所有內(nèi)部點(diǎn)上的通量，減少計(jì)算量。

如果略去分布函數(shù)中的高階展開(kāi)項(xiàng)，對(duì)應(yīng)于二階GKS，式(6)可簡(jiǎn)化為：

fc=(C4+C5akuk+C6A+C4a2y)g0

(10)

如果直接對(duì)時(shí)間積分，可以達(dá)到單步二階時(shí)間精度。實(shí)際上，由于該分布函數(shù)包含了一階時(shí)間導(dǎo)數(shù)項(xiàng)，可以結(jié)合兩步四階離散方法，從而高效地達(dá)到四階時(shí)間精度。通過(guò)結(jié)合新型的緊致重構(gòu)框架，可以充分發(fā)揮三階GKS和二階GKS的優(yōu)勢(shì)，構(gòu)造出更加高效的高精度格式。

1.2 SCFV-GKS

SCFV方法具有混合方法的思想，其放寬了子單元數(shù)的限制，同時(shí)將重構(gòu)模板拓展到面相鄰單元。這樣能夠有效避免SV方法劃分子單元的瓶頸。如圖1所示，為了滿(mǎn)足三階和四階緊致重構(gòu)的需求，可以將每個(gè)單元(也稱(chēng)為主單元)均勻劃分為4個(gè)子單元。子單元界面分為兩種類(lèi)型，當(dāng)界面(實(shí)線(xiàn))兩側(cè)為不同主單元時(shí)，將其稱(chēng)為I-型界面，而對(duì)于主單元內(nèi)部的界面(虛線(xiàn))則稱(chēng)為II-型界面。

圖1 三階SCFV方法的子單元?jiǎng)澐?/p>

原始的SCFV方法采用加權(quán)最小二乘法重構(gòu)單元內(nèi)的分布，由此得到的多項(xiàng)式只能保證主單元的守恒性，其內(nèi)部子單元的守恒性無(wú)法直接保證。Pan等采用了平移修正的方式對(duì)子單元守恒性進(jìn)行修正，但這樣會(huì)在子單元界面上額外引入間斷，導(dǎo)致通量計(jì)算量的增大。

我們首先通過(guò)約束最小二乘重構(gòu)對(duì)SCFV方法進(jìn)行改進(jìn)。在重構(gòu)目標(biāo)單元的多項(xiàng)式時(shí)，將目標(biāo)單元的子單元平均值作為約束條件，并通過(guò)標(biāo)準(zhǔn)的拉格朗日乘子方法求解待定系數(shù)矩陣[42]。由此得到的多項(xiàng)式能夠直接保證子單元的守恒性，而無(wú)需進(jìn)行額外的修正。在此基礎(chǔ)上，將緊致三階SCFV與三階GKS相結(jié)合，發(fā)展的了具有單步時(shí)空三階精度的SCFV-GKS。緊致三階重構(gòu)模板僅包含目標(biāo)主單元的4個(gè)子單元以及面相鄰主單元內(nèi)的9個(gè)子單元。子單元平均值按照式(4)分別更新。

為計(jì)算子單元界面上的數(shù)值通量，可以基于式(6)，直接沿子單元界面切向和時(shí)間方向積分，從而不需要空間高斯積分點(diǎn)以及時(shí)間方向的R-K推進(jìn)。值得注意的是，基于對(duì)SCFV方法的改進(jìn)，僅I-型界面上保留了間斷。而II-型界面上的分布則是連續(xù)的，因此可以基于簡(jiǎn)化的氣體分布函數(shù)式(9)進(jìn)行計(jì)算，從而有效降低通量的計(jì)算量。同時(shí)，在光滑區(qū)域，由于重構(gòu)是針對(duì)大單元進(jìn)行的，不需要對(duì)每個(gè)子單元分別進(jìn)行重構(gòu)，同樣可以減小計(jì)算量?？傊?，SCFV-GKS不僅能保證緊致性，計(jì)算效率也能得到顯著提高。

為了實(shí)現(xiàn)激波捕捉，可以通過(guò)激波探測(cè)器標(biāo)記問(wèn)題單元，并在問(wèn)題單元上實(shí)施基于逐級(jí)重構(gòu)的WENO限制器[43]，該限制器能夠在保證光滑區(qū)域高精度的同時(shí)，有效抑制激波附近的數(shù)值振蕩。由于限制后的多項(xiàng)式無(wú)法直接保證子單元的守恒性，為簡(jiǎn)單起見(jiàn)，仍采用原始SCFV方法中的平移修正來(lái)間接保證子單元的守恒性，實(shí)現(xiàn)子單元尺度上的間斷捕捉。SCFV-GKS的具體細(xì)節(jié)見(jiàn)文獻(xiàn)[44]。

1.3 基于CPR框架的高精度GKS

CPR框架的表達(dá)式如下：

(11)

其利用單元內(nèi)設(shè)置的若干求解點(diǎn)可以插值得到守恒量和通量在單元內(nèi)的分布，從而可以直接對(duì)通量函數(shù)多項(xiàng)式求散度，進(jìn)而更新每個(gè)求解點(diǎn)上的守恒變量Qij。考慮到單元界面兩側(cè)可能存在間斷，還需在每個(gè)單元界面上設(shè)置若干通量點(diǎn)，計(jì)算耦合相鄰單元影響的通量值，并由此對(duì)單元內(nèi)的通量分布進(jìn)行修正，即式(11)中的修正項(xiàng)δij。三階CPR的求解點(diǎn)和通量點(diǎn)分布如圖2所示，為減少計(jì)算量，求解點(diǎn)與通量點(diǎn)均采用Lobatto點(diǎn)，兩者位置重合。

圖2 三階CPR的求解點(diǎn)(圓形)和通量點(diǎn)分布(方形)

我們首先基于CPR的緊致三階重構(gòu)，結(jié)合三階GKS通量求解器，發(fā)展了單步時(shí)空三階精度的CPR-GKS。充分利用三階GKS的多維特性，基于單元界面中點(diǎn)，以及單元形心分別構(gòu)造了時(shí)空二階分布的氣體分布函數(shù)，即式(6)和式(9)，進(jìn)而得到通量在單元界面以及單元內(nèi)的時(shí)空二階分布。由此可以同時(shí)確定單元界面上所有通量點(diǎn)，以及單元內(nèi)所有求解點(diǎn)上的通量值。傳統(tǒng)CPR則需要逐點(diǎn)計(jì)算數(shù)值通量。同時(shí)，由于通量包含了時(shí)間演化，因此可以沿時(shí)間積分單步更新求解點(diǎn)值，而不需要結(jié)合多步R-K方法。此外，由于GKS中自動(dòng)耦合了無(wú)黏和黏性影響，毋須像傳統(tǒng)CPR那樣額外處理黏性通量。格式的具體細(xì)節(jié)見(jiàn)文獻(xiàn)[45]。

為進(jìn)一步提高格式精度，我們基于CPR的緊致四階重構(gòu)，結(jié)合二階GKS通量及其時(shí)間導(dǎo)數(shù)，借助兩步四階離散，發(fā)展了兩步四階GKS，只需要兩個(gè)子步即可達(dá)到時(shí)空四階精度。而傳統(tǒng)CPR則通常需要結(jié)合四步或五步R-K方法。相比于單步四階GKS，采用兩步四階GKS能夠有效簡(jiǎn)化氣體分布函數(shù)的構(gòu)造，降低通量的計(jì)算量。

兩步四階CPR-GKS的高效性，來(lái)源于對(duì)GKS通量中包含的豐富時(shí)空信息的充分利用。類(lèi)似地，單步三階CPR-GKS利用了GKS通量的多維時(shí)空特性，從而顯著提高了計(jì)算效率。Xu等利用界面上的演化解直接計(jì)算下一時(shí)刻的宏觀量，從而提高重構(gòu)緊致性[34-38]。實(shí)際上，如果充分利用四階GKS通量的時(shí)空信息，也有望建立更為高效的數(shù)值格式。

雖然基于CPR框架能夠高效地實(shí)現(xiàn)緊致高階重構(gòu)，但在高速流動(dòng)問(wèn)題中的激波捕捉技術(shù)尚待進(jìn)一步研究。通過(guò)借鑒SCFV方法在激波捕捉上的優(yōu)勢(shì)，我們進(jìn)一步發(fā)展了CPR/SCFV混合算法。通過(guò)激波探測(cè)器區(qū)分流場(chǎng)中的光滑區(qū)域和包含間斷的區(qū)域，前者采用CPR求解，而后者則采用SCFV求解。這樣既能夠保證光滑區(qū)域的高精度，同時(shí)兼顧激波捕捉的高分辨率和強(qiáng)穩(wěn)健性。傳統(tǒng)CPR限制器通常直接針對(duì)單元內(nèi)的高階多項(xiàng)式進(jìn)行限制，難以在亞網(wǎng)格尺度上捕捉間斷。而在CPR/SCFV混合算法中，通過(guò)將激波探測(cè)器標(biāo)記的問(wèn)題單元?jiǎng)澐譃橐唤M子單元，并且針對(duì)每個(gè)子單元分別進(jìn)行限制。如圖3所示，對(duì)于三階CPR，將問(wèn)題單元均勻劃分為9個(gè)子單元，這樣能夠在保證原有內(nèi)自由度的同時(shí)，避免子單元?jiǎng)澐值膹?fù)雜性。類(lèi)似的，四階CPR則將問(wèn)題單元均勻劃分為16個(gè)子單元。具體細(xì)節(jié)見(jiàn)文獻(xiàn)[46]。

圖3 三階CPR的子單元?jiǎng)澐?/p>

通過(guò)在子單元界面上引入間斷，激波的厚度能夠被有效地限制在大單元尺度范圍內(nèi)，從而實(shí)現(xiàn)子單元尺度上的激波捕捉。需要指出的是，我們?cè)诿總€(gè)子單元上采用的是二階TVD限制器[10]，該限制器能夠較好地抑制數(shù)值振蕩，并且實(shí)現(xiàn)較為簡(jiǎn)單，但引入的數(shù)值耗散較大，也在一定程度上增加了對(duì)激波探測(cè)器的依賴(lài)性。在后續(xù)的研究中，可以考慮針對(duì)子單元實(shí)施高階限制器。

高階格式在處理曲面邊界時(shí)，網(wǎng)格對(duì)物理邊界的逼近程度會(huì)嚴(yán)重影響數(shù)值計(jì)算精度，有必要采用曲邊網(wǎng)格?；贑PR 框架發(fā)展的高精度GKS，能夠比較方便地處理曲邊網(wǎng)格。只需通過(guò)等參變換將物理空間的高階網(wǎng)格單元變換到計(jì)算空間的標(biāo)準(zhǔn)單元上，并在標(biāo)準(zhǔn)單元實(shí)施CPR 算法即可。同時(shí)，由一維CPR 方法通過(guò)張量積的形式將其推廣到六面體網(wǎng)格，從而建立適合復(fù)雜外形的三維CPR-GKS。

2 數(shù)值驗(yàn)證及對(duì)比分析

通過(guò)典型算例對(duì)前面提到的幾種格式進(jìn)行驗(yàn)證，考察不同格式的精度、效率以及激波捕捉性能。為了說(shuō)明與傳統(tǒng)通量求解器的區(qū)別，還與傳統(tǒng)CPR進(jìn)行了一定的對(duì)比。二維問(wèn)題中直接基于單元尺度計(jì)算的CFL數(shù)，SCFV-GKS取為0.2，CRP-GKS以及傳統(tǒng)CPR均取0.1，三維問(wèn)題中CFL數(shù)取值為0.05。

2.1 可壓縮Couette流動(dòng)

考察在兩無(wú)限長(zhǎng)平行平板之間的流動(dòng)，上板溫度為T(mén)1= 1，運(yùn)動(dòng)速度為U1= 0.5，相應(yīng)的馬赫數(shù)為Ma=0.5，下板靜止且絕熱。雷諾數(shù)設(shè)為Re=500。黏性系數(shù)μ與溫度呈線(xiàn)性關(guān)系μ=μ1T/T1，普朗特?cái)?shù)為Pr=1。為了簡(jiǎn)單起見(jiàn)，邊界條件按照解析解給定。需要說(shuō)明的是，雖然流動(dòng)定常，這里仍按非定常顯式推進(jìn)，計(jì)算至足夠長(zhǎng)的時(shí)間(t=300)后評(píng)估各格式所得結(jié)果的誤差與所用CPU時(shí)間。圖4給出了速度的2-范數(shù)誤差隨網(wǎng)格尺度變化的曲線(xiàn)。其中傳統(tǒng)CPR的無(wú)黏通量采用HLLC格式，黏性通量采用BR2格式，并且結(jié)合四步R-K方法進(jìn)行時(shí)間推進(jìn)?？梢钥闯鏊懈袷蕉歼_(dá)到了設(shè)計(jì)的精度階數(shù)，絕對(duì)誤差主要受重構(gòu)階數(shù)影響。相比于三階格式，四階格式具有顯著更小的絕對(duì)誤差，并且隨著網(wǎng)格加密，絕對(duì)誤差下降得更快。

圖4 可壓縮Couette流動(dòng)：速度誤差隨網(wǎng)格尺度的變化

為進(jìn)一步考察格式的計(jì)算效率，圖5給出了絕對(duì)誤差隨計(jì)算所用CPU時(shí)間變化的曲線(xiàn)。雖然在相同網(wǎng)格尺度下，四階格式的計(jì)算量大于三階格式，但由于絕對(duì)誤差明顯更小，因此四階格式相比于三階格式具有明顯的效率優(yōu)勢(shì)。對(duì)于傳統(tǒng)CPR，雖然傳統(tǒng)Riemann求解器的計(jì)算量較小，但在黏性問(wèn)題中還需要額外處理黏性項(xiàng)，并且由于需要更多的子步，因此總的計(jì)算量與CPR-GKS相當(dāng)。總之，基于CPR框架發(fā)展的兩步四階GKS具有最顯著的優(yōu)勢(shì)。需要指出的是，本算例為光滑流動(dòng)問(wèn)題，因此沒(méi)有施加限制器。對(duì)于高速流動(dòng)問(wèn)題，施加限制器往往帶來(lái)較大的計(jì)算量，由于CPR-GKS能夠采用更少的子步數(shù)達(dá)到高階時(shí)間精度，因此其計(jì)算效率相比于傳統(tǒng)CPR還能夠進(jìn)一步提高。此外，考慮到三階CPR-GKS可以比四階格式取更大的CFL數(shù)，此時(shí)效率會(huì)有一定的提升，但仍然顯著低于四階格式。不同格式CFL數(shù)取值大小的影響仍有待進(jìn)一步研究。

圖5 可壓縮Couette流動(dòng)：速度誤差隨CPU時(shí)間的變化

2.2 雙馬赫反射

馬赫數(shù)Ma=10的正激波經(jīng)過(guò)30°斜坡形成的雙馬赫反射是典型的二維強(qiáng)激波問(wèn)題，可用于驗(yàn)證格式的穩(wěn)健性和激波捕捉能力。初始流場(chǎng)為:

計(jì)算時(shí)間為t=0.2，網(wǎng)格尺度為h=1/160。圖6給出了不同格式計(jì)算得到的三波點(diǎn)附近的密度等值線(xiàn)圖，其中包含30條在2.0～22.5之間均勻分布的密度等值線(xiàn)。可以看出，兩步四階CPR-GKS的計(jì)算結(jié)果最好，激波的分辨率最高，并且在滑移線(xiàn)上捕捉到了一系列小尺度渦結(jié)構(gòu)。相比于單步三階CPR-GKS，單步三階SCFV-GKS在滑移線(xiàn)上捕捉到的渦結(jié)構(gòu)更加明顯，但間斷略顯更寬。需要注意的是，與文獻(xiàn)[45]不同，本文的三階CPR-GKS由于在間斷區(qū)域結(jié)合了SCFV，每個(gè)單元內(nèi)包含了更多的內(nèi)自由度，因此對(duì)間斷的分辨率比SCFV-GKS更高。但是，兩者在子單元上采用的限制器是不同的：CPR-GKS采用的是二階TVD限制器，而SCFV-GKS則是HR-WENO，因此在滑移線(xiàn)附近，由于部分單元被標(biāo)記為間斷單元，高階限制器使得SCFV-GKS得到的渦結(jié)構(gòu)更加明顯。計(jì)算結(jié)果表明三種高精度GKS都能夠在激波捕捉中兼具高分辨率和強(qiáng)穩(wěn)健性，并且兩步四階CPR-GKS的優(yōu)勢(shì)更加明顯。

(a) 單步三階SCFV-GKS

圖7給出了兩步四階CPR-GKS計(jì)算得到的密度等值線(xiàn)的全局視圖，斜激波下方的等值線(xiàn)也比較光滑，說(shuō)明格式能夠有效抑制數(shù)值振蕩。圖8進(jìn)一步給出了由兩步四階CPR-GKS計(jì)算得到的馬赫桿附近的密度等值線(xiàn)，且同時(shí)展示了主單元和子單元的分布。可以看出激波厚度被成功限制在主單元的尺度內(nèi)，因此能夠?qū)崿F(xiàn)子單元尺度上的激波捕捉，驗(yàn)證了CPR/SCFV混合算法對(duì)間斷的高分辨率。

圖7 雙馬赫反射：密度等值線(xiàn)的全局視圖(t=0.2)，兩步四階CPR-GKS

圖8 雙馬赫反射：馬赫桿附近的等值線(xiàn)(t=0.2)，兩步四階CPR-GKS

2.3 三維超聲速圓球繞流

超聲速來(lái)流繞過(guò)圓球時(shí)，流場(chǎng)中包含復(fù)雜的波系和旋渦結(jié)構(gòu)，對(duì)于數(shù)值方法具有較大的挑戰(zhàn)性。來(lái)流馬赫數(shù)取為Ma=1.2，雷諾數(shù)為Re=1000。計(jì)算網(wǎng)格如圖9所示，共88 128個(gè)單元，圓球表面附近的法向最小網(wǎng)格尺度為0.0365，拉伸率約為1.15，圓球表面為無(wú)滑移絕熱壁。

圖9 三維超聲速圓球繞流：計(jì)算網(wǎng)格

圖10給出了t=400的瞬時(shí)密度分布云圖，其中可以看到圓球前面的弓形激波、后部的膨脹波及流動(dòng)分離形成的回流區(qū)，以及回流區(qū)末端的再壓縮波。在尾跡區(qū)則可以觀察到有大尺度的渦脫落。圖11給出了t=400 的瞬時(shí)Q判據(jù)等值面圖，其中Q=0.001，并以馬赫數(shù)著色。從圖中可以看到，四階CPR-GKS清晰地捕捉到了尾跡區(qū)的渦結(jié)構(gòu)，且流場(chǎng)存在明顯的對(duì)稱(chēng)面。圖12給出了阻力系數(shù)隨時(shí)間變化的曲線(xiàn)，可以看出流動(dòng)結(jié)構(gòu)的演化過(guò)程仍然具有明顯的周期性，但此時(shí)已出現(xiàn)較小的隨機(jī)擾動(dòng)。由阻力系數(shù)曲線(xiàn)得到的頻譜如圖13所示，可以觀察到兩個(gè)主導(dǎo)的頻率分別為0.039 和0.145，此外還存在幾個(gè)更小的峰值頻率，并且都集中在低頻區(qū)域。統(tǒng)計(jì)得到的阻力時(shí)均值為1.055，計(jì)算結(jié)果與Nagata等給出的結(jié)果1.054吻合得較好[47]。

圖10 三維超聲速圓球繞流：密度分布云圖(t=400)

圖11 三維超聲速圓球繞流：Q判據(jù)等值面圖(t=400)

圖12 三維超聲速圓球繞流：阻力系數(shù)隨時(shí)間變化的曲線(xiàn)

圖13 三維超聲速圓球繞流：阻力頻譜

3 結(jié) 論

本文介紹了近年來(lái)課題組在非結(jié)構(gòu)網(wǎng)格高效高精度氣體動(dòng)理學(xué)格式上的研究進(jìn)展。利用介觀GKS通量求解器，結(jié)合基于內(nèi)自由度的高效緊致重構(gòu)方法，發(fā)展了單步時(shí)空三階精度的SCFV-GKS和CPR-GKS，以及時(shí)空四階精度的CPR-GKS。利用GKS通量的多維空間分布，結(jié)合單元內(nèi)部的連續(xù)分布，減少通量的計(jì)算量，同時(shí)結(jié)合通量的時(shí)間演化特性減少或消除時(shí)間方向的推進(jìn)子步，從而有效提高了計(jì)算效率。混合光滑區(qū)CPR和間斷區(qū)域SCFV發(fā)展的兩步四階CPR-GKS，很好地兼顧了高精度、高分辨率和強(qiáng)穩(wěn)健性。進(jìn)一步結(jié)合曲邊網(wǎng)格技術(shù)發(fā)展了六面體網(wǎng)格上的四階CPR-GKS，有效增強(qiáng)了格式在復(fù)雜流動(dòng)中的實(shí)用性能。

在此基礎(chǔ)上，通過(guò)典型算例對(duì)幾種格式進(jìn)行驗(yàn)證和對(duì)比分析，考察了不同格式的精度、效率以及激波捕捉性能。結(jié)果表明，其中結(jié)合SCFV的兩步四階CPR-GKS具有顯著的優(yōu)勢(shì)?；跉怏w動(dòng)理論的GKS求解器包含豐富的時(shí)空演化信息，充分利用這些信息，結(jié)合高效的緊致重構(gòu)方法，能夠發(fā)展兼具高精度、高效率和強(qiáng)穩(wěn)健性的新型數(shù)值方法，為新一代CFD軟件開(kāi)發(fā)提供有力支撐。

在后續(xù)的研究工作中，還需要對(duì)高階GKS在三維流動(dòng)問(wèn)題，特別是包含強(qiáng)激波的流動(dòng)問(wèn)題中的計(jì)算性能做進(jìn)一步的驗(yàn)證。同時(shí)，本文涉及到的均為顯式格式，為增強(qiáng)對(duì)實(shí)際工程問(wèn)題的模擬能力，有必要開(kāi)展隱式高精度格式研究。此外，為模擬工程實(shí)際中面臨的高雷諾數(shù)湍流問(wèn)題，可以基于拓展BGK方程的跨尺度演化解發(fā)展新型的可壓縮湍流多尺度模擬方法，并應(yīng)用于典型高雷諾數(shù)工程湍流精細(xì)模擬。

致謝:感謝清華探索100高性能計(jì)算平臺(tái)、并行科技高性能計(jì)算支持服務(wù)為本文工作提供部分計(jì)算機(jī)時(shí)。