基于非線(xiàn)性耦合本構(gòu)關(guān)系的改進(jìn)邊界條件

江中正，趙文文，袁震宇，陳偉芳

浙江大學(xué) 航空航天學(xué)院，杭州 310027

高超聲速稀薄流動(dòng)一直以來(lái)都是計(jì)算流體力學(xué)研究中的熱點(diǎn)與難點(diǎn)問(wèn)題。由于在遠(yuǎn)離局部平衡態(tài)流域，譬如激波間斷處、靠近物面努森層內(nèi)和尾流膨脹區(qū)，沒(méi)有足夠的分子碰撞去平衡氣體輸運(yùn)行為，氣體流動(dòng)表現(xiàn)出極強(qiáng)的非平衡特點(diǎn)。采用一階Maxwell-Smoluchowski (M-S) 滑移邊界條件的線(xiàn)性牛頓黏性/傅里葉熱傳導(dǎo)本構(gòu)模型(NSF)在這些流域中的模擬能力十分有限，無(wú)法準(zhǔn)確地刻畫(huà)非平衡流中所表現(xiàn)出的強(qiáng)非線(xiàn)性特點(diǎn)。為了準(zhǔn)確模擬高超聲速稀薄非平衡流動(dòng)，國(guó)內(nèi)外眾多學(xué)者開(kāi)展了大量研究[1]，提出了各種行之有效的理論和數(shù)值方法，包括粒子仿真方法——蒙特卡洛直接數(shù)值模擬(Direct Simulation Monte Carlo、DSMC)[2]、線(xiàn)性化玻爾茲曼方法[3]、模型方程方法(BGK、ES-BKG、Shakhov)[4]、離散速度法(Discrete Velocity Method, DVM)[5]、擴(kuò)展動(dòng)力學(xué)方法(Chapmann-Enskog)[6]、Burnett方程[7]、Grad矩方法[8]及其正則化理論[9]和全局雙曲化理論[10]。近年興起的統(tǒng)一氣體動(dòng)理論格式(Unified Gas-Kinetic Scheme, UGKS)[11]和氣體運(yùn)動(dòng)論統(tǒng)一算法(Gas Kinetic Unified Algorithm, GKUA)[12]等氣體動(dòng)理學(xué)方法也為多尺度非平衡流動(dòng)模擬提供了十分有效的手段。

為了彌補(bǔ)線(xiàn)性納維-斯托克斯-傅里葉(Navier-Stokes-Fourier, NSF)本構(gòu)模型的精度不足，同時(shí)克服粒子類(lèi)方法和玻爾茲曼求解方法計(jì)算效率較低的缺陷，近年來(lái)一套非線(xiàn)性耦合本構(gòu)關(guān)系模型(Nonlinear Coupled Constitutive Relations, NCCR)在由B.C. Eu提出的嚴(yán)格滿(mǎn)足熱力學(xué)第二定律的廣義流體動(dòng)力學(xué)方程(Generalized Hydrodynamic Equations, GHE)[13-14]的基礎(chǔ)上被發(fā)展起來(lái)。該理論自提出之后得到了深入的研究和驗(yàn)證[15-19]，在非平衡流區(qū)域表現(xiàn)出優(yōu)于Navier-Stokes方程的特點(diǎn)。然而，針對(duì)該理論的壁面邊界條件研究為數(shù)不多，通常做法是直接采用一階M-S邊界條件開(kāi)展計(jì)算[20-21]。在一階M-S邊界條件中，滑移速度或跳躍溫度是速度或溫度一階梯度的線(xiàn)性比值，因此難以將NCCR模型的非線(xiàn)性特征表現(xiàn)出來(lái)，最終導(dǎo)致在壁面處人為降低了模型精度。

本文在Maxwell散射基礎(chǔ)上開(kāi)展了滑移邊界理論研究，將努森層滑移跳躍效應(yīng)通過(guò)與高階非守恒量(應(yīng)力和熱流)表征的非平衡程度聯(lián)系起來(lái)，提出一套與NCCR模型本身精度相適應(yīng)的高階邊界條件，同時(shí)針對(duì)二維圓柱單原子氣體繞流和平板繞流問(wèn)題給出了不同稀薄程度下不同滑移邊界條件NSF與NCCR模型的模擬結(jié)果，通過(guò)與DSMC的壁面結(jié)果進(jìn)行比較，驗(yàn)證該邊界條件模擬物面附近非平衡流動(dòng)的能力。

1 非線(xiàn)性耦合本構(gòu)理論

1.1 廣義流體動(dòng)力學(xué)方程

考慮外力影響的針對(duì)稀疏單組分單原子氣體的玻爾茲曼方程表示為

(1)

式中：f表示包含時(shí)間、空間、速度七維相空間的速度分布函數(shù)；F為體積力；碰撞積分表示為

(2)

其中：g=|v-v1|表示粒子的相對(duì)速度；b為粒子碰撞參數(shù)，亦稱(chēng)為瞄準(zhǔn)距離；φ為碰撞平面的方位角；v為分子運(yùn)動(dòng)速度。

為了推導(dǎo)廣義流體動(dòng)力學(xué)方程，在這里首先定義一個(gè)廣義速度矩：

Φ(abc…l)=〈mCaCbCc…Clf〉

(3)

