稀薄氣體流動(dòng)非線性耦合本構(gòu)關(guān)系研究進(jìn)展

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

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

引 言

為了更好地研究跨流域多尺度問題, 一般可以根據(jù)Knudsen數(shù)將流域劃分為不同的區(qū)間。Knudsen數(shù)是由丹麥物理學(xué)家Martin Knudsen首次提出的無量化參數(shù), 其定義為分子平均自由程與流動(dòng)特征長度的比值, 以下簡記為Kn。按照Kn的大小, 流動(dòng)可以劃分為4種流域: 連續(xù)流(Kn<0.01), 滑移流(0.0110)。對(duì)于連續(xù)流域的大量流動(dòng)問題, 傳統(tǒng)的Navier-Stokes(N-S)方程已經(jīng)擁有足夠的建模精度。然而, 對(duì)于高超聲速跨流域多尺度流動(dòng)問題, N-S方程的建模精度會(huì)大幅降低。因?yàn)? 其服從Newton黏性定律和Fourier熱傳導(dǎo)定律模型的線性本構(gòu)模型不足以描述非平衡領(lǐng)域的高度非線性問題。為了彌補(bǔ)這一缺陷, 含有高階非線性項(xiàng)的本構(gòu)方程相繼提出。包括基于Chapman-Enskog (CE) 2階與3階展開得到的Burnett方程[1]與Super Burnett方程[2],以及在Burnett方程的基礎(chǔ)上, 進(jìn)行簡化處理得到的一系列Burnett類方程[3-6]。除了基于 CE 展開得到的高階方程組, Grad[7]基于Hermite多項(xiàng)式展開構(gòu)造了Grad 13矩分布函數(shù), 并推導(dǎo)了最初的Grad 13 矩方程。為了克服Grad矩方程雙曲特性的缺陷, 學(xué)者們提出了正則化的矩方程[8-11]。雖然與傳統(tǒng)的N-S方程相比, Burnett方程和矩方程在多尺度流動(dòng)的模擬中均能夠顯現(xiàn)出一定優(yōu)勢(shì)。然而, 這類宏觀方程的數(shù)值穩(wěn)定性以及數(shù)值求解精度依然不能夠滿足實(shí)際的工程需求。

為了提供一種可靠并能夠?qū)崿F(xiàn)穩(wěn)定計(jì)算的高階流體動(dòng)力學(xué)模型, Eu[12-15]從廣義流體動(dòng)力學(xué)的基本理論出發(fā)并結(jié)合非平衡集成方法, 提出了一組廣義流體動(dòng)力學(xué)方程組(generalized hydrodynamic equations, GHE)。在建模過程中, Eu首先通過構(gòu)造一套指數(shù)型非平衡態(tài)分布函數(shù), 巧妙地將非平衡態(tài)到平衡態(tài)的熵增特性與宏觀非守恒量的耗散過程緊密結(jié)合起來。在此基礎(chǔ)上, 將該指數(shù)型分布函數(shù)代入Boltzmann碰撞項(xiàng)進(jìn)行建模。最終, 得到了嚴(yán)格滿足熱力學(xué)熵增定律的廣義流體動(dòng)力學(xué)方程組。然而, 該方程組的高度非線性與強(qiáng)耦合性, 使其僅限于研究一維線性化問題[16]。例如, 聲波在雙原子氣體, 如N2, H2, D2, HD中的吸收與色散問題。為了將GHE拓展至實(shí)際問題的研究中, Eu通過提出絕熱假設(shè)和封閉理論分別對(duì)高階非守恒量的時(shí)間導(dǎo)數(shù)項(xiàng)和高階矩的通量項(xiàng)進(jìn)行簡化處理。最終, 得到簡化的廣義流體動(dòng)力學(xué)方程組, 并成功運(yùn)用到單原子[17]和雙原子氣體[18-19]的一維激波結(jié)構(gòu)模擬中。當(dāng)計(jì)算Mach數(shù)從1.2～10變化時(shí)均能得到光滑激波解。為了構(gòu)建一套更易于求解的純代數(shù)形式的方程組, Myong[20-21]通過忽略簡化的廣義流體動(dòng)力學(xué)方程組中非守恒量的梯度項(xiàng)與熱流和速度梯度的點(diǎn)乘項(xiàng), 建立了非線性耦合本構(gòu)關(guān)系(nonlinear coupled constitutive relations, NCCR)模型。相比于傳統(tǒng)N-S方程的線性本構(gòu)模型, NCCR方程中應(yīng)力與熱流的表達(dá)形式具有高非線性與強(qiáng)耦合性的特征。雖然NCCR為捕捉非平衡流動(dòng)的非線性特征提供了堅(jiān)實(shí)的理論基礎(chǔ), 但是相比于傳統(tǒng)的線性本構(gòu)體系, 其求解難度大大增加。為了有效地求解NCCR方程, Myong提出了一種解耦求解算法。該方法的主要思想是將一個(gè)多維問題分解為各個(gè)方向獨(dú)立的一維問題, 最終以線性疊加的方式得到NCCR方程的收斂解。解耦方法已經(jīng)成功地將NCCR方程組運(yùn)用到高超聲速非平衡流動(dòng)問題中, 如一維激波結(jié)構(gòu)[20]、 平板繞流[22]等。在此研究基礎(chǔ)上, Myong通過增加附加體積應(yīng)力的非線性演化方程, 將NCCR方程進(jìn)一步拓展到雙原子氣體的研究中。有效地模擬了雙原子氣體一維激波結(jié)構(gòu)以及二維鈍頭繞流等問題[21], 并取得了與實(shí)驗(yàn)或DSMC方法較為吻合的結(jié)果。

