考慮分子轉(zhuǎn)動自由度的離散統(tǒng)一氣體動理學(xué)格式

姚博，張創(chuàng)，郭照立

華中科技大學(xué) 煤燃燒國家重點(diǎn)實(shí)驗(yàn)室，武漢 430074

多尺度流動在工程應(yīng)用中廣泛存在，例如臨近空間飛行器、微機(jī)電系統(tǒng)(Micro-Electro-Mechanical System,MEMS)等問題中，同一個計(jì)算域內(nèi)會同時存在連續(xù)流、滑移流甚至是自由分子流等不同流域，局部努森數(shù)可能相差數(shù)個量級[1]。在稀薄流中，基于連續(xù)介質(zhì)假設(shè)的Navier-Stokes方程不再成立。常規(guī)求解Boltzmann方程的計(jì)算格式往往采用算子分裂方法[2-4]，存在時間步長小于平均碰撞時間、網(wǎng)格尺寸小于分子平均自由程的限制。而多尺度問題中連續(xù)流的存在，使得直接使用這類方法計(jì)算整個流場消耗相當(dāng)巨大。

多尺度流動問題的直觀解決方法是將不同求解器耦合起來。在局部努森數(shù)較大的計(jì)算區(qū)間中求解Boltzmann方程，在局部努森數(shù)較小的計(jì)算區(qū)間中求解Navier-Stokes方程[5-6]。但準(zhǔn)確劃分不同流域、合理地將不同流域耦合起來都很困難。近年來，基于動理學(xué)理論發(fā)展出的一類漸進(jìn)收斂(Asymptotic Preserving, AP)格式[7]，使得同一個算法處理整個多尺度問題成為可能。其中，Xu等提出的統(tǒng)一氣體動理學(xué)格式(Unified Gas Kinetic Scheme,UGKS)[8-10]采用有限體積形式，在計(jì)算網(wǎng)格界面通量時采用控制方程在時間上的局部積分解，耦合了粒子間的碰撞與遷移過程，能夠準(zhǔn)確描述連續(xù)流到稀薄流流域的流動問題。Guo等提出的離散統(tǒng)一氣體動理學(xué)格式(Discrete Unified Gas Kinetic Scheme,DUGKS)[11-13]也采用有限體積形式，繼承了UGKS的主要優(yōu)點(diǎn)，但與UGKS也存在明顯的不同。在構(gòu)造界面通量時，DUGKS采用的是沿特征線積分得到的數(shù)值解，計(jì)算效率更高[14]。

在實(shí)際的多尺度流動問題中，氣體往往是多原子氣體，如空氣的最主要成分是氧氣和氮?dú)?種雙原子氣體。氧氣和氮?dú)獾霓D(zhuǎn)動特征溫度分別為2.07 K和2.88 K，因此即使在室溫下其轉(zhuǎn)動能也能被充分激發(fā)出來。而由于振動特征溫度具有較高的值(氧氣與氮?dú)夥謩e為2 256 K和2 719 K)，只有達(dá)到較高的溫度振動能才會被激發(fā)[15-17]。Rykov模型[18]就是考慮了分子平動和轉(zhuǎn)動自由度的動理學(xué)方程，能適用于較大范圍的雙原子氣體流動問題研究。近年來，不少工作基于Rykov模型方程研究了雙原子氣體流動問題[19]。例如Rykov等[20-22]先后計(jì)算了雙原子氣體的正激波結(jié)構(gòu)和微管道內(nèi)的流動；Liu等[23]構(gòu)造出了基于Rykov模型方程的統(tǒng)一UGKS，計(jì)算的正激波結(jié)構(gòu)、外掠平板等多個算例與直接蒙特卡羅(Direct Simulation Monte Carlo,DSMC)和實(shí)驗(yàn)結(jié)果吻合良好；Wu等[24]將Rykov模型方程中的彈性碰撞項(xiàng)用Boltzmann碰撞項(xiàng)替換，并采用快速譜方法進(jìn)行求解。

本文構(gòu)造出了基于Rykov模型方程的離散統(tǒng)一氣體動理學(xué)格式，從單原子氣體拓展為計(jì)算考慮了分子轉(zhuǎn)動自由度的雙原子氣體，能夠有效模擬雙原子氣體從連續(xù)到稀薄的多尺度流動問題。

1 Rykov模型方程