式中：角括號(hào)代表對(duì)分子速度的積分；C為分子熱運(yùn)動(dòng)速度。對(duì)式(3)進(jìn)行時(shí)間求導(dǎo)，得到廣義矩的時(shí)間變化率，然后將玻爾茲曼方程式(1)代入，并定義碰撞耗散項(xiàng)為

Λ(Φ)(abc…kl)=〈mCaCbCc…ClC(f,f1)〉

(4)

便得到廣義速度矩的演化方程為

(5)

若將關(guān)心的宏觀物理量分為守恒量(ρ,u,E)和非守恒量(Π,Q)，均可以通過(guò)統(tǒng)一的微觀式來(lái)表達(dá)，即

Φ(k)=〈h(k)f〉

(6)

式中：h(k)為兩組物理量的分子微觀表達(dá)式；宏觀物理量(密度、動(dòng)量、內(nèi)能、應(yīng)力張量和熱流)分別表示為

(7)

其對(duì)應(yīng)的分子微觀表達(dá)式為

(8)

其中：m為分子質(zhì)量；數(shù)學(xué)符號(hào)[A](2)表示二階無(wú)跡對(duì)稱(chēng)張量，即

(9)

由于守恒量是碰撞不變量，不存在碰撞耗散項(xiàng)，因此由廣義速度矩的演化方程式(5)在不考慮體積力影響下可以得到3大守恒方程為

(10)

(11)

(12)

為了推導(dǎo)非守恒量的演化方程以及封閉非守恒系統(tǒng)的分子碰撞耗散項(xiàng)式(4)，Eu[13, 22]構(gòu)建了非平衡態(tài)分布函數(shù)的形態(tài)定義，在數(shù)學(xué)和物理的層面保證熵和Calortropy變化的非負(fù)屬性。這個(gè)非平衡態(tài)分布函數(shù)充當(dāng)了連接熵增源項(xiàng)和宏觀非守恒量耗散演化過(guò)程之間的橋梁。該分布函數(shù)的單組分單原子氣體分子的形態(tài)定義為

(13)

式中：μ為標(biāo)準(zhǔn)化因子，用以平衡氣體分子數(shù)密度n與該模型分布函數(shù)的關(guān)系，如

其中：Xk為與宏觀量有關(guān)的系數(shù)，其作用跟Grad矩方程中的系數(shù)相似；T、kB分別為溫度、Boltzmann常數(shù)。利用Eu分布函數(shù)形態(tài)定義的形式計(jì)算熵增項(xiàng)，可以得到熵增項(xiàng)的表達(dá)式為

(14)

Eu認(rèn)為能量耗散是系統(tǒng)的分子間碰撞導(dǎo)致的，微觀上的分子碰撞又會(huì)導(dǎo)致非守恒量耗散演化的宏觀過(guò)程，而熵增特性正是對(duì)能量耗散的一種直接體現(xiàn)。式(14)建立了耗散項(xiàng)、宏觀非守恒量和熵增項(xiàng)的聯(lián)系，這正是Eu修正矩方法的基石，也是和Grad矩方法有所區(qū)別的關(guān)鍵。通過(guò)對(duì)熵增項(xiàng)式(14)進(jìn)行累積量展開(kāi)[23]，并截取一階近似，可得

(15)

當(dāng)趨于平衡態(tài)流動(dòng)時(shí)，κ趨于在近平衡態(tài)發(fā)展的Rayleigh-Onsager耗散函數(shù)，即

(16)

式中：d、p、η、λ分別為分子直徑、壓力、氣體黏性系數(shù)和熱傳導(dǎo)系數(shù)。由式(16)可見(jiàn)Eu的熵增項(xiàng)累積量式(15)滿(mǎn)足熵非負(fù)屬性。在此基礎(chǔ)上通過(guò)式(14)求得碰撞耗散項(xiàng)，最終可以得到一組針對(duì)單組分單原子氣體分子的非守恒量演化方程為

(17)

式中：ψ(4)和ψ(5)為更高階矩；cp為定壓比熱。式(17)結(jié)合3大守恒方程(式(10)～式(12))，即為Eu的廣義流體動(dòng)力學(xué)方程。

1.2 非線(xiàn)性耦合本構(gòu)關(guān)系

Eu的非守恒量演化方程仍舊是開(kāi)放的包含了更高階矩的偏微分方程系統(tǒng)。為了封閉該方程組，Eu和Alghoul[24]提出了另一種不同于Grad封閉的基本封閉原則：

ψ(4)=ψ(5)=0

(18)

(19)

為方便對(duì)比，這里給出線(xiàn)性的牛頓黏性和傅里葉傳熱導(dǎo)定律：

(20)

2 邊界條件

