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

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

引 言

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

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

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

1 NCCR方程理論基礎

1.1 廣義速度矩

對于流動問題, 通常可以將流場中的宏觀物理量分為守恒量和非守恒量, 宏觀量與微觀量的對應關系可以表示為

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

式中, 角括號代表對氣體分子在整個相速度空間的積分,Φ(k)代表兩組宏觀物理量,h(k)為這些宏觀量對應的微觀表達式。

對于守恒量, 其宏觀與微觀定義如下

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

式中,ρ,U和ρe分別為氣體的宏觀密度、 速度和內能,m,ξ和Hrot分別表示粒子的微觀質量, 速度和轉動Hamilton算子(僅針對雙原子氣體)。

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

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

其中,I代表單位張量。

1.2 非守量的演化方程

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

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

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

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

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

(1)

(2)

(3)

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

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

1.3 高階矩封閉與絕熱假設

(1)高階矩封閉

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

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

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

(2)絕熱假設

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

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

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

其中,

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

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

2 非線性耦合本構關系求解方法

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

2.1 解耦求解方法

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

2.2 耦合求解方法

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

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

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

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

F(x)=0

(4)

(1)不動點迭代法

對于雙原子NCCR本構關系模型, 首先嘗試采用不動點迭代法進行求解和分析。非線性方程組的一般形式(4)可改寫成等價的不動點形式[23]為

x=G(x)

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

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

圖1 不動點迭代法收斂情況[23] Fig. 1 Convergence of fixed-point iteration method[23]

(2)Newton迭代法

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

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

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

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

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

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

(5)

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

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

然后, 在舊時刻的迭代值x(k-1)基礎上, 通過迭代增量Δ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ù)值求解與驗證

3.1 一維激波結構

一維激波結構是驗證氣體數(shù)值模型精度的經典流動算例。其不僅能夠表征高Mach數(shù)流動中的熱力學非平衡特點, 同時能夠將邊值問題轉變?yōu)槌踔祮栴}, 從而避免了邊界條件的處理。因此, 對于NCCR方程的研究首先從一維激波結構[25]的模擬出發(fā)。

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

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

圖5 氮氣激波結構的無量綱密度厚度[25]Fig. 5 Comparison of dimensionless density thickness of Nitrogen shock structure[25]

3.2 二維圓柱繞流

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

圖6 氮氣圓柱繞流速度云圖[23]Fig. 6 Comparison of velocity contours for Nitrogen gas[23]

圖7 駐點線速度分布[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 背風區(qū)Mach數(shù)云圖對比(H=50 km)[26]Fig. 10 Contour lines of Mach number of after-body flows(H=50 km)[26]

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

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

4 改進模型NCCR+

(6)

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

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

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

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

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

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

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

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

5 NCCR方程耦合轉動非平衡效應

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

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

當f(γ)等于γ-1時表示由動理學理論得到的相應系數(shù), 當其等于γ時表示由連續(xù)介質理論得到的相應系數(shù)。上式中Zrot表示轉動碰撞數(shù), 其具體表達式如下

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

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

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

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

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

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

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

圖21 轉動溫度云圖[34]Fig. 21 Rotational temperature contours[34]

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

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

6 NCCR與熱化學非平衡的耦合計算

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

圖 23對比了81 km高度采用NCCR與N-S方程計算所得的溫度云圖。可以看到, 在考慮真實氣體效應情況下, N-S方程和NCCR方程預測的結果差別較大, 尤其是激波脫體距離。

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

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

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

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

7 總結與展望

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