尺度自適應的離散統(tǒng)一氣體動理學格式及在可壓縮流動中的應用

許 丁, 孫 祥, 劉 欣

(1. 西安交通大學 航天航空學院, 機械結構強度與振動國家重點實驗室, 西安 710049; 2. 陜西省先進飛行器服役環(huán)境與控制重點實驗室, 西安 710049)

0 引 言

近來在計算流體力學(Computational Fluid Dynamics, CFD)中，數(shù)值模擬寬范圍努森數(shù)Kn的多尺度流動問題成為一個熱點。隨著Kn的變化，流動對應領域及流動特征尺度不盡相同，對傳統(tǒng)CFD方法提出了挑戰(zhàn)。通常根據(jù)Kn對氣體流動領域進行了如下劃分[1]：在Kn<0.001時屬于連續(xù)流動領域，傳統(tǒng)基于連續(xù)介質假設的N-S方程和黏性壁面處的無滑移邊界條件是適用的，流動對應宏觀動力學特征尺度，如繞流物體特征長度，這正是傳統(tǒng)CFD方法蓬勃發(fā)展并得到廣泛應用的領域；當0.00110時流動屬于自由分子流領域，非平衡效應很強，流動對應微觀動理學特征尺度，即分子平均自由程和平均碰撞松弛時間。上述對流動領域的劃分屬于比較粗線條的，一個全場定義的Kn不足以精確刻畫流場當?shù)孛奎c的流動特征。最近陳杰等[2]提出一個反映氣體局部稀薄效應判據(jù)的工作值得關注。

盡管對于Kn較大的過渡流及自由分子流，Bird[3]提出的Direct Simulation Monte Carlo(DSMC)方法是一種有效方法，但是對于連續(xù)流動，DSMC所消耗的計算時間及計算資源將異常巨大。因此發(fā)展適用于寬范圍Kn數(shù)變化、跨流域的統(tǒng)一算法對解決實際工程當中的多尺度流動是一個關鍵問題。

另一方面，直接求解分布函數(shù)滿足的Boltzmann方程是氣體動理學方法的研究思路，近來一些氣體動理學格式也應運而生。如格子Boltzmann方法(Lattice Boltzmann Method, LBM)[4-6]、氣體動理學統(tǒng)一算法(Gas Kinetic Unified Algorithms, GKUA)[7-11]、氣體動理學格式(Gas Kinetic Scheme, GKS)[12-16]、統(tǒng)一氣體動理學格式(Unified Gas Kinetic Scheme, UGKS)[17-22]、離散統(tǒng)一氣體動理學格式(Discrete Unified Gas Kinetic Scheme, DUGKS)[23-25]等。這些方法與傳統(tǒng)基于N-S方程的CFD方法相比，可以獲得更多的流動信息，且已有大量成功的應用[26-33]。

本文從DUGKS格式出發(fā)，指出DUGKS在處理碰撞項時，采用顯式和隱式的簡單算術平均并未充分考慮不同流動領域碰撞項對應物理過程的不同，尤其沒有充分考慮兩個時間尺度即計算時間步長和碰撞松弛時間之間的相對大小對粒子碰撞過程宏觀表現(xiàn)的影響。為此，本文將對碰撞項的離散采用顯式和隱式的加權平均，且權系數(shù)緊密依賴于計算時間步長和當?shù)嘏鲎菜沙跁r間之比，即當?shù)嘏瓟?shù)。進一步通過物理上的分析與數(shù)學上的推演，表明該權系數(shù)可以根據(jù)當?shù)亓鲃犹卣鞒叨冗M行自適應調節(jié)，從而將原始DUGKS格式進一步發(fā)展成為具有尺度自適應特性的離散統(tǒng)一氣體動理學格式(Scale Adaptive Discrete Unified Gas Kinetic Scheme, SADUGKS)。

1 Boltzmann-Shakhov模型方程及其沿特征線離散的一般形式

為了克服Boltzmann-BGK模型對應普朗特數(shù)Pr≡1與實際氣體不相符，本文基于Boltzmann- Shakhov模型來構造計算格式。