在高超聲速稀薄流動(dòng)模擬中，氣固表面分子碰撞與作用機(jī)理是一個(gè)十分重要而且復(fù)雜的影響因素。入射氣體分子在物體表面附近與物面反射分子相互碰撞聚集，這個(gè)聚集層通常被稱(chēng)為努森層，其厚度與平均分子自由程l是同一個(gè)量級(jí)。當(dāng)來(lái)流并不稀薄，即來(lái)流平均分子自由程l→0時(shí)，努森層厚度趨于零，在微觀分子層面表征為入射分子被物體表面吸附，在宏觀層面即為最常用的附著壁和等溫壁條件。然而，當(dāng)來(lái)流稀薄程度增大，物面努森層效應(yīng)逐漸顯著，附著壁和等溫壁邊界條件不再成立時(shí)，氣體流動(dòng)與壁面之間會(huì)呈現(xiàn)不連續(xù)的現(xiàn)象，即引起速度滑移和溫度跳躍現(xiàn)象。氣固表面分子作用機(jī)理在準(zhǔn)確刻畫(huà)努森層流動(dòng)乃至整個(gè)宏觀稀薄流動(dòng)中具有十分重要的作用。為此，有研究提出兩種常見(jiàn)的氣固分子碰撞模型：Maxwell散射理論和Langmuir表面吸附理論。作者前期研究工作[19]表明，盡管Langmuir吸附理論為氣固碰撞提供了一種吸附等溫物理模型，但由此衍生出的Langmuir滑移邊界[28]并不能很好地預(yù)測(cè)稀薄流壁面處的滑移速度和跳躍溫度。因此，本文主要基于Maxwell散射理論開(kāi)展邊界條件研究。該理論中Maxwell首先引入兩個(gè)簡(jiǎn)單假設(shè)模型：鏡面反射和完全漫反射。在來(lái)流反射分子中σ部分假定為完全漫反射，(1-σ)部分為鏡面反射，σ為介于0～1之間的比例系數(shù)。

常見(jiàn)的一階Maxwell滑移條件有3種形式：Maxwell滑移條件、Gokcen滑移條件、Lockerby滑移條件。Maxwell[29]在1879年基于分子運(yùn)動(dòng)論和小努森數(shù)假設(shè)，推導(dǎo)了一階泰勒展開(kāi)的滑移條件為

(21)

式中：u為靠近物面處氣體切向速度；uw為物面速度；σu為速度適應(yīng)系數(shù)；y為壁面法向；x為壁面切向；右端第2項(xiàng)則是在非等溫壁下的熱蠕動(dòng)效應(yīng)。

對(duì)于溫度跳躍效應(yīng)，Smoluchowski[30]在1898年也推出了相類(lèi)似的溫度邊界條件：

(22)

式中：T為靠近物面處氣體溫度；Tw為壁溫；σT為溫度適應(yīng)系數(shù)；γ為比熱比；Pr為普朗特?cái)?shù)。式(21)和式(22)也被稱(chēng)為M-S邊界條件。M-S邊界條件結(jié)合Navier-Stokes方程，既有效地增強(qiáng)了Navier-Stokes方程模擬滑移稀薄流的能力又保證了較小的計(jì)算量，在工程上得到廣泛應(yīng)用。

然而，G?k?en和MacCormack[31]研究認(rèn)為基于小努森數(shù)一階展開(kāi)得到的M-S邊界條件模擬努森數(shù)范圍十分有限，在大的局部努森數(shù)區(qū)域預(yù)測(cè)的壁面結(jié)果與粒子仿真方法的結(jié)果有很大差異。因此，他們通過(guò)研究壁面切向動(dòng)量損失和引入努森層處的速度，提出另一種通用的滑移條件:

(23)

式中:a代表速度或者溫度；下標(biāo)l表示一個(gè)平均分子自由程處的物理量，w表示物面處的物理量。該條件被認(rèn)為在小努森數(shù)時(shí)能夠還原Maxwell條件，在較大努森數(shù)條件下能得到更優(yōu)的壁面結(jié)果。但是等式右端依然用到線(xiàn)性的NSF本構(gòu)關(guān)系(應(yīng)力張量隨應(yīng)變張量線(xiàn)性變化)，實(shí)際上這與努森層內(nèi)流動(dòng)速度分布的非線(xiàn)性特點(diǎn)不符。因此G?k?en邊界條件的可用范圍值得商榷。

此外，Lockerby[32]針對(duì)切應(yīng)力均勻分布平板壁面的Cercignani線(xiàn)性玻爾茲曼速度解進(jìn)行曲線(xiàn)擬合，通過(guò)引入一個(gè)壁面函數(shù)式(24)來(lái)修正努森層黏性ηnew=ηΨ-1的方式來(lái)修正線(xiàn)性的速度梯度項(xiàng)。

(24)

同時(shí)Lockerby等的研究[33]回歸了不考慮熱蠕動(dòng)Maxwell邊界條件的一般形式，認(rèn)為右端項(xiàng)采用應(yīng)力形式比應(yīng)變率形式更具有普適性。

(25)

(26)

在Maxwell一般形式(式(25)和式(26))中，采用由壁面函數(shù)修正的有效黏性ηnew=ηΨ-1以及線(xiàn)性的牛頓黏性定律和傅里葉熱傳導(dǎo)定律，便可得到Lockerby邊界條件。Lockerby通過(guò)進(jìn)一步研究認(rèn)為式(25)前面仍需要乘上系數(shù)0.8。該邊界條件中通過(guò)構(gòu)造非線(xiàn)性函數(shù)去修正線(xiàn)性的應(yīng)力應(yīng)變關(guān)系的做法，為嘗試捕捉努森層的非線(xiàn)性物理量分布特點(diǎn)提供了一種新的思路。