然而, 解耦思想忽略了真實(shí)流動(dòng)中各方向的耦合特性, 同時(shí)弱化了NCCR方程本身固有的耦合特征。在實(shí)際的研究過程中也發(fā)現(xiàn), 解耦求解方法對(duì)于三維復(fù)雜流動(dòng)問題存在收斂困難以及計(jì)算不穩(wěn)定的缺陷。尤其在尾部膨脹流域, 極易出現(xiàn)負(fù)密度等非物理現(xiàn)象。為了充分地發(fā)揮NCCR方程的耦合特性, 提出了一種耦合迭代算法[23-24]。該方法充分結(jié)合了不動(dòng)點(diǎn)迭代法與Newton迭代法各自的優(yōu)勢(shì), 無論在強(qiáng)壓縮區(qū)或膨脹區(qū)均能得到NCCR方程的收斂解。大量的數(shù)值驗(yàn)證結(jié)果表明, 耦合算法不僅在激波結(jié)構(gòu)[25]、 圓柱繞流[23]、 鈍頭繞流[26]等標(biāo)準(zhǔn)模型上得到了有效的驗(yàn)證, 同時(shí), 也能夠成功地拓展到高超聲速復(fù)雜工程問題的流場結(jié)構(gòu)和氣動(dòng)力熱預(yù)測(cè)中, 如空心擴(kuò)張圓管流動(dòng)、 HTV-2飛行器[26]和Apollo返回艙等。為了進(jìn)一步提高NCCR方程計(jì)算氣動(dòng)力、 熱的精度, 結(jié)合耦合求解算法提出了關(guān)于NCCR方程的高精度邊界條件[27]。關(guān)于NCCR方程的拓展研究還有很多, 例如, 在NCCR的本構(gòu)體系中考慮了非定常項(xiàng)用于研究氣動(dòng)加熱問題[28-29]; 在氣體動(dòng)理學(xué)格式[30-31]的框架下, 求解NCCR方程, 提出的拓展氣體動(dòng)理學(xué)格式[32]等。然而, Myong提出的傳統(tǒng)NCCR方程對(duì)熱通量演化方程的簡化處理, 在一定程度上降低了原始方程的計(jì)算精度。為了使NCCR發(fā)揮出更高的計(jì)算精度, 提出了一種改進(jìn)的NCCR+方程[33], 并同樣采用耦合求解的方法對(duì)NCCR+方程進(jìn)行了有效的求解。通過模擬Apollo返回艙再入過程迎風(fēng)面強(qiáng)激波壓縮流和不同擴(kuò)張拐角高超聲速膨脹流, 表明NCCR+方程的確能夠在強(qiáng)激波壓縮區(qū)域和膨脹區(qū)域提高傳統(tǒng)NCCR方程的計(jì)算精度,具有更強(qiáng)的非平衡流動(dòng)模擬能力。在這些研究的基礎(chǔ)上, 考慮到高超聲速流動(dòng)中不僅存在多尺度效應(yīng), 同時(shí)也存在顯著的熱力學(xué)與熱化學(xué)非平衡效應(yīng), 于是將NCCR方程的非線性效應(yīng)與多物理效應(yīng)充分地結(jié)合, 構(gòu)建了NCCR方程的雙溫度計(jì)算模型[34], 以及NCCR方程的熱化學(xué)反應(yīng)的耦合計(jì)算模型[35]。這些模型已成功地運(yùn)用于高超聲速復(fù)雜外形的計(jì)算, 用于揭示稀薄氣體效應(yīng)與真實(shí)氣體效應(yīng)的耦合作用機(jī)理。本文將圍繞NCCR方程的理論研究與工程應(yīng)用兩方面, 對(duì)近年來NCCR方程的研究進(jìn)展進(jìn)行綜述性回顧。

1 NCCR方程理論基礎(chǔ)

1.1 廣義速度矩

對(duì)于流動(dòng)問題, 通常可以將流場中的宏觀物理量分為守恒量和非守恒量, 宏觀量與微觀量的對(duì)應(yīng)關(guān)系可以表示為

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

式中, 角括號(hào)代表對(duì)氣體分子在整個(gè)相速度空間的積分,Φ(k)代表兩組宏觀物理量,h(k)為這些宏觀量對(duì)應(yīng)的微觀表達(dá)式。

對(duì)于守恒量, 其宏觀與微觀定義如下

Φ(1)=ρ,Φ(2)=ρU,Φ(3)=ρe