在D維空間中，Boltzmann-Shakhov模型[24, 34]如下：

(1)

其中，f=f(xi,t,uj,ηk,ξl)為D維空間中t時刻、位置x=(x1,...,xD)處、速度為u=(u1,…,uD)的粒子分布函數(shù)，η=(η1,...,η3-D)為宏觀速度恒為零的方向上的粒子運動速度，ξ=(ξ1,…,ξK)為粒子內自由度，K為粒子內自由度的維度，與粒子種類相關；τ為碰撞松弛時間，與動力學黏性系數(shù)μ及壓力p相關，τ=μ/p；fs為Shakhov平衡態(tài)分布函數(shù)，可表示為Maxwell平衡態(tài)分布函數(shù)及Pr數(shù)修正兩部分。

為節(jié)約內存開銷和計算量，將f對內自由度和宏觀速度恒為零的速度方向上積分得到兩個約化分布函數(shù)：

(2)

(3)

宏觀守恒量Q=(ρ,ρU,ρE)T、應力張量τ及熱流q可由分布函數(shù)求矩得到：

(4)

(5)

(6)

其中c=u-U，geq為Maxwell平衡態(tài)分布：

(7)

對式(1)積分可得到約化分布函數(shù)g及h的控制方程：

(8)

(9)

其中，gs和hs分別為g和h的Shakhov平衡態(tài)分布函數(shù)，二者均可表示為Maxwell平衡態(tài)及Pr修正兩部分：

(10)

(11)

(12)

heq=(K+3-D)RTgeq

(13)

(14)

根據(jù)粒子碰撞過程滿足質量、動量與能量守恒定律，約化分布函數(shù)g及h滿足相應相容關系：

(15)

式(8)與式(9)形式上一致，可統(tǒng)一表示為：

(16)

其中φ=g或h。對于式(16)，t+s時刻的分布函數(shù)可形式上表示為[12]：

φ(xi,t+s,uj)=

e-s/τφini(xi-uis,uj)

(17)

其中：t為初始時刻，s為時間發(fā)展長度，x′i=xi-ui(s-t′)為粒子運動軌跡；φini為初始t時刻的氣體分布函數(shù)，即φini(xi,uj)=φ(xi,t,uj)。

盡管Boltzmann-Shakhov方程(16)基于微觀觀點建立，其所蘊含的流動尺度卻不限于微觀動理學尺度，這一點可以從式(17)清楚看出：式(17)右邊第一項積分項反映粒子之間相互大量碰撞作用的積累效果，該項已被證明可以看作為宏觀連續(xù)流動動力學模型[12, 35]，對應尺度為宏觀動力學尺度，如繞流物體特征長度；式(17)右邊第二項反映粒子的遷移運動，對應尺度為微觀動理學尺度，即分子平均自由程和平均碰撞松弛時間。式(17)將宏觀動力學尺度和微觀動理學尺度有機地結合到一起，可以很好反映多尺度流動特性。但是式(17)并不能直接使用，因為右邊第一項積分項所涉及的φs是未知的，需要滿足相容關系式(15)。因此從這個角度上來說式(17)只是式(16)形式上的解。

另一方面，認識到Boltzmann-Shakhov方程(16)自身已將宏觀動力學尺度和微觀動理學尺度有機地結合到一起，因此可以直接從Boltzmann-Shakhov方程(16)出發(fā)來構造多尺度計算格式。這里將式(16)沿特征線進行離散，可得到：

(18)

其中：

φs,ini(xi,uj)=φs(xi,t,uj)

(19)

分別代表初始時刻t的平衡態(tài)分布函數(shù)和碰撞項。式(18)右端的碰撞項采用了顯式和隱式加權平均處理。α為隱式部分對應的權系數(shù)，其為碰撞松弛時間τ和時間發(fā)展長度s的比值τ/s的函數(shù)，具體形式待定。式(18)可進一步整理為：

φ(xi,t+s,uj)=?1φs(xi,t+s,uj)+

?2φs,ini(xi-uis,uj)+?3φini(xi-uis,uj)

(20)

其中，系數(shù)?1、?2、?3為：

(21)