此外，精確描述速度滑移和溫度跳躍現(xiàn)象的另一條途徑是構(gòu)造二階乃至更高階的邊界條件，Cercignani[34]、Beskok[35]、Deissler[36]、Hsia[37]等在Maxwell一階邊界條件基礎(chǔ)上做了富有成效的研究。尤其Deissler通過(guò)物理推導(dǎo)手段以及Hsia和Domoto通過(guò)實(shí)驗(yàn)手段，均指出了單純通過(guò)對(duì)Maxwell條件進(jìn)行泰勒展開(kāi)的數(shù)學(xué)方法來(lái)構(gòu)造高階邊界條件會(huì)存在二階乃至偶數(shù)階滑移項(xiàng)反向過(guò)度修正的問(wèn)題。為了區(qū)別傳統(tǒng)的二階M-S邊界條件，現(xiàn)將偶數(shù)項(xiàng)修正之后的Hsia-Domoto邊界條件表示為

(27)

(28)

對(duì)比Maxwell邊界條件一般形式(式(25)和式(26))和不考慮熱蠕動(dòng)的常見(jiàn)形式(式(21)和式(22))，其主要的差異在線(xiàn)性的NSF本構(gòu)模型(式(20))。參考Lockerby利用有效黏性修正真實(shí)氣體黏性的非線(xiàn)性做法，對(duì)Maxwell邊界一般形式采取另一種非線(xiàn)性修正方法：在保留真實(shí)氣體黏性系數(shù)前提下結(jié)合NCCR本構(gòu)模型，利用NCCR的應(yīng)力和熱流修正NSF的應(yīng)力與熱流項(xiàng)。考慮到前期研究工作[19]定量指出二階修正項(xiàng)在準(zhǔn)確預(yù)測(cè)努森層滑移邊界值的重要性，在提出的修正邊界條件也引入了二階修正項(xiàng)(式(27)和式(28)) 的影響，并忽略熱蠕動(dòng)項(xiàng)，具體表示為

(29)

(30)

式(29)和式(30)中的應(yīng)力和熱流通過(guò)在壁面求解NCCR方程得到。此外，為了防止滑移速度和跳躍溫度變化劇烈而導(dǎo)致高階邊界條件在迭代計(jì)算過(guò)程中的不穩(wěn)定，Lockerby[38]和包福兵[39]在壁面處采用邊界松弛方法進(jìn)行計(jì)算，并取得較好的計(jì)算穩(wěn)定性。因此，本文也采取相同的邊界松弛方式:

(31)

式中：us和Tj分別為滑移速度和跳躍溫度；上標(biāo)r、n+1、n分別表示松弛之后、新時(shí)刻、上一時(shí)刻的值；Rf為迭代松弛因子，根據(jù)具體算例選取不同的值。

3 數(shù)值算法

選取無(wú)窮遠(yuǎn)來(lái)流參數(shù)對(duì)3大守恒方程和非線(xiàn)性耦合本構(gòu)模型(NCCR)進(jìn)行無(wú)量綱化，即

(32)

式中:下標(biāo)∞表示無(wú)窮遠(yuǎn)來(lái)流參數(shù)。為了簡(jiǎn)化表達(dá)，下文均忽略上標(biāo)星號(hào)，最終得到無(wú)量綱化的流動(dòng)控制方程為

(33)

式中：守恒量U、無(wú)黏通量Fc和黏性通量Fv分別為

(34)

針對(duì)單原子氣體的非線(xiàn)性耦合本構(gòu)關(guān)系，有

(35)

式中：E為比總能；下標(biāo)0代表線(xiàn)性的NSF應(yīng)力和熱流，上標(biāo)^表示參數(shù)有如下關(guān)系：

(36)

其中：γ代表參考值。

3.1 空間離散

在有限體積法框架下離散求解流動(dòng)控制方程式(33)。為了克服AUSM+的數(shù)值振蕩和AUSMD的“紅玉”現(xiàn)象，Kim等[40]結(jié)合這兩種格式的優(yōu)勢(shì)并且考慮到計(jì)算效率和魯棒性，通過(guò)引入修正的壓力權(quán)函數(shù)，提出了AUSMPW+通量計(jì)算格式。該格式計(jì)算效率較高，魯棒性好，而且激波捕捉能力強(qiáng)，在高速流動(dòng)中得到了廣泛應(yīng)用，因此在本文的無(wú)黏通量計(jì)算中也采用AUSMPW+格式。AUSMPW+格式的通量表達(dá)式為

(37)

(38)

式中：

(Δ-)i=qi-qi-1

(Δ+)i=qi+1-qi

其中重構(gòu)主要針對(duì)原始變量q進(jìn)行，ε是防止分母為零的一個(gè)小量，而黏性通量均采用中心差分離散。

3.2 時(shí)間推進(jìn)

本文采用LU-SGS (Lower-Upper Symmetric Gauss-Seidel)隱式時(shí)間推進(jìn)方法，該方法具有較好的計(jì)算穩(wěn)定性和收斂性，是最常用的隱式時(shí)間推進(jìn)方法之一。對(duì)控制方程式(33)進(jìn)行時(shí)間一階隱式離散后可得

[Ι+Δt(DξA+DηB+DζC)]ΔQn=-ΔtRHS

(39)

式中：Dξ、Dη、Dζ為微分算子；RHS為殘差項(xiàng)；A、B和C為無(wú)黏雅克比矩陣。通過(guò)矩陣分裂和LU分解得到LU-SGS方法的表達(dá)式為

LD-1UΔQ=-RHS

(40)

式中：

(41)