式中,ρ,U和ρe分別為氣體的宏觀密度、 速度和內(nèi)能,m,ξ和Hrot分別表示粒子的微觀質(zhì)量, 速度和轉(zhuǎn)動(dòng)Hamilton算子(僅針對(duì)雙原子氣體)。

非守恒量的宏觀與微觀定義如下

Φ(4)=Π,Φ(5)=Δ,Φ(6)=Q

其中,I代表單位張量。

1.2 非守量的演化方程

由于守恒量的演化方程嚴(yán)格滿足質(zhì)量、 動(dòng)量、 能量守恒關(guān)系, 且有明確統(tǒng)一的表達(dá)形式, 不再贅述。而非守恒量的演化方程, 由于建模方式的不同, 最終的表達(dá)形式不盡相同。為了對(duì)非守恒量的演化方程進(jìn)行建模, 同時(shí)將非平衡態(tài)到平衡態(tài)的熵增特性與宏觀非守恒量的耗散過程緊密結(jié)合起來, Eu提出了一種指數(shù)型非平衡態(tài)分布函數(shù)[12]。在無外力作用下, 單組分氣體分子的非平衡態(tài)分布函數(shù)表達(dá)形式為

其中,T,kB分別代表溫度和Boltzmann常數(shù),μ是歸一化因子,X(k)是與宏觀量有關(guān)的系數(shù), 其表達(dá)形式為

X(4)=-Φ(4)/(2p)=-Π/(2p)

X(5)=-3Φ(5)/(2p)=-3Δ/(2p)

Eu基于該非平衡態(tài)分布函數(shù), 并通過碰撞累積量展開的方式對(duì)Boltzmann的碰撞項(xiàng)建模。并在此基礎(chǔ)上, 結(jié)合廣義速度矩的演化方程, 最終得到了非守恒量的演化方程, 如下

(1)

(2)

(3)

其中,ψ(1),ψ(2),ψ(3),ψ(p)代表高階矩;γ′=(5-3γ)/2;q(κ)是非線性耗散項(xiàng), 其具體表達(dá)形式為sinh(κ)/κ; sinh(κ)是雙曲正弦函數(shù), sinh(κ)=(eκ-e-κ)/2,κ通過Rayleigh-Onsager耗散方程給定, 其表達(dá)形式為

從其表達(dá)形式可以看出,κ是一個(gè)非常核心的參數(shù), 該無量綱參數(shù)綜合了與氣體相關(guān)的物理信息: (1)與氣體屬性相關(guān)的信息, 其中包括分子質(zhì)量(m)、 直徑(d)、 剪切黏性系數(shù)(η)、 體積黏性系數(shù)(ηb)和熱傳導(dǎo)系數(shù)(λ); (2)與流場宏觀物理量相關(guān)的信息, 其中包括溫度(T)和壓力(p); (3)與非守恒量相關(guān)的信息, 包括應(yīng)力(Π)、 附加體積應(yīng)力(Δ)和熱流(Q)。

1.3 高階矩封閉與絕熱假設(shè)

(1)高階矩封閉

通過式(1)～(3)不難發(fā)現(xiàn), 雖然得到了高階矩的演化方程, 但是廣義流體動(dòng)力學(xué)方程組依然是一個(gè)有待封閉的系統(tǒng)。因?yàn)楦唠A矩依然沒有明確的表達(dá)形式， 所以為了封閉方程組必須對(duì)高階矩進(jìn)行建模。Eu為了封閉上述方程組, 提出了一套有別于Crad封閉的簡單封閉方法[17], 即假設(shè)高階矩為零

ψ(1)=ψ(2)=ψ(3)=0

這樣的封閉方式, 恰好能與等式右端具有2階碰撞精度的q(κ)保持一致。雖然Eu與Myong對(duì)高階矩的封閉方法有各自不同的看法, 但是基于兩種封閉理論的簡化效果是一致的, 即不考慮高階矩在演化方程中的作用。

(2)絕熱假設(shè)

其實(shí)NCCR方程的表達(dá)形式與N-S方程既有聯(lián)系, 又有區(qū)別。一方面該方程組保留了N-S本構(gòu)方程中的線性部分, 另一方面也有其特有的非線性效應(yīng)。首先, 在NCCR方程中的應(yīng)力Π和附加體積應(yīng)力Δ與速度的各梯度分量, 不再是線性關(guān)系。同樣地, 熱流Q與溫度梯度也不再線性相關(guān)。另外, 非線性耗散項(xiàng)q(κ)的引入不僅增強(qiáng)了方程組的非線性程度, 同時(shí)又將非守恒變量耦合在一起。然而, 在傳統(tǒng)的N-S方程中, 可以很明顯地發(fā)現(xiàn), 應(yīng)力和熱流不是直接相關(guān)的變量。為了便于對(duì)比, 將N-S方程的本構(gòu)模型, 羅列如下

1.4 NCCR方程的無量綱化與歸一化

為了更好地實(shí)現(xiàn)NCCR方程的理論求解, 在實(shí)際的應(yīng)用過程中通常采用如下的無量綱體系對(duì)控制方程與非線性耦合本構(gòu)模型進(jìn)行無量綱化處理[37-39]