式(20)為Boltzmann-Shakhov方程沿特征線離散的一般形式，其右端前兩項和平衡態(tài)分布函數(shù)相關反映粒子之間的碰撞作用，將和式(17)右端第一項積分項對應；式(20)右端第三項和初始非平衡態(tài)分布函數(shù)相關反映粒子的遷移效應，將和式(17)右端第二項相對應。通過對比式(17)和式(20)中的非平衡態(tài)分布函數(shù)相關項，并令?3=e-s/τ，可確定權系數(shù)α為：

(22)

式(22)表明權系數(shù)α數(shù)學上直接依賴于時間比值τ/s，物理上的討論將在后面第3節(jié)給出。確定了α后將其帶入到式(21)，可得系數(shù)?1、?2和?3具體形式：

(23)

為了進一步分析式(20)的數(shù)學物理特性，這里將從兩種極限流動的角度來討論，首先將式(20)整理為以下形式：

φ(xi,t+s,uj)=φs(xi,t+s,uj)+

β1[φini(xi-uis,uj)-φs,ini(xi-uis,uj)]+

β2[φs,ini(xi-uis,uj)-φs(xi,t+s,uj)]

(24)

其中：

β1=e-s/τ,β2=(1-e-s/τ)τ/s

(25)

對于Kn→0的宏觀連續(xù)流動問題，物理上在一個時間步長內粒子經(jīng)歷的碰撞次數(shù)足夠多，即τ?s，數(shù)學上β1→0，β2→τ/s，式(24)可化為：

φ(xi,t+s,uj)≈φs(xi,t+s,uj)+

≈φs(xi,t+s,uj)-τDtφs(xi,t+s,uj)

(26)

其中Dt=?t+uj?xj。對式(8)進行Chapman-Enskog展開至O(τ)(對應為宏觀N-S方程)同樣可以得到式(26)[12]，因此說明在宏觀連續(xù)流動區(qū)域式(20)與N-S方程一致, 此時流動尺度為宏觀動力學尺度。

對于Kn?0的稀薄流動或微尺度流動問題，物理上τ?s，從而β1→1，β2→1，式(24)化為：

φ(xi,t+s,uj)→φini(xi-uis,uj)

(27)

上式其實就是自由分子流運動的解，此時流動對應微觀動理學尺度。

由上述分析可知，Boltzmann-Shakhov方程沿特征線離散的一般形式(20)將宏觀動力學尺度和微觀動理學尺度有機地結合到一起，其中引入的系數(shù)?1、?2和?3依賴于比值τ/s且可根據(jù)當?shù)亓鲃宇I域和流動尺度的不同進行自適應調節(jié)，比值τ/s在調節(jié)中起到了關鍵作用。

2 尺度自適應的離散統(tǒng)一氣體動理學格式

尺度自適應的離散統(tǒng)一氣體動理學格式(Scale Adaptive Discrete Unified Gas Kinetic Scheme, SADUGKS)可視為在原始DUGKS格式[24]的基礎上進行改進，二者主要區(qū)別在于權系數(shù)α的形式不同。

基于有限體積法框架，將式(16)在網(wǎng)格單元Vj及時間步長(tn,tn+Δt)內積分，得：

(28)

(29)

對式(28)中通量的時間積分項采用中點公式處理：

(30)

對式(28)中碰撞項的時間積分采用加權梯形公式：

(31)

其中α(Δt)為式(22)給出的權系數(shù):

(32)

需要指出的是在原始DUGKS格式中[24]，式(31)中的權系數(shù)取為αDUGKS≡0.5，即代表碰撞項的顯式和隱式算術平均。

將式(30)和式(31)代入到式(28)中，可得：

(33)

從式(33)可以看出該式右端包含了通量和碰撞項的隱式部分，因此在使用該式進行迭代遞推之前必須先解決通量和碰撞項的隱式離散問題，下面對其分別進行處理。

(34)

(35)

上式也等價于:

(36)

其中系數(shù)φ1和φ2為，

(37)

則式(34)可化為：

(38)

(39)

(40)

由式(4)、(15)和(35)可知:

(41)

結合式(38)位于界面xi,b處，tn+r時刻的宏觀守恒量為:

