基于修正N-S方程本構(gòu)關(guān)系的氣體動(dòng)理論耦合方法

皮興才，李志輝，2，*，彭傲平，吳俊林，蔣新宇

(1. 中國空氣動(dòng)力研究與發(fā)展中心 超高速空氣動(dòng)力研究所，綿陽 621000；2. 國家計(jì)算流體力學(xué)實(shí)驗(yàn)室，北京 100191)

0 引 言

近年來，飛行高度在40～70 km、馬赫數(shù)在5～20的近空間飛行器的研發(fā)受到重視[1]。在此環(huán)境下飛行普遍面臨著局部稀薄效應(yīng)導(dǎo)致的連續(xù)介質(zhì)假設(shè)失效的問題。因此，考慮氣體分子輸運(yùn)統(tǒng)計(jì)分布函數(shù)演化特性的介觀方法獲得了廣泛的關(guān)注。Boltzmann方程是介觀方法的理論基礎(chǔ)[2]。在研究應(yīng)用中，為了定量化表征處理Boltzmann方程復(fù)雜的碰撞項(xiàng)，很多簡化的模型方程被提出，包括BGK模型方程[3]、Shakhov模型方程[4]、ES-BGK模型方程[5]、Rykov模型方程[6]等，所能模擬的氣體流動(dòng)介質(zhì)也從最初的單原子理想氣體發(fā)展到多原子、多組分，考慮轉(zhuǎn)動(dòng)、振動(dòng)激發(fā)及熱化學(xué)非平衡效應(yīng)的跨流域空氣動(dòng)力學(xué)模擬體系日趨完善[7-17]。

模型方程是描述氣體分子速度分布函數(shù)關(guān)于時(shí)間、空間及速度相空間七個(gè)維度的偏微分方程。由于方程的數(shù)值求解涉及到對七個(gè)維度的離散，提高計(jì)算效率一直是亟待解決的難題。針對BGK模型、ESBGK模型方程，Mieusses[18]通過構(gòu)造具有保守恒特性的離散平衡分布函數(shù)解析表達(dá)形式，建立了嚴(yán)格守恒的數(shù)值方法，使得該方法能以較少的離散速度點(diǎn)實(shí)現(xiàn)方程求解。Chen、Huang等[19-20]同樣立足于減少離散速度空間的大小，采用叉樹數(shù)據(jù)結(jié)構(gòu)對UGKS的離散速度空間進(jìn)行優(yōu)化以提高計(jì)算效率。進(jìn)一步地，2016年，江定武[21]針對顯式UGKS計(jì)算效率低下的問題，利用CFD領(lǐng)域的經(jīng)典數(shù)值方法構(gòu)造了隱式格式，并采用MPI+OpenMP并行計(jì)算策略以改善UGKS的計(jì)算效率及收斂速度。在耦合方法方面，2012年Alessandro等[22]采用DSMC方法對單原子ES-BGK方程進(jìn)行數(shù)值求解，并通過區(qū)域剖分實(shí)現(xiàn)了與Euler方程的耦合求解，并將其應(yīng)用到一維間斷問題分析模擬。2014年Rovenskaya等[23]采用有限差分法與離散速度坐標(biāo)法結(jié)合，開展了Shakhov模型方程與NS模型方程的耦合求解，并將其應(yīng)用于低速稀薄氣體流動(dòng)分析，實(shí)現(xiàn)了任意壓比下的薄板稀薄繞流模擬。最近，基于多重網(wǎng)格技術(shù)及采用宏觀方程預(yù)測平衡分布函數(shù)的全隱UGKS成功應(yīng)用于跨流域稀薄流模擬，大幅提升了傳統(tǒng)UGKS的計(jì)算效率[24]。