其中,

與此同時(shí), 為了將NCCR方程寫成極簡練的形式, 引入如下歸一化方式

最終, 可以將歸一化處理后的NCCR方程寫成

2 非線性耦合本構(gòu)關(guān)系求解方法

NCCR模型盡管在GHE方程基礎(chǔ)上得到很大的簡化, 但是由于高階非守恒量的強(qiáng)非線性與強(qiáng)耦合性特征, 難以對(duì)其直接進(jìn)行數(shù)值求解。NCCR目前的求解方法主要有兩種, 解耦求解方法與耦合求解方法。

2.1 解耦求解方法

Myong提出的NCCR模型的非線性代數(shù)方程組解耦求解方法[21-22], 實(shí)則將一個(gè)耦合系統(tǒng)分裂為單一方向的解耦系統(tǒng)。并在此基礎(chǔ)上對(duì)每個(gè)方向獨(dú)立求解, 最終通過簡單線性迭代的方式給出NCCR方程的收斂解。Myong的解耦算法[40]將三維問題近似分裂為x,y,z這3個(gè)方向的一維問題來處理。在三維有限控制體的一個(gè)面上(如x面), 由熱力學(xué)力(ux,vx,wx,Tx)引起的應(yīng)力和熱流分量(Πxx,Πxy,Πxz,Δ,Qx), 可以近似成兩個(gè)分裂求解器的線性疊加。其中一個(gè)求解器求解壓縮膨脹流動(dòng)(ux,0,0,Tx), 另一個(gè)計(jì)算剪切流動(dòng)(0,vx,0,0)和(0,0,wx,0)。

2.2 耦合求解方法

NCCR模型的解耦求解算法雖然在一維激波結(jié)構(gòu)和二維簡單外形問題上初步展示它的能力。然而, 這種分裂思想忽略真實(shí)流動(dòng)中3個(gè)方向的相互影響, 弱化了本構(gòu)關(guān)系中非守恒量之間固有的耦合關(guān)系。因此, 本節(jié)重點(diǎn)闡述如何通過一種耦合思想完整地求解非線性耦合本構(gòu)關(guān)系模型。首先, 在剪切應(yīng)力無跡和對(duì)稱關(guān)系條件下, 應(yīng)力滿足如下關(guān)系

因此, 適用于雙原子氣體分子的NCCR本構(gòu)方程可表達(dá)成如下非線性方程組的一般形式

式中, 每個(gè)函數(shù)fi可以看作是將一個(gè)九維實(shí)數(shù)空間里的自變量映射到一維實(shí)數(shù)空間的關(guān)系。若將含有9個(gè)未知量的9個(gè)非線性方程表示成一個(gè)矢量函數(shù)形式, 則

那么, 雙原子NCCR模型的一般形式可簡化為

F(x)=0

(4)

(1)不動(dòng)點(diǎn)迭代法

對(duì)于雙原子NCCR本構(gòu)關(guān)系模型, 首先嘗試采用不動(dòng)點(diǎn)迭代法進(jìn)行求解和分析。非線性方程組的一般形式(4)可改寫成等價(jià)的不動(dòng)點(diǎn)形式[23]為

x=G(x)

雖然不動(dòng)點(diǎn)迭代收斂速度只有1階, 甚至越接近解的區(qū)域迭代收斂速度越慢, 但其對(duì)迭代初值依賴性不大。因此, 不需要給出一個(gè)充分接近解的初值。不動(dòng)點(diǎn)迭代最大的難點(diǎn)在于如何構(gòu)造一個(gè)滿足不動(dòng)點(diǎn)迭代收斂定理的多變量顯式迭代式。所以, 在構(gòu)建不動(dòng)點(diǎn)迭代式之前, 首先將無量綱歸一化處理后的雙原子NCCR模型轉(zhuǎn)換成如下形式

考慮到以上方程形式依然復(fù)雜, 可在Rayleigh-Onsager耗散函數(shù)形式的基礎(chǔ)上將上式合并成一個(gè)表達(dá)式

圖1 不動(dòng)點(diǎn)迭代法收斂情況[23] Fig. 1 Convergence of fixed-point iteration method[23]

(2)Newton迭代法

由不動(dòng)點(diǎn)迭代收斂性分析可知, 給雙原子NCCR模型構(gòu)造一個(gè)在所有范圍均能完全收斂的不動(dòng)點(diǎn)迭代式是比較困難的。因此, 繼續(xù)嘗試采用另一種求解算法——Newton迭代法。首先, 將NCCR模型轉(zhuǎn)換成另一種更廣義的形式[23]

G(x)=x-A(x)-1·F(x)

其中,A(x)是從九維實(shí)空間映射到一維實(shí)空間的函數(shù)矩陣, 在G的不動(dòng)點(diǎn)p處非奇異。在非線性方程組的Newton迭代法中,A(x)的一個(gè)合適形式是Jacobi矩陣J(x)。從x(0)初值開始, 迭代過程不斷演化推進(jìn)求解。對(duì)于k≥1, 有如下關(guān)系式

x(k)=G(x(k-1))