其中：A±、B±和C±為近似雅克比矩陣；ρ(A)、ρ(B)和ρ(C)為無(wú)黏雅克比矩陣的譜半徑。式(40)可通過(guò)向前掃描、向后掃描和更新3個(gè)步驟進(jìn)行求解。一般在式(39)中對(duì)黏性通量進(jìn)行顯式時(shí)間離散，但在黏性影響較大的區(qū)域容易造成計(jì)算不穩(wěn)定。參考Navier-Stokes方程求解方法，本文也通過(guò)修改近似雅克比矩陣特征值，采用譜半徑代替黏性雅克比矩陣，對(duì)NCCR黏性項(xiàng)進(jìn)行了隱式處理，以ξ方向?yàn)槔?

(42)

3.3 本構(gòu)模型求解

NCCR本構(gòu)模型式(35)是應(yīng)力和熱流的隱式表達(dá)，在3個(gè)方向上各項(xiàng)相互耦合影響，并不能像NSF本構(gòu)關(guān)系簡(jiǎn)單地通過(guò)守恒量梯度一步顯式求解，需要額外求解步驟。為此，Myong[17-18]提出的分裂算法通過(guò)忽略一些耦合項(xiàng)，將流動(dòng)分解成3個(gè)方向之間互不影響的一維問(wèn)題，然后進(jìn)行線(xiàn)性疊加。前期研究[20-21]表明，分裂算法容易造成計(jì)算不穩(wěn)定，尤其在尾流的模擬中表現(xiàn)明顯，同時(shí)分裂算法忽略的一些耦合項(xiàng)也人為降低了模型的計(jì)算精度。因此，本文采用完整耦合的求解思路，利用最速下降法對(duì)NCCR本構(gòu)模型進(jìn)行求解，有效保證模型的計(jì)算精度。最速下降法不過(guò)分依賴(lài)初值的選取，是求解單組分單原子NCCR模型的有效算法。對(duì)于NCCR模型，可以表示成如下非線(xiàn)性方程組的一般形式：

F(x)=0

(43)

(44)

(45)

如何選取合適的步長(zhǎng)αk，使得

(46)

是最速下降法的關(guān)鍵。最直接的解析方法是構(gòu)造有關(guān)步長(zhǎng)α的單值函數(shù)，即

(47)

通過(guò)對(duì)α求導(dǎo)解一元方程的根來(lái)尋找函數(shù)極小值。但是這種做法計(jì)算成本太高，本文采用一維搜索近似得到每步較優(yōu)步長(zhǎng)來(lái)進(jìn)行迭代。選取3個(gè)步長(zhǎng)α1=0<α2=α3/2<α3，使得h(α3)

(48)

4 計(jì)算結(jié)果與分析

4.1 圓柱繞流

首先研究一組在不同稀薄來(lái)流程度下的高超聲速單原子氣體圓柱繞流問(wèn)題，該算例的DSMC數(shù)據(jù)以及Gokcen和Lockerby邊界條件的NSF數(shù)據(jù)均來(lái)自文獻(xiàn)[41]。圓柱直徑為12英寸，如圖1所示。無(wú)窮遠(yuǎn)來(lái)流速度為U∞=2 624.1 m/s，馬赫數(shù)為10。采用等溫壁模型，壁溫為500 K。不同稀薄來(lái)流下的氣體密度由表1給出，其中Kn∞表示采用無(wú)窮遠(yuǎn)來(lái)流條件計(jì)算的克努森數(shù)。

圖1 幾何外形定義Fig.1 Definition of geometry

Kn∞Mass density/(kg·m-3)Number density/(particles·m-3)0.0021.408×10-42.124×10210.012.818×10-54.247×10200.055.636×10-68.494×10190.251.127×10-61.699×1019

單原子氣體為氬氣，比熱比為1.67，氣體常數(shù)為208.16，普朗特?cái)?shù)為0.67，其黏性計(jì)算采用變徑硬球(Variable Hard Sphere, VHS)模型：

(49)

(50)

式中：黏性指數(shù)s為0.734；參考分子直徑dr為3.595×10-10m；參考溫度Tr為1 000 K；kB取1.380 658×10-23J·K-1；m為氬氣的分子質(zhì)量。

本文研究重點(diǎn)是量化比較線(xiàn)性(NSF)與非線(xiàn)性(NCCR)兩個(gè)宏觀本構(gòu)模型在幾組不同邊界條件影響下預(yù)測(cè)的圓柱物面壓力、摩阻和熱流系數(shù)與DSMC結(jié)果之間的差異。因此物面無(wú)量綱壓力系數(shù)Cp、摩阻系數(shù)Cf和熱流系數(shù)Ch定義為

(51)

圖2～圖5分別給出了Kn=0.002,0.01,0.05,0.25 這4組條件下，壓力系數(shù)Cp、摩阻系數(shù)Cf和熱流系數(shù)Ch沿圓柱壁面隨著幾何角度θ的變化情況，并結(jié)合一階M-S、Gokcen、Lockerby邊界的NSF預(yù)測(cè)結(jié)果以及結(jié)合一階M-S和二階非線(xiàn)性修正邊界的NCCR結(jié)果，與DSMC粒子仿真結(jié)果(由MONACO軟件計(jì)算得到)進(jìn)行比較。來(lái)流Kn=0.002的流動(dòng)處于連續(xù)流域范疇。

圖2 物面系數(shù)曲線(xiàn)(Kn=0.002)Fig.2 Curves of surface coefficient (Kn=0.002)