在氣體動(dòng)理學(xué)理論框架下，N-S方程可以通過Boltzmann方程或模型方程的Chapman-Enskog展開一階近似。這體現(xiàn)了N-S方程與模型方程的理論可銜接性，也體現(xiàn)出N-S方程的質(zhì)量、動(dòng)量及能量守恒的基本特性即使在稀薄非平衡效應(yīng)加重的近連續(xù)流區(qū)依舊成立，而對于稀薄效應(yīng)更為嚴(yán)重、氣體分子速度分布函數(shù)呈現(xiàn)顯著非平衡特征的流動(dòng)，N-S方程所遵從的本構(gòu)關(guān)系難以正確表征模擬時(shí)，通常采用滑移邊界可解決邊界條件失效問題，但針對本構(gòu)關(guān)系失效卻是難以修正。這也是二階滑移邊界條件，如Cercignan模型[25-26]、Deissler模型[27]等的模擬精度相較于一階滑移邊界條件的提升有限的根本原因[28-29]。為了描述更復(fù)雜的本構(gòu)關(guān)系，肖洪[30]在廣義流體力學(xué)方程（GHE）的基礎(chǔ)上提出了非線性耦合本構(gòu)關(guān)系（NCCR），并在激波結(jié)構(gòu)、二維高超聲速稀薄繞流等問題模擬中獲得成功。Burnett方程、Grad矩方程同樣立足于通過理論建模的方式獲得非守恒量—應(yīng)力張量及熱流的輸運(yùn)關(guān)系，也取得了成功的應(yīng)用[31]。

為了從數(shù)值層面構(gòu)建本構(gòu)關(guān)系，拓寬N-S方程在稀薄氣體動(dòng)力學(xué)方面的模擬能力，規(guī)避復(fù)雜的宏觀非守恒量輸運(yùn)方程帶來的求解難題，本研究以多原子氣體ES-BGK模型方程為基礎(chǔ)，采用氣體動(dòng)理論統(tǒng)一算法（GKUA）[7-17]描述非線性本構(gòu)關(guān)系及壁面速度滑移、溫度跳躍，修正N-S方程的線性本構(gòu)關(guān)系及滑移邊界條件，拓展其在過渡流、滑移流域的應(yīng)用能力。在后續(xù)內(nèi)容中，本文首先通過Chapman-Enskog基于Maxwell平衡態(tài)分布漸近展開給出ES-BGK方程當(dāng)?shù)仄胶鈶B(tài)分布函數(shù)、本構(gòu)關(guān)系的一階近似，其次，構(gòu)造基于N-S方程本構(gòu)關(guān)系修正及滑移邊界耦合的數(shù)值求解方法，并通過可壓縮平板邊界層流動(dòng)、圓柱繞流經(jīng)典案例來進(jìn)行驗(yàn)證。

1 理論模型

經(jīng)典熱力學(xué)理論認(rèn)為分子內(nèi)能可以通過分子的自由度及溫度來表征?；谠摾碚?，Brull等[5,16]在單原子ES-BGK模型的基礎(chǔ)上通過引入內(nèi)能自由度eint作為分布函數(shù)的自變量，即分布函數(shù)f=f(t,x,ξ,eint)構(gòu)造了一種適用于多原子分子氣體的ES-BGK類模型方程：

式中，

其中，δ為轉(zhuǎn)動(dòng)自由度；Z為轉(zhuǎn)動(dòng)松弛參數(shù)；I為單位矩陣；b∈[-0.5,1) ；氣體分子熱運(yùn)動(dòng)速度C是氣體分子隨機(jī)速度ξ與宏觀速度u之差，其平衡態(tài)分布函數(shù)為：

式中，

對ES-BGK模型方程進(jìn)行基于Maxwell平衡態(tài)分布的Chapman-Enskog漸近分析?；诹烤V原理將式(1)改寫為如下形式[26,32]：

式中，無量綱參數(shù) ε是與Knudsen數(shù)同階的量，用于衡量松弛碰撞、對流及擴(kuò)散等物理過程的時(shí)間尺度。

對分布函數(shù)及當(dāng)?shù)仄胶鈶B(tài)分布函數(shù)進(jìn)行漸近展開，具有如下形式：

由Chapman-Enskog漸近分析需要遵循的約束條件：

可知f(0)=fM。將Ttr及Tint同樣進(jìn)行多尺度展開：

式中，ΔTtr及ΔTint表示非平衡效應(yīng)引起的平動(dòng)及轉(zhuǎn)動(dòng)溫度增量，它們滿足如下關(guān)系：

對于非守恒量—應(yīng)力張量及熱流進(jìn)行多尺度展開：

式中，

為了獲得當(dāng)?shù)仄胶夥植己瘮?shù)G的一階近似表達(dá)式，對式(2)中 det(Γ) 及 Γ-1進(jìn)行多尺度展開分析。根據(jù)前述關(guān)系，整理Γ表達(dá)式可得：