是的，李莉現(xiàn)在覺得很幸福，是那種能摸得著地的幸福。梅子說的許峰，是她記憶深處的人，如若不提，她幾乎忘記他了。

=x(k-1)-J(x(k-1))-1F(x(k-1))

(5)

雖然, Newton迭代法一般具有二次收斂速度, 但是過度依賴于一個(gè)充分準(zhǔn)確的初值, 而且Jacobi矩陣必須存在。此外, 式中的迭代需要在每一步計(jì)算Jacobi矩陣的逆矩陣, 這樣會(huì)存在計(jì)算量巨大和計(jì)算復(fù)雜性增加的缺陷。為了避免對(duì)Jacobi矩陣直接求逆, 式(5)的計(jì)算過程通過兩步方式代替實(shí)現(xiàn)。首先采用Gauss消去法求解線性方程組, 以獲得每一步的迭代增量Δx

J(x(k-1))Δx=-F(x(k-1))

然后, 在舊時(shí)刻的迭代值x(k-1)基礎(chǔ)上, 通過迭代增量Δx以獲得新的近似解x(k), 如

x(k)=x(k-1)+Δx

圖2 Newton迭代法收斂情況[23]Fig. 2 Convergence of Newton′s iteration method[23]

(3)混合迭代法

圖3 混合迭代算法收斂情況[23]Fig. 3 Convergence of hybrid iteration method[23]

3 NCCR數(shù)值求解與驗(yàn)證

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

一維激波結(jié)構(gòu)是驗(yàn)證氣體數(shù)值模型精度的經(jīng)典流動(dòng)算例。其不僅能夠表征高M(jìn)ach數(shù)流動(dòng)中的熱力學(xué)非平衡特點(diǎn), 同時(shí)能夠?qū)⑦呏祮栴}轉(zhuǎn)變?yōu)槌踔祮栴}, 從而避免了邊界條件的處理。因此, 對(duì)于NCCR方程的研究首先從一維激波結(jié)構(gòu)[25]的模擬出發(fā)。

圖4展示了通過NCCR方程求解氬氣激波結(jié)構(gòu)的無量綱密度厚度倒數(shù)結(jié)果。從圖中的對(duì)比結(jié)果可以發(fā)現(xiàn), 采用VHS模型且取黏性指數(shù)s=0.81, 能夠使得NCCR的解與實(shí)驗(yàn)結(jié)果吻合最好。然而, 采用Sutherland黏性的NCCR方程預(yù)測(cè)的結(jié)果與實(shí)驗(yàn)值相比有一定差距。圖5詳細(xì)對(duì)比了NCCR方程求解雙原子氣體得到的無量綱密度厚度倒數(shù)與實(shí)驗(yàn)結(jié)果。由密度厚度倒數(shù)曲線可見, 傳統(tǒng)NSF方程在不考慮體積黏性情況下預(yù)測(cè)結(jié)果最差。然而, NCCR方程在大Mach數(shù)范圍內(nèi)與實(shí)驗(yàn)測(cè)量值吻合最好。同時(shí), 也可以發(fā)現(xiàn), 不管是NSF方程還是NCCR方程, 體積黏性均有助于提高對(duì)雙原子氣體流動(dòng)的預(yù)測(cè)精度。

圖4 氬氣激波結(jié)構(gòu)的無量綱密度厚度[25]Fig. 4 Comparison of dimensionless density thickness of Argon shock structure[25]

圖5 氮?dú)饧げńY(jié)構(gòu)的無量綱密度厚度[25]Fig. 5 Comparison of dimensionless density thickness of Nitrogen shock structure[25]

3.2 二維圓柱繞流

為進(jìn)一步驗(yàn)證NCCR方程的模擬精度, 在有限體積計(jì)算框架下實(shí)現(xiàn)了NCCR方程對(duì)簡單二維外形的求解。圖6展示了NCCR方程模擬的高超聲速圓柱繞流無量綱速度云圖。圖7展示了由解耦算法與耦合算法預(yù)測(cè)的NCCR速度曲線與DSMC的對(duì)比結(jié)果。不難發(fā)現(xiàn), NCCR方程采用耦合算法預(yù)測(cè)的結(jié)果與DSMC在激波處十分吻合。然而, NCCR的解耦算法預(yù)測(cè)的結(jié)果與DSMC有較為明顯的偏差。這也意味著解耦算法忽略掉NCCR模型的重要流動(dòng)特征信息的確會(huì)降低方程的模擬精度。因此, 耦合收斂算法對(duì)NCCR模型求解的穩(wěn)定性和精度有著非常重要的影響。

圖6 氮?dú)鈭A柱繞流速度云圖[23]Fig. 6 Comparison of velocity contours for Nitrogen gas[23]

圖7 駐點(diǎn)線速度分布[23]Fig. 7 Comparison of velocity profiles along stagnation line[23]

3.3 類HTV-2飛行器

圖8 Mach數(shù)云圖(H=50 km)[26]Fig. 8 Mach number contours(H=50 km)[26]

圖9 Mach數(shù)云圖(H=90 km)[26]Fig. 9 Mach number contours(H=90 km)[26]