Q(xi,b,tn+r)

(42)

(43)

由相容關系式(15)同時結合定義式(43)可知:

(44)

將式(43)代入到式(33)，整理得：

(45)

再次指出的是在原始DUGKS格式中，從式(34)至式(39)，式(43)和式(45)中出現(xiàn)的所有權系數(shù)α恒取為αDUGKS≡0.5。

需要說明的是上述在推導SADUGKS中速度空間是連續(xù)的，實際在使用SADUGKS時還需將速度空間離散，其離散方法與現(xiàn)有DUGKS[23-24]格式相同。關于邊界條件的處理分兩類，針對宏觀連續(xù)流動問題，SADUGKS格式在固體壁面處采用了無滑移邊界條件；而對較大Kn流動問題，則采用了完全漫反射邊界條件，具體實現(xiàn)與DUGKS格式[23-24]亦相同，不再贅述。

3 尺度自適應的離散統(tǒng)一氣體動理學格式的性質

首先，作為DUGKS格式的一種改進，SADUGKS格式同樣具有漸近保持性[23-24]，這一點可以清楚地從式(26)看出。

其次，SADUGKS格式具有尺度自適應特性，能夠根據(jù)當?shù)亓鲃犹卣鞒叨鹊牟煌M行自適應調節(jié)。SADUGKS格式中，時間步長按照CFL條件確定：

(46)

其中，Δ為空間網(wǎng)格尺度，CFL為CFL數(shù)；|u|max為粒子最大離散速度，通常與聲速a同量級。同時，碰撞松弛時間τ可近似為：

(47)

τ/Δt～ls/Δ～Knlocal

(48)

其中Knlocal為當?shù)嘏瓟?shù)，反映當?shù)氐牧鲃訝顟B(tài)和流動尺度。在SADUGKS格式中，用來計算界面通量的式(40)和碰撞項加權離散的式(31)中皆出現(xiàn)權系數(shù)α，式(32)給出權系數(shù)α的定義式反映出α依賴于比值τ/Δt，即依賴于當?shù)嘏瓟?shù)，因此α取值的大小會隨著當?shù)亓鲃訝顟B(tài)和尺度的不同而進行自適應調節(jié)。

α→1,τ?Δt

α→1/2,τ?Δt

(49)

≈Ωφ(xi,t+Δt,uj)

(50)

可見數(shù)學結果與上述物理分析一致。

(51)

可見此時數(shù)學結果也與上述物理分析相一致。

在原始DUGKS格式中，本文所有出現(xiàn)權系數(shù)α的地方都取為常數(shù)αDUGKS≡1/2。通過上面的分析可以看出，這種常數(shù)取法主要是從數(shù)值離散的角度考慮，并未充分考慮在不同流動領域碰撞項對應物理過程的不同，尤其沒有考慮兩個時間尺度即計算時間步長和碰撞松弛時間之間的相對大小對粒子碰撞過程宏觀表現(xiàn)的影響。

根據(jù)上述分析可知，式(32)給出的權系數(shù)α緊密依賴于時間比值τ/Δt，并可根據(jù)當?shù)亓鲃訝顟B(tài)和流動尺度的不同進行自適應調節(jié)。當?shù)嘏瓟?shù)越小、越接近連續(xù)流動時α越接近于1；當?shù)嘏瓟?shù)越大、流動非平衡效應越強時α越接近于1/2。綜合來說，SADUGKS是一種尺度自適應的格式，適用于從連續(xù)流動到滑移流動、過渡流動及自由分子流所有流動領域的問題，在第4節(jié)將通過若干典型算例對這一點進行檢驗。

4 尺度自適應的離散統(tǒng)一氣體動理學格式在可壓縮流動中的應用

4.1 一維激波結構問題

本問題中將用SADUGKS格式研究一維激波結構問題，并將結果與文獻[36,37]的結果進行對比。氣體為氬氣，分子內自由度設置為K=0，普朗特數(shù)Pr=2/3，比熱比γ=5/3，黏性系數(shù)與溫度相關：μ∝Tw，其中w與分子模型相關。分子平均自由程λ與黏性系數(shù)直接相關[24]：