Rykov模型方程考慮了氣體分子內(nèi)部的轉(zhuǎn)動自由度，分別處理分子的彈性碰撞和非彈性碰撞[18]。氣體分子的速度分布函數(shù)記為f=f(t,x,ξ,η,e)，其中t為時間,x為空間坐標(biāo),ξ為分子在D維空間中的平動速度,η為分子速度在三維空間中的其他分量,e=M2/(2J)為分子轉(zhuǎn)動能，M為角動量，J為轉(zhuǎn)動慣量。宏觀量可以通過分布函數(shù)定義為

(1a)

(1b)

(1c)

(1d)

(1e)

(1f)

q=qr+qt

(1g)

式中:c為分子熱速度，c=ξ-U;k為玻爾茲曼常數(shù);n為分子數(shù)密度;m為分子質(zhì)量;U為流體速度;Tt和Tr分別為平動和轉(zhuǎn)動溫度;ε為分子轉(zhuǎn)動自由度對應(yīng)的平均能量;qt為與平動自由度對應(yīng)的熱流;qr為轉(zhuǎn)動能量傳輸產(chǎn)生的熱流;ρ為密度。

在不考慮外力項(xiàng)的情形下，Rykov模型方程中氣體分子速度分布函數(shù)f=f(t,x,ξ,η,e)的演化方程有如下形式：

(2)

式中:右側(cè)兩項(xiàng)分別表示非彈性碰撞項(xiàng)和彈性碰撞項(xiàng)。彈性碰撞項(xiàng)保持平動動能守恒，非彈性碰撞項(xiàng)則實(shí)現(xiàn)平動和轉(zhuǎn)動二者間能量的轉(zhuǎn)換。νt和νr分別為彈性和非彈性碰撞頻率；ft和fr分別為彈性碰撞項(xiàng)和非彈性碰撞項(xiàng)的平衡態(tài)：

(3a)

(3b)

(3c)

(3d)

式中:松弛時間τ=μt/pt,μt和pt分別為溫度Tt對應(yīng)的動力黏度和壓力;T為平衡態(tài)溫度;參數(shù)Z反映單一轉(zhuǎn)動碰撞中平動碰撞發(fā)生的次數(shù);D0為雙原子氣體的自擴(kuò)散系數(shù);ω0和ω1用于調(diào)整熱流qt和qr至合適的松弛值：

(4a)

(4b)

Rykov在文獻(xiàn)[18]中將分布函數(shù)對轉(zhuǎn)動能e求積，引入2個約化的分布函數(shù)，式(2)則轉(zhuǎn)化為2個關(guān)于約化后分布函數(shù)的方程。在D維問題中，分布函數(shù)的演化僅與D維的離散速度ξ有關(guān)而與η無關(guān)。為進(jìn)一步消除分布函數(shù)對冗余變量的依賴，降低計(jì)算消耗，可采用類似的處理方法引入3個約化的分布函數(shù)

(5a)

(5b)

(5c)

則宏觀量可表示為

(6a)

(6b)

(6c)

(6d)

(6e)

(6f)

q=qr+qt

(6g)

式中:R為氣體常數(shù)。

分布函數(shù)f0、f1和f2的控制方程可對式(2)積分得到

(7a)

(7b)

(7c)

(8a)

(8b)

(8c)

(8d)

(8e)

(8f)

(8g)

(9a)

(9b)

(9c)

控制方程式(7)可以改寫為

(10a)

(10b)

(10c)

2 數(shù)值格式

本節(jié)基于式(10)構(gòu)建考慮分子轉(zhuǎn)動自由度的離散統(tǒng)一氣體動理學(xué)格式。首先將其統(tǒng)一為

(11)

(12)

(13)

式中：Fn+1/2(ξ)為穿過網(wǎng)格界面的流量；Δt=tn+1-tn為時間步長；|Vj|和?Vj分別為控制體Vj的體積和表面積；dS為表面積；r為指向控制體表面外的單位法向量;φj和Ωj為單元上分布函數(shù)和碰撞項(xiàng)的平均值：

(14a)

(14b)

(15a)

(15b)

于是式(12)轉(zhuǎn)變?yōu)?/p>

(16)

由式(11)中碰撞項(xiàng)的守恒性可以得到

(17a)

(17b)

(17c)

(17d)

(17e)

(17f)

q=qr+qt

(17g)