式中：

將式(24)帶入式(29)，并保留d et(Γ) 及Γ-1的Taylor級數(shù)到O(ε)可得：

針對量熱完全氣體不考慮ΔTtr及 ΔTint對結(jié)果的影響，對當(dāng)?shù)仄胶鈶B(tài)G進(jìn)行Taylor開展，其O(1)及O(ε)階項(xiàng)如下所示：

進(jìn)一步地，將f及G的展開式(18)及式(19)代入式(17),可知：

等式左端等于[5]：

式中，

可得：

因此，利用式(9)可知分布函數(shù)的Chapman-Enskog漸近展開一階近似表達(dá)式為：

式中，

可見，ES-BGK模型的應(yīng)力張量σ及 熱流q的一階漸近展開結(jié)果與N-S方程的本構(gòu)關(guān)系相同。由此可知，ES-BGK方程與N-S方程在連續(xù)流域的物理描述上具有理論一致性。進(jìn)一步說明，N-S方程在描述稀薄氣體流動(dòng)問題時(shí)，其質(zhì)量、動(dòng)量及能量守恒的基本特性依舊有效，其線性本構(gòu)關(guān)系是方程失效的核心。這也為研究者拓展N-S方程的應(yīng)用范圍提供了理論支持?；谏鲜鲫P(guān)系，本研究提出一種基于修正NS方程本構(gòu)關(guān)系的耦合氣體動(dòng)理論的方法，以保證耦合方法在質(zhì)量、動(dòng)量及能量守恒性，拓展N-S方程在局部稀薄效應(yīng)氣體流動(dòng)問題模擬的能力。

2 數(shù)值方法

本耦合方法基于課題組氣體動(dòng)理論統(tǒng)一算法(GKUA)及開源代碼OpenCFD-EC2D的NS有限體積求解器構(gòu)建。作為原理驗(yàn)證研究，本文中耦合區(qū)域通過計(jì)算調(diào)整給定。其中，GKUA求解需要進(jìn)行N-S方程本構(gòu)關(guān)系修正的區(qū)域—壁面附近與平均分子自由程相當(dāng)?shù)膮^(qū)域，N-S求解器對全流區(qū)進(jìn)行守恒方程求解，如圖1所示。

圖1 求解域分區(qū)示意圖Fig. 1 Schematic view of the computational domain decomposition

氣體動(dòng)理論耦合策略大致可分如下幾步：

1）通過有限體積法全流域求解N-S方程；

2）通過GKUA方法求解壁面區(qū)域，耦合邊界處的分布函數(shù)根據(jù)特征線n指向給定分布函數(shù)的邊界值。其中， (n·ξ＞0)的離散分布函數(shù)采用外推方式給定；(n·ξ≤0)的離散分布函數(shù)由N-S方程計(jì)算給出的宏觀物理量，通過式（38）計(jì)算離散分布函數(shù)給定。壁面分布函數(shù)按完全漫反射條件給出；

3）求解GKUA計(jì)算域的應(yīng)力張量 σGKUA、熱流矢量qGKUA及 壁面速度滑移uslip及溫度跳躍Tjump條件；

4）利用獲得的 σGKUA、qGKUA替換GKUA域的NS計(jì) 算 結(jié)果 σNS及qNS， 利 用uslip及Tjump給 定N-S求解域邊界條件；

5）以修正后的本構(gòu)關(guān)系及壁面邊界進(jìn)行下一時(shí)間步N-S方程及GKUA求解。

為了降低對計(jì)算資源的需求，采用離散速度坐標(biāo)法直接數(shù)值求解Boltzmann模型方程通常對分布函數(shù)的自變量進(jìn)行約化[7-8]，引入變量：

對于二維流動(dòng)問題，引入下式對ES-BGK模型方程進(jìn)行進(jìn)一步的約化：

在空間格式構(gòu)造方面，本文采用二階NND格式對N-S方程及ES-BGK模型方程進(jìn)行界面重構(gòu)[7-9]，采用Steger-Warming格式進(jìn)行界面通量計(jì)算，并且采用LU-SGS隱式格式對時(shí)間項(xiàng)進(jìn)行離散以提高效率。

3 算例驗(yàn)證