圖10 背風(fēng)區(qū)Mach數(shù)云圖對(duì)比(H=50 km)[26]Fig. 10 Contour lines of Mach number of after-body flows(H=50 km)[26]

圖11 背風(fēng)區(qū)溫度云圖對(duì)比(H=90 km)[26]Fig. 11 Contour lines of Mach number of after-body flows(H=90 km)[26]

因此, 可以看出在連續(xù)流區(qū)域, NCCR方程的模擬結(jié)構(gòu)能夠恢復(fù)到N-S方程的結(jié)果。在N-S方程失效的非平衡區(qū)域, NCCR方程能夠體現(xiàn)出其差異性。

4 改進(jìn)模型NCCR+

(6)

若對(duì)方程(6)進(jìn)行無量綱化與歸一化處理, 可以將其表達(dá)形式寫為

為了驗(yàn)證NCCR+方程的模擬精度, 基于高超聲速擴(kuò)張拐角的內(nèi)流場, 進(jìn)行了詳細(xì)的對(duì)比分析。從圖 12的對(duì)比結(jié)果可以明顯地發(fā)現(xiàn), 在擴(kuò)張區(qū)域采用DSMC預(yù)測(cè)所得的溫度場與基于NCCR+模擬的數(shù)值結(jié)果吻合較好。為了進(jìn)行定量的分析研究, 截取x=0.15 mm溫度曲線, 并將NCCR+的預(yù)測(cè)結(jié)果與N-S, NCCR 和DSMC一同對(duì)比分析。從圖13的對(duì)比結(jié)果中均可以發(fā)現(xiàn), 擴(kuò)張區(qū)域由于強(qiáng)烈的稀薄氣體非平衡效應(yīng), N-S求解器的數(shù)值計(jì)算結(jié)果與DSMC相差甚遠(yuǎn)。而采用NCCR+方程預(yù)測(cè)所得的溫度曲線與DSMC的結(jié)果吻合得很好。而采用傳統(tǒng)NCCR方程預(yù)測(cè)所得的結(jié)果也與DSMC方法的模擬值比較接近。從以上分析結(jié)果中可以看到, NCCR+方程對(duì)于非平衡擴(kuò)張區(qū)域的模擬具有潛在的優(yōu)勢(shì)。同時(shí)也可以發(fā)現(xiàn)在傳統(tǒng)NCCR方法計(jì)算過程中, 對(duì)熱流本構(gòu)關(guān)系的簡化處理確實(shí)在一定程度會(huì)降低模型的計(jì)算精度。

圖12 DSMC和NCCR+溫度云圖對(duì)比[33]Fig. 12 Comparison of temperature contours between DSMC and NCCR+[33]

圖13 截面處溫度曲線對(duì)比[33]Fig. 13 Comparison of temperature profiles[33]

為了驗(yàn)證采用NCCR+方法, 模擬三維工程外形的能力, 分別對(duì)95 km和105 km的Apollo返回艙再入過程進(jìn)行了數(shù)值模擬與對(duì)比。從圖 14和15中的壓力對(duì)比云圖可以明顯發(fā)現(xiàn), 采用NCCR+預(yù)測(cè)得到的激波脫體距離明顯大于N-S方程, 更能體現(xiàn)出稀薄氣體非平衡效應(yīng)。

圖14 壓力云圖對(duì)比(H=95 km)[33]Fig. 14 Comparison of pressure contours (H=95 km)[33]

從圖 16和17物面摩阻系數(shù)對(duì)比曲線可以發(fā)現(xiàn), 在95 km與105 km高度, 基于NCCR+方程模擬得到的摩阻系數(shù)在迎風(fēng)區(qū)與背風(fēng)區(qū)均與DSMC結(jié)果吻合較好。而基于傳統(tǒng)N-S方程模擬得到的摩阻系數(shù)在背風(fēng)區(qū)結(jié)果明顯高于DSMC的預(yù)測(cè)值。同時(shí)在迎風(fēng)區(qū), NCCR+預(yù)測(cè)所得的熱流曲線比NCCR更接近于DSMC的預(yù)測(cè)結(jié)果。

圖17 物面摩阻系數(shù)對(duì)比(H=105 km)[33]Fig. 17 Comparison of surface friction coefficient(H=105 km)[33]

5 NCCR方程耦合轉(zhuǎn)動(dòng)非平衡效應(yīng)

針對(duì)高超聲速氣體而言, 氣體的壓縮膨脹效應(yīng)十分顯著, 體積黏性表征了氣體的壓縮膨脹效應(yīng)。因此, 對(duì)于體積黏性的研究也尤為必要。通常為了更好地描述體積黏性, 引入黏性系數(shù)比, 即體積黏性系數(shù)與剪切黏性系數(shù)的比值, 其表達(dá)形式如下