計(jì)算界面流量Fn+1/2的關(guān)鍵在于獲得界面處原始分布函數(shù)fn+1/2。DUGKS采用在半個時間步長s=Δt/2內(nèi)對式(16)沿特征線積分：

φ(tn+s,xb,ξ)-φ(tn,xb-ξs,ξ)=

(18)

(19a)

(19b)

則式(18)可以寫為

(20)

其中：

(xb-xj-ξs)·σj，(xb-ξs)∈Vj

(21)

(22a)

(22b)

(22c)

(22d)

(22e)

(22f)

q=qr+qt

(22g)

通過tn+1/2時刻界面處宏觀量可以得到的相應(yīng)的Rykov分布函數(shù)φR，依據(jù)式(19)得到原始的分布函數(shù)：

(23)

在上述推導(dǎo)中，分子運(yùn)動速度ξ被認(rèn)為是連續(xù)的。在實(shí)際計(jì)算中，需要將速度空間離散為速度集ξi(i=1,2,…,N)。離散速度集通常選為Gaussian-Hermite或Newton-Cotes等特定數(shù)值積分公式的節(jié)點(diǎn)坐標(biāo)，上述積分需要被替換為數(shù)值積分。例如，對于守恒變量有

(24a)

(24b)

(24c)

(24d)

由于在重構(gòu)界面分布函數(shù)時耦合了分子的碰撞和遷移過程，DUGKS具備漸進(jìn)收斂特性。DUGKS中的時間步長由Courant-Friedrichs-Lewy (CFL)條件確定[11]：

(25)

式中:ζ為CFL數(shù);Δx為最小網(wǎng)格尺寸;ξm為最大離散速度值;Um為最大流動速度。

對于定常流動問題，選取2個相鄰時間步中宏觀量場的平均變化比α作為計(jì)算收斂的判據(jù):

(26)

式中:W∈{ρ,U,T}，當(dāng)α小于給定的參考值時認(rèn)為計(jì)算達(dá)到收斂。

Rykov模型方程中參數(shù)Z與旋轉(zhuǎn)碰撞數(shù)Zrot之間存在如下關(guān)系[24]:

(27)

其中，當(dāng)計(jì)算雙原子時固定參數(shù)d=2。Zrot可采用Landau-Teller-Jeans轉(zhuǎn)動能量松弛模型[25]計(jì)算得到，即

(28)

3 數(shù)值驗(yàn)證

本節(jié)將采用數(shù)個雙原子氣體流動的問題驗(yàn)證基于Rykov模型方程的DUGKS方法的正確性和有效性，包括一維氮激波結(jié)構(gòu)、二維超聲速平板繞流、二維超聲速圓柱繞流和2個連通方腔中的跨流域流動等問題。

3.1 激波結(jié)構(gòu)

(29a)

(29b)

(29c)

(29d)

式中：Ma為上游馬赫數(shù)；Ma′為下游馬赫數(shù)。

計(jì)算域選為-25λ1≤x≤25λ1，并采用100個均勻分布的物理網(wǎng)格。速度空間取為-15≤ξ≤15，采用201個均勻分布的Newton-Cotes積分點(diǎn)。x≤0范圍內(nèi)的分布函數(shù)初始值為上游參數(shù)確定的平衡態(tài)分布，x>0范圍內(nèi)的分布函數(shù)初始值為下游參數(shù)確定的平衡態(tài)分布。平動溫度Tt和轉(zhuǎn)動溫度Tr的初始值與平衡態(tài)溫度T相等，即Tt1=Tr1=T1和Tt2=Tr2=T2。激波上下游邊界處分布函數(shù)分別固定為上下游宏觀量對應(yīng)的平衡態(tài)分布。計(jì)算中CFL數(shù)均取為0.95，收斂判據(jù)為熱流的平均變化比小于10-8，即αq<10-8。

DSMC方法中取ZDSMC=3.5，DUGKS和UGKS當(dāng)中Z=2.4，根據(jù)文獻(xiàn)[24]取參數(shù)δ=1/1.33,ω0=0.477以及ω1=1.919。分別比較了Ma=1.53、4.0、5.0和7.0 這4種馬赫數(shù)下的氮激波結(jié)構(gòu)。圖1給出不同馬赫數(shù)下，DUGKS與文獻(xiàn)[23]中UGKS和DSMC的氮激波結(jié)構(gòu)計(jì)算結(jié)果對比。圖1中均為歸一化后的宏觀量，以Q表示宏觀量，下標(biāo)u和d分別表示上游和下游值，歸一化后為Q′=(Q-Qu)/(Qd-Qu)?？梢钥吹剑珼UGKS與UGKS的計(jì)算結(jié)果在不同馬赫數(shù)下都吻合良好。二者計(jì)算的密度曲線與DSMC的結(jié)果吻合良好，溫度曲線則在激波上游位置上升得更早一些。產(chǎn)生這一差異的主要原因是DUGKS和UGKS采用的松弛模型只采用了單一的松弛時間，在靠近激波上游的位置將高速分子的能量平均分配給全部分子[12]。