為了驗(yàn)證構(gòu)建的基于修正本構(gòu)關(guān)系的氣體動(dòng)理論耦合方法，本文分別進(jìn)行了可壓縮平板邊界層流動(dòng)模擬及稀薄圓柱繞流模擬。其中，DSMC計(jì)算結(jié)果由DS2V軟件[33]計(jì)算得到。

3.1 可壓縮平板邊界層

本研究首先模擬了連續(xù)流域不同馬赫數(shù)可壓縮平板邊界層流動(dòng)，以驗(yàn)證介觀方法及傳統(tǒng)CFD在連續(xù)流域模擬的一致性。計(jì)算域x×y= [0,1 m]×[0,0.5 m]，網(wǎng)格量：100×100，壁面第一層網(wǎng)格高度：0.000 1 m。來流狀態(tài)為：高度H= 55 km，靜壓p= 42.5 Pa，靜溫T=206.7 K，密度ρ= 5.68×10-4kg/m3，馬赫數(shù)Ma= 2、4、8，各 狀 態(tài) 的 雷 諾 數(shù)Re分 別 為：2.22×104、4.45×104、8.89×104。壁面溫度等于來流靜溫。分別采用N-S及GKUA進(jìn)行了Ma= 2、4、8流場的數(shù)值模擬，并將其速度場和溫度場結(jié)果與文獻(xiàn)[34]給出的可壓縮層流邊界層方程的理論結(jié)果進(jìn)行了對比，見圖2。根據(jù)邊界層相似理論，其橫坐標(biāo) η定義如下：

對比可知，采用單一方法計(jì)算連續(xù)流域可壓縮平板邊界層流動(dòng)給出的速度場及溫度場結(jié)果與文獻(xiàn)吻合良好，驗(yàn)證了本文的GKUA及N-S的可靠性。

隨著高度的增加，稀薄效應(yīng)將逐漸顯現(xiàn)，壁面附近將出現(xiàn)速度滑移及溫度跳躍。本文對來流條件：高度H= 80 km，靜壓p= 1.05 Pa，靜溫T= 198.6 K，密度ρ= 1.845 8×10-5kg/m3，馬赫數(shù)Ma= 2，雷諾數(shù)Re= 7.9×102的可壓縮平板邊界層流動(dòng)進(jìn)行了模擬，壁面溫度等于來流靜溫。數(shù)值方法包括GKUA方法、無滑移NS方法、DSMC方法及采用前文所述理論構(gòu)建的耦合方法—“Kinetic Boundary”。其中，耦合方法的本構(gòu)關(guān)系修正區(qū)域給定為壁面附近2倍平均分子自由程高度范圍內(nèi)的區(qū)域。由圖3所示壁面物理量對比可以看出，構(gòu)建的耦合方法的壁面壓力計(jì)算結(jié)果與GKUA、N-S計(jì)算結(jié)果吻合良好，略大于DSMC結(jié)果；耦合方法計(jì)算的溫度跳躍結(jié)果和GKUA吻合良好，并和DSMC基本吻合；滑移速度及壁面剪切應(yīng)力與GKUA、DSMC計(jì)算結(jié)果吻合良好。并且，由于平板前緣的速度滑移現(xiàn)象較平板后部嚴(yán)重，采用耦合方法、GKUA及DSMC方法計(jì)算的壁面剪切應(yīng)力在平板前緣較N-S??；在平板后緣，各種結(jié)果相互吻合良好。

圖2 不同馬赫數(shù)下平板邊界層模擬結(jié)果Fig. 2 Simulation results of a flat-plate boundary layer at different Mach numbers

圖3 Ma = 2平板邊界層壁面物理量Fig. 3 Distributions of physical quantities along the streamwise direction on the surface of a flat-plate boundary layer at Ma = 2

3.2 稀薄圓柱繞流

為了進(jìn)一步驗(yàn)證采用氣體動(dòng)理論耦合方法的有效性，本文以高度H= 90 km，靜壓p= 0.183 Pa，靜溫T= 186.8 K，密度ρ= 3.416 31×10-6kg/m3，馬赫數(shù)Ma=2、4，其雷諾數(shù)Re分別為150、300作為來流條件，模擬了圓柱超聲速稀薄繞流，并對流場及壁面物理量進(jìn)行了對比分析。壁面溫度等于來流靜溫。其中，耦合方法的N-S方程本構(gòu)關(guān)系修正區(qū)域給定為壁面附近2倍平均分子自由程高度范圍內(nèi)的區(qū)域。圓柱半徑1 m，遠(yuǎn)場半徑10 m，網(wǎng)格量：周向60×徑向80，壁面第一層網(wǎng)格高度0.000 5 m。