由圖2可以看出，線(xiàn)性NSF模型和非線(xiàn)性NCCR模型的預(yù)測(cè)壁面壓力、摩阻、熱流系數(shù)結(jié)果與DSMC的結(jié)果非常吻合?？紤]到計(jì)算成本的因素，線(xiàn)性NSF模型仍然是連續(xù)流域模擬中既有效又高效的理論方法。而從另一方面看，本文基于NCCR模型提出的二階修正邊界條件在連續(xù)流域能有效回歸NSF的預(yù)測(cè)結(jié)果，體現(xiàn)了其在近平衡流附近的模擬能力。

圖3 物面系數(shù)曲線(xiàn)(Kn=0.01)Fig.3 Curves of surface coefficient (Kn=0.01)

由圖3可看出，當(dāng)Kn=0.01時(shí)，除了在θ=120°尾跡區(qū)附近的物面，結(jié)合一階M-S邊界NSF預(yù)測(cè)的摩阻系數(shù)稍微偏高之外，其余位置不同方法的計(jì)算結(jié)果依然吻合較好。Gokcen和Lockerby邊界比一階M-S邊界有效提升了NSF的摩阻預(yù)測(cè)能力。相比NSF本構(gòu)模型，結(jié)合一階M-S邊界和二階非線(xiàn)性修正邊界的NCCR都準(zhǔn)確模擬3組物面系數(shù)。來(lái)流Kn=0.01可以被認(rèn)為是無(wú)滑移流域的上限，是連續(xù)流和稀薄滑移流域的界限。盡管流動(dòng)在物面附近開(kāi)始出現(xiàn)連續(xù)性假設(shè)失效，預(yù)測(cè)的結(jié)果表明滑移邊界條件仍能夠有效拓展宏觀本構(gòu)模型的非平衡流動(dòng)模擬能力。

如圖4所示，隨著來(lái)流努森數(shù)繼續(xù)增加，不同方法壁面壓力系數(shù)的預(yù)測(cè)依然較為一致，但摩阻和熱流系數(shù)預(yù)測(cè)的結(jié)果出現(xiàn)明顯的差異。盡管來(lái)流Kn=0.05的流動(dòng)處于滑移流域，但尾跡壁面附近區(qū)域的非平衡效應(yīng)更加強(qiáng)烈，連續(xù)性失效更加嚴(yán)重。對(duì)于結(jié)合一階M-S邊界和Lockerby邊界的NSF結(jié)果，尾跡區(qū)的物面摩阻和熱流系數(shù)預(yù)測(cè)偏高，而Gokcen邊界相對(duì)更為準(zhǔn)確。與所有NSF結(jié)果相比，NCCR預(yù)測(cè)結(jié)果更加趨近于DSMC仿真結(jié)果，這表明非線(xiàn)性耦合本構(gòu)模型(NCCR)比線(xiàn)性的牛頓黏性和傅里葉熱傳導(dǎo)模型(NSF)對(duì)非平衡流動(dòng)模擬精度更高，能為有效解決連續(xù)性失效流域的模擬問(wèn)題提供更加可靠的理論手段。

隨著努森數(shù)進(jìn)一步增加，流動(dòng)離開(kāi)滑移流域進(jìn)入過(guò)渡流域。在來(lái)流Kn=0.25的流動(dòng)中氣體已經(jīng)非常稀薄，由于缺乏足夠的氣體粒子碰撞，非平衡效應(yīng)非常強(qiáng)烈。圖5的物面系數(shù)分布表明，在過(guò)渡流域模擬中，滑移邊界不再能有效地提升線(xiàn)性NSF本構(gòu)模型的預(yù)測(cè)能力。強(qiáng)烈的非平衡效應(yīng)不僅存在壁面努森層流動(dòng)中，而且?guī)缀鯏U(kuò)展到整個(gè)流場(chǎng)范圍。在過(guò)渡流域中，若遠(yuǎn)離物面流動(dòng)仍然使用線(xiàn)性的NSF本構(gòu)模型進(jìn)行模擬，將嚴(yán)重影響到物面系數(shù)的預(yù)測(cè)精度。NSF的壁面壓力和熱流系數(shù)均比DSMC的預(yù)測(cè)結(jié)果偏高，除了采用Gokcen邊界的NSF摩阻系數(shù)偏低外，其他滑移邊界的摩阻系數(shù)也普遍偏高，壁面摩阻峰值偏高約12.5%。相比NSF的預(yù)測(cè)結(jié)果，NCCR在壁面系數(shù)的預(yù)測(cè)中表現(xiàn)更優(yōu)，比如圓柱前部的壓力系數(shù)和后部的熱流系數(shù)更加趨近DSMC的結(jié)果。對(duì)于壁面摩阻系數(shù)的預(yù)測(cè)，一階M-S和非線(xiàn)性修正邊界的差異在過(guò)渡流域明顯開(kāi)始增大。隨著氣流加速流過(guò)圓柱，采用一階M-S滑移邊界的NCCR物面系數(shù)預(yù)測(cè)結(jié)果明顯偏離DSMC，在θ=100°處(尾流膨脹區(qū))偏離值達(dá)到最大，而結(jié)合非線(xiàn)性修正滑移邊界的NCCR摩阻預(yù)測(cè)在迎風(fēng)與背風(fēng)區(qū)域與DSMC的結(jié)果基本吻合?？梢?jiàn)，基于NCCR的非線(xiàn)性修正滑移邊界在壁面處保留非線(xiàn)性?xún)?yōu)勢(shì)，保證了NCCR理論模型的精度一致性，相比一階M-S滑移邊界極大提升了NCCR壁面系數(shù)的預(yù)測(cè)能力。