(52)

選取計算區(qū)域為-25λ1≤x≤25λ1，并采用100個均勻網(wǎng)格將其離散，每個網(wǎng)格單元大小Δx=0.5λ1，由于網(wǎng)格單元Δx和分子平均自由程λ1量級相當，計算將給出激波內部結構中物理量的變化。速度空間采用32×32半平面Gauss-Hermite積分離散[38]。初始時刻，x≤0處的流場采用激波上游條件及其Maxwell平衡態(tài)分布進行初始化，x>0處的流場采用激波下游條件及其Maxwell平衡態(tài)分布進行初始化。左右邊界均采用外推邊界條件,CFL數(shù)設置為0.5。

本文采用硬球模型w=0.5，圖1展示了Ma=1.2數(shù)值模擬結果，并與文獻[36,37]中的結果進行了對比。其中物理量的無量綱化方法如下：

(53)

另外，參照文獻[36]給出固定激波所在位置的方法，選取ρ(x0)=(ρ1+ρ2)/2處為坐標原點，將坐標進行無量綱化：

(54)

從圖1中可以看出，SADUGKS的結果與UGKS[37]及完整Boltzmann方程[36]都能很好地吻合。由于本問題是分辨激波內部結構，流動對應尺度為分子平均自由程尺度，按照前面關于權系數(shù)α的討論，此時α應接近0.5，從圖1中也可以清楚看出全場α范圍在0.503至0.505之間。

圖1 SADUGKS格式求解激波結構問題的密度、溫度、熱流、應力與權系數(shù)α分布(Ma=1.2)Fig.1 Distribution of density, temperature, heat flux, stress and α for shock structure problem (Ma=1.2)

4.2 從激波結構分辨格式到激波捕捉格式

在4.1小節(jié)問題中，由于所用網(wǎng)格單元Δx和分子平均自由程λ1量級相當，因此計算結果可以分辨激波內部結構。接下來本節(jié)將重點討論當網(wǎng)格單元大小變化時對數(shù)值結果的影響，以驗證SADUGKS格式的尺度自適應特性。

氣體參數(shù)、激波上下游初值等參數(shù)與4.1小節(jié)一致，CFL數(shù)取為0.9。速度空間采用32×32半平面Gauss-Hermite積分離散[38]。選取計算區(qū)域為-50Δx≤x≤50Δx，區(qū)域總長為L=100Δx，采用100個均勻網(wǎng)格將其離散，網(wǎng)格單元Δx選取了四種不同大小，分別是Δx=0.5λ1，λ1，10λ1和100λ1。圖2給出了無量綱密度、溫度以及權系數(shù)α在不同Δx時的空間分布曲線，其中橫坐標x*為采用激波上游的分子平均自由程λ1進行無量綱化后的空間位置，x*具體形式見式(54)，同時無黏Euler方程的精確解也畫在圖中作為對比。從圖2可以看出，當Δx較小，滿足Δx～λ1時，真實激波內部結構可以被分辨出來，此時SADUGKS表現(xiàn)為激波結構分辨格式(shock structure resolving scheme)；隨著Δx增加，如當Δx≥10λ1時，SADUGKS給出的結果將逼近無黏Euler方程的精確解，說明對于網(wǎng)格單元遠大于分子平均自由程時，激波被作為理想間斷面來處理，SADUGKS通過自身格式黏性最后得到的是數(shù)值激波，此時SADUGKS表現(xiàn)為激波捕捉格式(shock capturing scheme)。

這個例子說明SADUGKS可以根據(jù)比值Δx/λ1的變化而在激波結構分辨格式和激波捕捉格式之間做出自適應調整，這對于實際復雜工程問題獲得所用計算網(wǎng)格下的合理物理解是很有益的。因為一個問題在求解前進行網(wǎng)格劃分時并不知道流場物理量的局部細節(jié)特性，所以劃分的網(wǎng)格帶有較強的人為因素，但是SADUGKS格式自身可以通過一套統(tǒng)一的算法自動獲得流場在所用網(wǎng)格尺度下的合理物理解。