圖4 Ma = 2圓柱稀薄繞流流場Fig. 4 Rarefied flow fields around a circular cylinder at Ma = 2

圖5 圓柱稀薄繞流壁面物理量Fig. 5 Distributions of physical quantities along the circumferential direction on the surface of a circular cylinder in rarefied flow

圖4給出了氣體動(dòng)理論耦合方法計(jì)算的流場馬赫數(shù)、靜壓、密度云圖的與DSMC計(jì)算結(jié)果的對比情況。圖5給出了圓柱壁面物理量計(jì)算結(jié)果對比，其中，橫坐標(biāo)θ= 180°對應(yīng)于圓柱迎風(fēng)面駐點(diǎn)。從圖4可以看出，本文構(gòu)建的耦合方法獲得的流場與DSMC獲得的流場在圓柱迎風(fēng)側(cè)吻合良好，在圓柱背風(fēng)回流區(qū)存在一定的差異。結(jié)合圖5的壁面剪切應(yīng)力及滑移速度對比圖可看出，本文構(gòu)建的耦合方法獲得的回流區(qū)小于DSMC獲得的分離區(qū)，是導(dǎo)致圓柱背風(fēng)區(qū)流場云圖存在差異的主要原因。

從圖5可知，各方法獲得的壁面壓力吻合良好。壁面速度滑移結(jié)果顯示：DSMC、GKUA及氣體動(dòng)理論耦合方法在圓柱迎風(fēng)一側(cè)(θ∈[90°,180°])吻合良好；在背風(fēng)區(qū)，由于計(jì)算的分離區(qū)的大小差異導(dǎo)致結(jié)果存在一定偏差，但總體規(guī)律基本吻合。進(jìn)一步地，雖然圖中給出的Ma= 4流動(dòng)無量綱壁面速度滑移小于Ma= 2流動(dòng)，但其有量綱的絕對量是大于Ma2流動(dòng)的壁面速度滑移。壁面剪切應(yīng)力對比結(jié)果顯示，DSMC、GKUA及氣體動(dòng)理論耦合方法計(jì)算結(jié)果基本吻合，而無滑移NS計(jì)算結(jié)果和其他結(jié)果存在一定差異。并且，該差異整體上隨壁面滑移速度的增大而增加，流動(dòng)分離點(diǎn)(θ≈ 40°)附近的流動(dòng)較復(fù)雜，其差異最大。壁面熱流計(jì)算結(jié)果顯示無滑移N-S結(jié)果略大于其他方法計(jì)算結(jié)果，但各方法計(jì)算結(jié)果總體吻合。

4 結(jié) 論

本研究通過對多原子氣體ES-BGK方程進(jìn)行Chapman-Enskog漸近分析，構(gòu)建了分布函數(shù)、應(yīng)力張量、熱流矢量的漸近一階表達(dá)式，證明了ES-BGK方程與NS方程在連續(xù)流域的物理描述上具有理論一致性。在此基礎(chǔ)上，提出通過對稀薄效應(yīng)明顯的流動(dòng)區(qū)域進(jìn)行應(yīng)力張量及熱流矢量修正，構(gòu)造了氣體動(dòng)理論耦合方法，避免了離散速度坐標(biāo)法引入的質(zhì)量、動(dòng)量及能量守恒型誤差對耦合方法計(jì)算結(jié)果的影響。

進(jìn)一步地，本研究通過引入約化變量，構(gòu)造了分布函數(shù)一階表達(dá)式的二維約化形式，并通過有限體積法對二維耦合方程進(jìn)行數(shù)值求解。通過二維可壓縮平板邊界層及圓柱繞流算例的數(shù)值模擬及結(jié)果對比，初步驗(yàn)證了本研究提出氣體動(dòng)理論耦合方法的有效性，預(yù)期進(jìn)一步研究發(fā)展，可得到一套適于近連續(xù)滑移過渡流區(qū)高效可靠的氣體動(dòng)理論耦合方法。

致謝：感謝中國科學(xué)院力學(xué)研究所李新亮研究員團(tuán)隊(duì)提供的幫助。