針對(duì)不同種類的氣體, 由于氣體屬性不同,fb的取值其實(shí)也不相同。Prangsma等[42]利用聲波吸收實(shí)驗(yàn)在77 K下測(cè)得氖氣的ηb=0, 從實(shí)驗(yàn)層面證明了單原子氣體體積黏性等于零的理論預(yù)測(cè)。同時(shí), 在溫度為293 K的測(cè)試環(huán)境下, 分別測(cè)量出了N2, CO, CH4和 CD4, 它們的fb值分別為 0.73, 0.55, 1.33 和1.17。這也表明在常溫下, 這些氣體的體積黏性系數(shù)與剪切黏性系數(shù)基本在同一量級(jí)。在之前NCCR方程的研究中, 針對(duì)雙原子N2而言,fb的值統(tǒng)一近似取為0.8。為了給NCCR方程提供雙原子氮?dú)鈌b更加精確的取值方式, 通過連續(xù)性理論和分子動(dòng)理學(xué)理論給出了fb的取值的精確表達(dá)形式, 具體如下

當(dāng)f(γ)等于γ-1時(shí)表示由動(dòng)理學(xué)理論得到的相應(yīng)系數(shù), 當(dāng)其等于γ時(shí)表示由連續(xù)介質(zhì)理論得到的相應(yīng)系數(shù)。上式中Zrot表示轉(zhuǎn)動(dòng)碰撞數(shù), 其具體表達(dá)式如下

如圖18所示, 在300 K左右, 相比于連續(xù)介質(zhì)理論下的體積黏性系數(shù), 基于動(dòng)理學(xué)理論得到的黏性系數(shù)比fb與實(shí)驗(yàn)結(jié)果更加吻合。同時(shí), 在圖中也可以看到, 氮?dú)夥肿芋w積黏性系數(shù)與剪切黏性系數(shù)的比值隨著溫度的升高而逐漸增大。在600 K左右時(shí), 近似等于0.8, 這也給出了傳統(tǒng)NCCR理論在高超聲速建模過程中黏性系數(shù)比取值為0.8的合理性區(qū)間。

圖18 黏性系數(shù)比[34]Fig. 18 Viscosity ratio[34]

雙原子氣體除了體積黏性效應(yīng), 在高M(jìn)ach數(shù)流動(dòng)問題中, 由于轉(zhuǎn)動(dòng)溫度向平動(dòng)溫度松弛過程的相對(duì)弛豫時(shí)間往往較長, 因此轉(zhuǎn)動(dòng)非平衡的效應(yīng)尤為明顯。為了更好地描述高超稀薄環(huán)境下的轉(zhuǎn)動(dòng)非平衡流動(dòng)過程, 在NCCR方程的基礎(chǔ)上加入了轉(zhuǎn)動(dòng)能松弛源項(xiàng),最終實(shí)現(xiàn)了可以考慮雙溫度非平衡效應(yīng)的NCCR模型[34]。

圖19比較了帶轉(zhuǎn)動(dòng)能松弛源項(xiàng)的N-S方程和NCCR方程計(jì)算所得的Mach數(shù)云圖。從圖中可以看到, 在滑移流區(qū), 即Kn=0.05, 兩者略有差異; 當(dāng)Kn=0.25到達(dá)過渡流區(qū), 流場非平衡程度顯著, 采用NCCR方法預(yù)測(cè)的激波脫體距離明顯大于N-S方程預(yù)測(cè)的結(jié)果, 其對(duì)非平衡效應(yīng)的捕捉能力更為突出。

(a) Kn=0.05 (b) Kn=0.25圖19 Mach數(shù)對(duì)比云圖[34]Fig. 19 Comparison of Mach number contours[34]

如圖 20所示, 在Kn=0.05的滑移流區(qū), NCCR與N-S方程預(yù)測(cè)駐點(diǎn)線溫度曲線差異明顯。相比于N-S方程預(yù)測(cè)的結(jié)果, NCCR方程預(yù)測(cè)的激波層耗散性更強(qiáng), 并且在平動(dòng)和轉(zhuǎn)動(dòng)溫度分布上均與DSMC模型更吻合。

激波與激波邊界層相互作用[41]是高超聲速航天器繞流問題中一個(gè)常見的現(xiàn)象。高超聲速尖壁前緣附近的激波與邊界層相互作用, 會(huì)對(duì)平動(dòng)能與轉(zhuǎn)動(dòng)能的交換產(chǎn)生強(qiáng)烈的非平衡效應(yīng)。本算例采用NCCR的雙溫度模型模擬了高超聲速平板繞流, 并與文獻(xiàn)[43]中Tsuboi等的實(shí)驗(yàn)測(cè)量結(jié)果和DSMC的數(shù)值模擬結(jié)果進(jìn)行了比較。如圖 21所示, 在研究過程中采用含雙溫模型的NCCR方程, 給出了平板周圍的轉(zhuǎn)動(dòng)溫度Trot。

圖21 轉(zhuǎn)動(dòng)溫度云圖[34]Fig. 21 Rotational temperature contours[34]

圖22展示了采用雙溫度NCCR計(jì)算的平動(dòng)和轉(zhuǎn)動(dòng)溫度與DSMC和實(shí)驗(yàn)數(shù)據(jù)[39]對(duì)比結(jié)果。從對(duì)比圖中可以看到, 雙溫度NCCR方程的計(jì)算結(jié)果與實(shí)驗(yàn)測(cè)量值吻合很好, 甚至在某些區(qū)域比DSMC方法更接近于實(shí)驗(yàn)測(cè)量結(jié)果。