SADUGKS格式在激波結構分辨格式和激波捕捉格式之間的切換得益于權系數(shù)α的引入，從圖2權系數(shù)α的變化曲線可以清楚看出：當Δx?λ1時，α取值接近1，流動可認為是連續(xù)流動；當Δx～λ1時，α取值接近0.5，流動趨于稀薄流動，流動的非平衡效應增強。需要強調的是，通常所說的連續(xù)流動和稀薄流動都是一個相對的概念，和選用的網(wǎng)格單元尺度與分子平均自由程尺度相對大小有關。對同一個流動，網(wǎng)格單元尺度選用的不同，可獲得流動的信息量不同，描述流動的方法也是有區(qū)別的。網(wǎng)格單元尺度越小，觀察到的流動就越接近于一個個運動的微觀粒子，獲得的信息量就越豐富；而當網(wǎng)格單元尺度遠大于分子平均自由程尺度時，可獲得的流動信息是大量微觀粒子運動的統(tǒng)計平均，此時就是連續(xù)介質框架下的描述方法。通過上述討論說明SADUGKS是一種尺度自適應的氣體動理學格式，在該格式中，計算網(wǎng)格發(fā)揮著積極作用，最終可得到計算網(wǎng)格所對應尺度下的流動物理信息。

作為對比，對于傳統(tǒng)在連續(xù)介質框架下基于N-S方程建立的CFD方法來說，N-S方程建立后被作為一個數(shù)學偏微分方程來處理，網(wǎng)格進行加密的主要目的是減小N-S方程離散誤差，網(wǎng)格的作用未被充分發(fā)揮出來。另外，從物理上考慮，N-S方程是基于連續(xù)介質框架得到的，N-S方程蘊含的流動物理尺度是遠大于分子平均自由程的，因此即便將網(wǎng)格加密到分子平均自由程的量級，也只是獲得N-S方程作為一個數(shù)學偏微分方程更為精確的數(shù)值解，但所獲得的結果不可能超過N-S方程物理上成立的范圍，即無法給出在分子平均自由程尺度上的流動物理信息。

圖2 SADUGKS格式求解激波結構問題的密度、溫度與權系數(shù)α隨網(wǎng)格大小變化，CFL=0.9, Ma=1.2Fig.2 Density, temperature and α profiles of the shock structure with different cell sizes，CFL=0.9, Ma=1.2

4.3 從自由分子流領域到連續(xù)流領域

Sod問題為一維Riemann問題中的一個經(jīng)典問題，其無黏解包含激波、膨脹波及接觸間斷多種波系結構。該問題被廣泛用來檢驗數(shù)值格式捕捉間斷的能力，其初始條件如下：

(55)

本文這里考慮在上述初值下，引入氣體黏性后的流動。氣體考慮為空氣，分子內自由度設置為K=2，普朗特數(shù)Pr=0.72，比熱比γ=1.4。氣體黏性系數(shù)與溫度滿足μ=μ0(T/T0)w，其中w與分子模型相關，這里采用硬球模型w=0.5，μ0為參考黏性系數(shù)，T0為參考溫度。參考分子平均自由程λ0與參考黏性系數(shù)直接相關：

(56)

其中ρ0為參考密度。本文中，選用左側(x<0)初值ρ1、p1、T1分別作為參考密度、壓力和溫度。

μ0=10時的密度、溫度、速度以及此時流場中式(32)給出的權系數(shù)α隨空間分布如圖3(a)所示。在此黏性系數(shù)下，左側初值對應分子平均自由程λ0=12.77，比值λ0/Δx=1277，此時流動對應為自由分子流。從圖2中可以看出，SADUGKS結果與UGKS結果基本完全重合，且這種遠偏離平衡態(tài)的流動與無黏Euler解有著明顯的區(qū)別：一方面，由于流動更為稀薄，分子平均自由程遠大于網(wǎng)格單元，因此空間網(wǎng)格可以分辨包括激波在內的各種波系內部結構；另一方面，流動黏性很大，耗散作用增強，各種波系被顯著光滑抹平，物理量在整個空間上連續(xù)變化。圖中給出的權系數(shù)α在整個空間幾乎是一條平直線，值為0.5，這一點和前面第3節(jié)中提到的對于大Kn數(shù)的流動α→1/2結論一致。當黏性系數(shù)μ0減小到1時，左側初值對應分子平均自由程λ0=1.277，比值λ0/Δx=127.7，從圖3(b)可以看出，相比μ0=10時物理量分布變化不大，同樣此時權系數(shù)α在整個空間都接近0.5。