圖1 不同馬赫數(shù)下DUGKS、UGKS和DSMC的氮激波結(jié)構(gòu)計(jì)算結(jié)果對比Fig.1 Comparison of nitrogen shock structure results of DUGKS, UGKS, and DSMC under different Mach numbers

3.2 超聲速平板繞流

飛行器在臨近空間飛行中以高超聲速飛行時，流場中存在激波-激波、激波-壁面的相互作用，飛行器表面附近會出現(xiàn)高熱流和高壓強(qiáng)的現(xiàn)象。

本算例取自Tsuboi和Matsumoto[26]的風(fēng)洞試驗(yàn)run34，Xu等[27]基于多溫度模型方程的氣體動理學(xué)格式，Liu等[23]基于Rykov模型方程的統(tǒng)一氣體動理學(xué)格式，先后對該問題進(jìn)行了數(shù)值模擬。在run34試驗(yàn)中，噴嘴出口馬赫數(shù)Ma=4.89，總溫T0=670 K，總壓p0=983 Pa，噴嘴出口溫度Te=116 K，密度ρe=6.15×10-5kg/m3。平板由金屬銅制造而成，試驗(yàn)中通過水冷實(shí)現(xiàn)恒溫290 K。試驗(yàn)氣體為氮?dú)?，氣體常數(shù)R=297 J·kg-1·K-1，熱度系數(shù)ω=0.75。黏性與溫度之間關(guān)系為

(30)

式中：μ0=1.656×10-5Pa·s,T0=273.5 K。計(jì)算中采用硬球分子模型，平均自由程λ與黏性的關(guān)系為

(31)

圖2 超聲速平板繞流計(jì)算網(wǎng)格Fig.2 Computational mesh for supersonic flow around flat plate

圖3(a)～圖3(d)分別為計(jì)算得到的密度、平均溫度、平動溫度以及轉(zhuǎn)動溫度云圖?？梢钥吹剑D(zhuǎn)動溫度Tr升高滯后于平動溫度Tt。這一現(xiàn)象可以解釋為：高速的雙原子氣體與壁面碰撞后，宏觀動能轉(zhuǎn)換為分子熱運(yùn)動動能，使平動溫度升高。經(jīng)過分子間的相互作用，再傳遞給轉(zhuǎn)動溫度對應(yīng)的內(nèi)能。在平板上端x=5 mm和x=20 mm處沿坐標(biāo)y方向上的溫度曲線在圖4中給出，計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)、UGKS解總體吻合良好。

圖3 DUGKS計(jì)算的平板繞流密度、平均溫度、平動溫度和轉(zhuǎn)動溫度云圖Fig.3 Contours of DUGKS results around flat plate for density, equilibrium temperature, translational temperature, and rotational temperature

圖4 x=5 mm和x=20 mm處沿y方向溫度分布曲線對比Fig.4 Temperature plots along y direction at x=5 mm and x=20 mm

3.3 超聲速圓柱繞流

本算例進(jìn)一步測試了氮?dú)獾某曀兮g體繞流問題，算例來源于文獻(xiàn)[23]。來流氮?dú)獾乃俣萓∞=1 684.5 m/s，溫度T∞=273 K，分子數(shù)密度n∞=1.299 44×1021/m3，動力黏度μ∞=1.657 88×10-5Pa·s，取熱度系數(shù)ω=0.75。圓柱表面為恒定溫度Tw=273 K?？捎?jì)算出，來流氮?dú)怦R赫數(shù)Ma=5.0，努森數(shù)Kn=0.1。

圖5 超聲速圓柱繞流計(jì)算網(wǎng)格Fig.5 Computational mesh of supersonic flow around a cylinder