圖22 截面處溫度曲線對(duì)比[34]Fig. 22 Comparison of temperature profiles[34]

6 NCCR與熱化學(xué)非平衡的耦合計(jì)算

目前, 直接模擬Monte Carlo方法(direct simula-tion Monte Carlo, DSMC)[44]是唯一能夠用于模擬過渡流化學(xué)反應(yīng)的有效方法。但現(xiàn)有DSMC方法對(duì)于近連續(xù)熱化學(xué)非平衡流條件下的流動(dòng)求解仍不能同時(shí)具備高效、 工程可用的特點(diǎn)[45]。因此, 本節(jié)主要論述NCCR方程與化學(xué)反應(yīng)動(dòng)力學(xué)的耦合計(jì)算模型[35,46]。其中, 化學(xué)反應(yīng)模型主要采用Gupta 7組元化學(xué)反應(yīng)動(dòng)力學(xué)模型[47], 包括5種中性空氣組分N2, O2, N, O和NO以及兩種帶電粒子e-和NO+。

圖 23對(duì)比了81 km高度采用NCCR與N-S方程計(jì)算所得的溫度云圖?？梢钥吹? 在考慮真實(shí)氣體效應(yīng)情況下, N-S方程和NCCR方程預(yù)測(cè)的結(jié)果差別較大, 尤其是激波脫體距離。

圖23 N-S與NCCR溫度云圖對(duì)比(81 km)[35]Fig. 23 Comparison of temperature contours between N-S and NCCR (81 km)[35]

除此之外, 將基于理想氣體的NCCR方程與基于真實(shí)氣體效應(yīng)的NCCR方程的計(jì)算結(jié)果進(jìn)行了對(duì)比分析。圖 24和25分別比較了71 km和 81 km 高度下對(duì)稱面溫度分布云圖。從對(duì)比云圖中可以發(fā)現(xiàn), 在考慮真實(shí)氣體效應(yīng)的情況下, 預(yù)測(cè)的激波脫體距離明顯小于理想氣體模型預(yù)測(cè)的結(jié)果。因?yàn)樵谡鎸?shí)氣體效應(yīng)的影響下, 由于化學(xué)產(chǎn)物的生成, 鈍錐前緣的等效密度升高。因此, 稀薄非平衡效應(yīng)減弱, 激波脫體距離變小。這一現(xiàn)象是符合物理規(guī)律的。同時(shí), 可以明顯地發(fā)現(xiàn), 在考慮真實(shí)氣體效應(yīng)的情況下預(yù)測(cè)的流場溫度遠(yuǎn)遠(yuǎn)低于理想氣體效應(yīng)下的預(yù)測(cè)結(jié)果, 這一現(xiàn)象也是符合物理規(guī)律的。

圖24 理想氣體與真實(shí)氣體效應(yīng)下溫分布云圖對(duì)比(H=71 km)[35]Fig. 24 Comparison of temperature contours between perfect gas effect and real gas effect(H=71 km)[35]

圖25 理想氣體與真實(shí)氣體效應(yīng)下溫分布云圖對(duì)比(H=81 km)[35]Fig. 25 Comparison of temperature contours between perfect gas effect and real gas effect (H=81 km)[35]

7 總結(jié)與展望

本文主要從理論研究與工程應(yīng)用兩個(gè)方面綜述回顧了NCCR方程的發(fā)展歷程。在理論研究方面, 主要分析了耦合無分裂算法求解NCCR方程的有效性與必要性。該耦合方法結(jié)合了不動(dòng)點(diǎn)迭代和Newton迭代的各自優(yōu)勢(shì), 彌補(bǔ)了Myong分裂思路對(duì)NCCR方程流動(dòng)的缺陷?；谠摲椒? 成功地將NCCR方程的求解從一維拓展至三維流動(dòng), 為推向工程應(yīng)用提供了堅(jiān)實(shí)的基礎(chǔ)。為了進(jìn)一步提高NCCR方程的求解精度, 提出了改進(jìn)的NCCR+方法, 并通過高超聲速擴(kuò)張拐角流動(dòng)與Apollo返回艙的再入流場, 驗(yàn)證NCCR+的確具有更有效的稀薄非平衡效應(yīng)捕捉能力。為了更加符合高超聲速流動(dòng)熱力學(xué)與熱化學(xué)非平衡的特性, 提出了基于非線性耦合本構(gòu)關(guān)系的多溫度模型, 以及考慮化學(xué)反應(yīng)源項(xiàng)的NCCR方程, 并用于探究稀薄氣體效應(yīng)與真實(shí)氣體效應(yīng)耦合作用下的流動(dòng)機(jī)理。雖然NCCR方程的研究取得了一定進(jìn)展, 但仍有很大的發(fā)展空間, 將對(duì)其進(jìn)行不斷深入的研究, 進(jìn)一步拓展NCCR方程的應(yīng)用范圍。