圖4 物面系數(shù)曲線(xiàn)(Kn=0.05)Fig.4 Curves of surface coefficient (Kn=0.05)

圖5 物面系數(shù)曲線(xiàn)(Kn=0.25)Fig.5 Curves of surface coefficient (Kn=0.25)

4.2 平板繞流

平板阻力預(yù)測(cè)是經(jīng)典的流體力學(xué)問(wèn)題，其外形簡(jiǎn)單，流動(dòng)信息豐富，是組成復(fù)雜飛行器結(jié)構(gòu)最基本的幾何元素，因此在低速不可壓、高速可壓縮等領(lǐng)域得到廣泛深入的研究[42-45]。對(duì)于高超聲速平板繞流問(wèn)題，激波與邊界層干擾嚴(yán)重，黏性作用加劇，在平板前緣區(qū)域存在非常強(qiáng)烈的稀薄非平衡效應(yīng)，這些疊加的復(fù)雜效應(yīng)都會(huì)影響表面摩阻的準(zhǔn)確預(yù)測(cè)，因此也是一個(gè)邊界條件研究的經(jīng)典驗(yàn)證算例。

本節(jié)以零厚度零迎角平板高超聲速繞流進(jìn)行數(shù)值模擬研究，進(jìn)一步分析結(jié)合非線(xiàn)性修正邊界條件的NCCR模型在復(fù)雜流動(dòng)的預(yù)測(cè)能力。流域左邊界和上邊界是無(wú)窮遠(yuǎn)來(lái)流，來(lái)流溫度為300 K，右邊界為出口，下邊界為X/L=0～1的零厚度平板，平板壁溫為300 K，其余為對(duì)稱(chēng)邊界。為了減少雙原子氣體在高速流動(dòng)中轉(zhuǎn)動(dòng)自由度激發(fā)帶來(lái)的影響，該算例依然采用單原子氬氣，氣體屬性均采用4.1節(jié)的描述，重點(diǎn)研究了Ma=10、Kn=0.002和Ma=5、Kn=0.01兩組平板繞流算例。壁面采用全漫反射條件。

圖6和圖7分別展示了兩組來(lái)流條件下(Ma=10、Kn=0.002和Ma=5、Kn=0.01)由NSF本構(gòu)模型結(jié)合一階M-S滑移邊界和NCCR模型結(jié)合非線(xiàn)性修正邊界條件模擬的馬赫數(shù)Ma、局部克努森數(shù)KnGLL、無(wú)量綱壓力p/p∞和密度ρ/ρ∞云圖情況。局部克努森數(shù)是由Boyd等[46]提出，用來(lái)具體刻畫(huà)當(dāng)?shù)亓饔虻倪B續(xù)性失效情況，即

(52)

該參數(shù)的方向?qū)?shù)采取當(dāng)?shù)靥荻确较颍琎可以根據(jù)具體的情況選取密度、溫度或者壓力。Boyd等[46]認(rèn)為0.05可以作為當(dāng)?shù)剡B續(xù)性失效的判據(jù)。本文基于當(dāng)?shù)孛芏忍荻葋?lái)求局部克努森數(shù)。由圖6和圖7中兩個(gè)算例的局部克努森數(shù)云圖可以看出，整個(gè)流域除了平板前緣、激波和壁面附近外，失效區(qū)域和程度都不大，因此由NSF和NCCR模擬的流場(chǎng)大部分都相似，這也可以通過(guò)圖中的馬赫數(shù)、壓力和密度云圖的對(duì)比看出。當(dāng)流動(dòng)的連續(xù)性失效程度不大，稀薄非平衡效應(yīng)不明顯時(shí)，結(jié)合非線(xiàn)性修正滑移邊界條件的NCCR能夠回歸傳統(tǒng)的NSF解。從圖6和圖7中的馬赫數(shù)云圖還能夠非常清晰地分辨在平板前緣之后的激波與邊界層融合層以及之后緊鄰的高超聲速黏性干擾區(qū)，在Ma=5、Kn=0.01的算例還能觀察到在尾部激波與邊界層之間的無(wú)黏區(qū)，而在Ma=10、Kn=0.002的算例中，由于馬赫數(shù)增大，激波角減少，邊界層增大，其無(wú)黏區(qū)范圍不明顯。在融合層，激波尚未充分發(fā)展，而邊界層相對(duì)較厚，此時(shí)激波與邊界層相互融合，難以區(qū)分，這里也會(huì)存在比較顯著的速度滑移和溫度跳躍的稀薄現(xiàn)象；而之后的高超聲速黏性干擾區(qū)，激波已經(jīng)充分發(fā)展，尤其在Ma=10時(shí)非?？拷吔鐚油饩?，邊界層則迅速增長(zhǎng)，黏性效應(yīng)增強(qiáng)，與激波相互影響著波后的無(wú)黏區(qū)和平板表面的熱流和摩阻分布。

圖6 馬赫數(shù)、局部克努森數(shù)、無(wú)量綱壓力和密度云圖對(duì)比(Ma=10、Kn=0.002)Fig.6 Comparison of detailed contour of Mach number, local Knudsen number, dimensionless pressure and density (Ma=10，Kn=0.002)