圖6(a)～圖6(f)分別為DUGKS計(jì)算圓柱繞流的密度、x方向速度、y方向速度、平均溫度、平動溫度以及轉(zhuǎn)動溫度云圖。圖7顯示的是在y=0 m處沿x方向的參數(shù)曲線?？梢园l(fā)現(xiàn)，DUGKS計(jì)算的密度、速度曲線都與UGKS和DSMC結(jié)果吻合良好，而DUGKS和UGKS二者的溫度曲線比DSMC的計(jì)算結(jié)果上升地更早一些，這與前面的一維激波結(jié)構(gòu)的研究結(jié)論是一致的。

圖6 DUGKS計(jì)算的超聲速圓柱繞流流場云圖Fig.6 Contours of supersonic flow around a cylinder of DUGKS results

圖7 超聲速圓柱繞流在y=0 m位置的流場分布Fig.7 Flow distributions along line y=0 m for supersonic flow around a cylinder

3.4 兩個連通方腔中的跨流域流動

在太空探索中，火箭燃料罐爆裂或者宇航倉泄漏時，泄漏到太空中的氣體會經(jīng)歷從連續(xù)流到自由分子流等不同流域的變化。高效地模擬跨流域流動對預(yù)測航天器泄漏等實(shí)際問題具有指導(dǎo)意義。本算例參照文獻(xiàn)[28]構(gòu)造出了2個方腔中的跨流域流動，如圖8所示。方腔及管道中氣體為氮?dú)?。方腔邊長為L=1 m，管道長為L，寬為H=L/8。初始狀態(tài)下，左、右兩個方腔分別維持在不同壓力，管道中間處的隔板將兩側(cè)氣體分開。兩側(cè)氣體溫度均為273 K，左側(cè)壓力pL=48.16 Pa，右側(cè)壓力pR=0.004 816 Pa。方腔及管道壁面均為恒定溫度Tw=273 K。可計(jì)算得到，左、右兩側(cè)的努森數(shù)分別為KnL=0.001和KnR=10，左側(cè)處于連續(xù)流域，右側(cè)接近自由分子流域。在t=0 s時刻，移去管道中的隔板，左側(cè)的氣體在壓力的作用下迅速噴入右側(cè)。

圖8 兩個連通方腔跨流域流動示意圖Fig.8 Schematic diagram of multiscale flow of two connected cavities

圖9 兩個連通方腔跨流域流動的計(jì)算網(wǎng)格Fig.9 Computational mesh of multiscale flow of two connected cavities

圖10 t/tc=1時馬赫數(shù)及壓力分布云圖Fig.10 Contours of Mach number and pressure at t/tc=1

圖11 t/tc=0.09，0.5，1.0和2.0時刻沿中軸線(y=0 m)溫度分布曲線Fig.11 Temperature distributions along line y=0 m at time t/tc=0.09, 0.5, 1.0, and 2.0

4 結(jié) 論

基于考慮分子轉(zhuǎn)動自由度的離散統(tǒng)一氣體動理學(xué)格式，對激波結(jié)構(gòu)、超聲速平板繞流以及超聲速圓柱繞流等算例進(jìn)行計(jì)算，得到以下結(jié)論：

1) 基于Rykov動理學(xué)模型方程并采用Landau-Teller-Jeans轉(zhuǎn)動能量松弛模型的離散統(tǒng)一氣體動理學(xué)格式，可以反映出雙原子氣體非平衡流動時轉(zhuǎn)動和平動自由度對應(yīng)的能量轉(zhuǎn)換過程，可用于計(jì)算雙原子氣體的多尺度流動問題。

2) DUGKS計(jì)算的結(jié)果與UGKS、DSMC的解以及實(shí)驗(yàn)值吻合良好。其中，激波結(jié)構(gòu)算例中馬赫數(shù)范圍從1.5變化到7.0，超聲速外掠平板以及超聲速圓柱繞流分別測試了DUGKS在帶有尖端的平板以及鈍體中的計(jì)算性能，驗(yàn)證了考慮分子轉(zhuǎn)動自由度的離散統(tǒng)一氣體動理學(xué)格式的準(zhǔn)確性及穩(wěn)定性。

本文僅構(gòu)造了考慮分子轉(zhuǎn)動自由度的離散統(tǒng)一氣體動理學(xué)格式并驗(yàn)證了一維及二維算例，同時考慮分子振動自由度并將其應(yīng)用于實(shí)際的復(fù)雜三維計(jì)算有待進(jìn)一步研究。