(a) μ0=10.0

(b) μ0=1.0

(c) μ0=0.1

(d) μ0=10-3

(e) μ0=10-5圖3 SADUGKS格式求解Sod問題的密度、溫度、速度和權系數(shù)α分布Fig.3 Distribution of density, temperature, velocity and α for Sod problem by SADUGKS

隨著黏性系數(shù)μ0減小到0.1時，左側初值對應分子平均自由程λ0=0.1277，比值λ0/Δx=12.77，相比前兩個較大黏性系數(shù)的工況來說，圖3(c)給出的物理量分布已出現(xiàn)一些變化，尤其速度峰值已明顯增加，并且權系數(shù)α在整個空間也出現(xiàn)起伏，不再是平直直線了。

當黏性系數(shù)μ0減小到10-3時，左側初值對應分子平均自由程λ0=1.277×10-3，比值λ0/Δx=0.1277，在此時的網(wǎng)格下已不能分辨激波內部結構了，從圖3(d)中可見，SADUGKS結果逐漸接近無黏Euler解。權系數(shù)α在整個空間已出現(xiàn)顯著變化，從最左側的0.5293逐漸減小到右側的0.5033。

當黏性系數(shù)μ0進一步減小到10-5時，左側初值對應分子平均自由程λ0=1.277×10-5，比值λ0/Δx=1.277×10-3，此時流動屬于連續(xù)流領域，并且由于黏性非常小，從圖3(e)可以看出，SADUGKS的結果幾乎與無黏Euler解完全重合。當前空間網(wǎng)格尺度遠大于分子平均自由程尺度，對應網(wǎng)格不再能分辨激波內部流動細節(jié)，激波表現(xiàn)為宏觀物理量的間斷。此時，SADUGKS表現(xiàn)為激波捕捉格式。從權系數(shù)α分布圖中可以看到α隨空間顯著變化，從左側的0.9716逐漸減小到右側的0.7660，這反映當?shù)嘏瓟?shù)和當?shù)靥卣鞒叨入S空間的變化。

本算例計算網(wǎng)格單元大小不變，改變μ0實現(xiàn)分子平均自由程λ0的變化，驗證了SADUGKS格式的尺度自適應特性。對于實際工程問題，使用一種能夠根據(jù)當?shù)亓鲃映叨冗M行自適應調節(jié)的統(tǒng)一算法顯得尤為重要，而SADUGKS格式正是這樣一種尺度自適應的算法，能夠反映出當?shù)亓鲃游锢硖匦浴Ｇ腋鶕?jù)當?shù)亓鲃映叨炔煌?，SADUGKS格式可在激波結構分辨格式和激波捕捉格式之間進行自適應切換。

4.4 高速可壓縮Couette流動問題

最后將SADUGKS格式應用于高速可壓縮Couette流動問題，流動參數(shù)與采用IP方法[39]進行計算的文獻一致。流動示意圖如圖4所示，兩無限大平行平板之間距離為1，兩板之間充滿氬氣，分子內自由度K=0，普朗特數(shù)Pr=2/3，比熱比γ=5/3，黏性系數(shù)與溫度相關μ=μ0(T/T0)w，采用變徑硬球模型w=0.81。兩板溫度固定為273 K，下板靜止，上板在自身平面內勻速運動，速度為U2=300 m/s，對應馬赫數(shù)Ma約為0.97。流體與壁面之間存在剪切作用，同時黏性耗散作用會使得流場內部溫度升高。通過改變初始氬氣密度以改變流動努森數(shù)Kn，本文共選取五組Kn=0.01，0.1，1.0，10和100，涵蓋了從連續(xù)流動到自由分子流的流動領域。不同Kn下流動將表現(xiàn)出不同的溫度及速度分布。