圖7 馬赫數(shù)、局部克努森數(shù)、無(wú)量綱壓力和密度云圖對(duì)比(Ma=5、Kn=0.01)Fig.7 Comparison of detailed contour of Mach number, local Knudsen number, dimensionless pressure and density (Ma=5,Kn=0.01)

圖8 表面摩擦阻力系數(shù)曲線(xiàn)(Ma=10、Kn=0.002)Fig.8 Curves of surface friction coefficient (Ma=10，Kn=0.002)

圖9 表面摩擦阻力系數(shù)曲線(xiàn)(Ma=5、Kn=0.01)Fig.9 Curves of surface friction coefficient (Ma=5,Kn=0.01)

準(zhǔn)確預(yù)測(cè)超聲速流動(dòng)的摩阻一直是氣動(dòng)研究的重難點(diǎn)。為了定量分析基于NCCR的非線(xiàn)性修正邊界的能力，本文具體比較了由DSMC、結(jié)合一階M-S的NSF、結(jié)合一階M-S和非線(xiàn)性修正邊界的NCCR預(yù)測(cè)平板單側(cè)表面摩阻系數(shù)分布的情況，如圖8和圖9所示。DSMC的結(jié)果采自文獻(xiàn)[44]。在Ma=10、Kn=0.002的X/L≥0.05和Ma=5、Kn=0.01的X/L≥0.2平板范圍，NSF、NCCR和DSMC三者預(yù)測(cè)的摩阻非常吻合，說(shuō)明這些區(qū)域仍能滿(mǎn)足連續(xù)性介質(zhì)假設(shè)，稀薄效應(yīng)不明顯。而在這些范圍之前平板前緣和融合層區(qū)域，三者預(yù)測(cè)的結(jié)果產(chǎn)生顯著的差異，NSF即使結(jié)合一階M-S邊界也無(wú)法準(zhǔn)確預(yù)測(cè)該處的摩阻，過(guò)高預(yù)估了該值。相比結(jié)合一階M-S的NCCR解，NCCR結(jié)合非線(xiàn)性修正滑移邊界條件對(duì)壁面摩阻分布和極值大小的預(yù)估精度有著比較明顯的提升，更加吻合DSMC的結(jié)果。從圖6和圖7的局部克努森數(shù)云圖可以發(fā)現(xiàn)，該處平板能夠達(dá)到KnGLL=0.5以上，處于過(guò)渡流域，稀薄非平衡效應(yīng)顯著，因此結(jié)合滑移邊界的NSF能力明顯不足，無(wú)法勝任對(duì)于該處模擬的需求，需要精度更高的理論模型結(jié)合高精度的邊界條件進(jìn)行預(yù)測(cè)。另外，在該處的摩阻分布，還能觀察到一個(gè)有趣的現(xiàn)象，表面摩擦阻力系數(shù)有個(gè)先增后降的過(guò)程，即在平板某處產(chǎn)生一個(gè)摩阻峰值。這個(gè)現(xiàn)象實(shí)際上在一定程度反映了平板前緣的非平衡效應(yīng)。當(dāng)?shù)乜伺瓟?shù)較大，壁面速度滑移明顯，因此會(huì)引起壁面切向速度的法向梯度的減少，從而降低了表面摩阻系數(shù)，這就是造成在X/L=0處摩阻系數(shù)非常小的原因。

5 結(jié) 論

本文針對(duì)最新發(fā)展的非平衡氣體動(dòng)力學(xué)理論——非線(xiàn)性耦合本構(gòu)關(guān)系模型(NCCR)開(kāi)展了邊界條件研究，提出了一套在壁面處與NCCR模型精度保持一致的非線(xiàn)性修正滑移邊界條件。在有限體積的計(jì)算框架下，結(jié)合LU-SGS隱式時(shí)間技術(shù)和AUSMPW+現(xiàn)代通量計(jì)算格式，采用耦合求解NCCR模型的思路，對(duì)不同稀薄程度的高超聲速單原子氣體圓柱繞流和平板繞流問(wèn)題進(jìn)行數(shù)值模擬，重點(diǎn)比較了采用不同滑移邊界的NSF和NCCR模型對(duì)物面壓力、摩阻與熱流系數(shù)的預(yù)測(cè)情況。

1) 結(jié)合滑移邊界條件的NSF只能有效拓展到滑移流域的模擬，NCCR模型可準(zhǔn)確模擬的非平衡流域范圍比NSF更寬。

2) 結(jié)合一階M-S邊界的NCCR模型在非平衡效應(yīng)強(qiáng)烈的努森層流動(dòng)以及激波與邊界層的平板融合干擾區(qū)的模擬中會(huì)弱化NCCR的非線(xiàn)性?xún)?yōu)勢(shì)，摩阻預(yù)測(cè)明顯偏離DSMC仿真結(jié)果。而本文的非線(xiàn)性修正邊界條件有效增強(qiáng)了NCCR模型在過(guò)渡流域和激波/邊界層融合干擾區(qū)的非線(xiàn)性模擬能力，在壁面流動(dòng)中保證了與流動(dòng)模型精度的一致性。

3) 在稀薄效應(yīng)影響下，NSF預(yù)測(cè)的摩阻和熱流結(jié)果偏高，而NCCR模擬結(jié)果更趨近于DSMC粒子仿真結(jié)果。