兩板間布置數(shù)目為31的均勻網(wǎng)格，速度空間采用28×28半平面Gauss-Hermite積分離散[38]。入口與出口采用周期性邊界條件，上下壁面采用完全漫反射邊界條件[19]。CFL數(shù)設置為0.5。

圖4 高速可壓縮Couette流動示意圖Fig.4 Schematic diagram for the compressible Couette flow

不同Kn下的速度剖面如圖5所示，從圖中可以看出，當Kn=0.01時，流動可認為仍滿足連續(xù)介質假設，但是上下壁面已出現(xiàn)輕微的速度滑移現(xiàn)象。隨著Kn不斷增加，非平衡效應變強，壁面處的滑移效應逐漸變得顯著，尤其當Kn=100時，板間氣體速度分布趨于均勻，約為150 m/s，和上、下板速度遠不相同。

圖5 SADUGKS計算所得不同Kn下的速度剖面圖Fig.5 Velocity profiles for the compressible Couette flow at different Knudsen numbers

圖6所示為不同努森數(shù)下的溫度剖面。當Kn=0.01時，流動在壁面已出現(xiàn)輕微溫度跳躍現(xiàn)象。由于上板克服壁面黏性剪切做功最終轉化為熱能，并逐漸向流道中間集中，造成溫度峰值出現(xiàn)在中心對稱線上。隨著Kn不斷增加，非平衡效應變強，壁面處的溫度跳躍現(xiàn)象更為明顯，板間氣體的溫度梯度逐漸減小。尤其當Kn=100時，板間流動趨于自由分子流，板間氣體溫度幾乎相同，約為308 K，和上、下板溫度遠不相同。

為了反映氣體Pr數(shù)對結果的影響，這里考慮了Pr=1時Kn=0.01和Kn=0.1兩種工況，并和前面Pr=2/3的對應結果進行了對比，見圖7給出的溫度分布曲線，可以清楚看出Pr=1時的結果與Pr=2/3及文獻IP方法[39]的結果出現(xiàn)了較大的偏差。從式(10～14)可以看出，Boltzmann-Shakhov模型中取Pr=1時退化為Boltzmann-BGK模型，這說明為了正確反映氣體Pr數(shù)的影響從而獲得準確溫度場，基于Boltzmann-Shakhov模型構造計算格式是適當?shù)摹?/p>

無論對于速度還是溫度，在所有Kn情況下，SADUGKS格式所得結果均能與文獻IP方法[39]的結果很好地吻合，表明了SADUGKS格式對于寬Kn流動計算的有效性和尺度自適應特性。

圖6 SADUGKS計算所得不同Kn下的溫度剖面圖Fig.6 Temperature profiles for the compressible Couette flow at different Knudsen numbers

圖7 Pr=1和Pr=2/3的結果對比(Kn=0.01和Kn=0.1)Fig.7 Temperature profiles for the compressible Couette flow at different Knudsen numbers with Pr=1 and Pr=2/3

5 結 論

文中提出了Boltzmann方程沿特征線離散的一般形式，式中碰撞項的離散采用顯式和隱式加權平均的方法，權系數(shù)α依賴于當?shù)嘏鲎菜沙跁r間和計算時間步長的比值。通過權系數(shù)α的引入，對文獻現(xiàn)有離散統(tǒng)一氣體動理學格式(DUGKS)進行了改進，提出了具有尺度自適應特性的離散統(tǒng)一氣體動理學格式(SADUGKS)。對SADUGKS進行分析，明確權系數(shù)α與反映當?shù)亓鲃映叨鹊漠數(shù)嘏瓟?shù)直接相關，表明SADUGKS格式能夠根據(jù)當?shù)亓鲃映叨冗M行自適應調節(jié)。SADUGKS格式的自適應特性和有效性通過典型可壓縮流動進行了檢驗，取得了較好的計算結果，表明SADUGKS格式適用于寬范圍Kn變化、跨流域的多尺度流動問題的數(shù)值模擬。后續(xù)將選用更為復雜的流動問題對SADUGKS做進一步的檢驗和驗證。