朱海濤，李巖，蘭子奇

中國航空研究院數(shù)字仿真中心，北京 100012

準確預測高超聲速飛行器周圍的氣動熱力學環(huán)境是設計、發(fā)展和評估新型高超聲速飛行器的關(guān)鍵[1]。由于飛行速度很大，在黏性和強激波壓縮作用下，氣體巨大的動能轉(zhuǎn)化為內(nèi)能，致使高超聲速飛行器所處的氣動熱力學環(huán)境異常苛刻。在地面復現(xiàn)這種高溫高焓環(huán)境非常困難，進而很難通過風洞試驗較為準確地預測高超聲速飛行的氣動熱力學參數(shù)，而且成本高昂。高溫氣體物理化學特性的深入研究和高性能計算技術(shù)的快速發(fā)展使得先進的計算流體力學方法成為預測高超聲速飛行器氣動熱力學特性不可或缺的重要手段[2]。

數(shù)值模擬技術(shù)很早便被應用于高超聲速流場氣動熱力學計算。早在20世紀80年代，C.Park等[3]對高溫高焓空氣的化學反應速率、熱力學參數(shù)、輸運性質(zhì)等物理、化學特性進行了研究和總結(jié)，建立了考慮真實氣體效應的三維熱化學非平衡Navier-Stokes 方程組。此后，高超聲速繞流場數(shù)值模擬技術(shù)得以快速發(fā)展；其中具有代表性的是:M. P.Netterfield 等[4]通過求解二溫度熱化學非平衡Navier-Stokes 方程組獲得了球模型的熱化學非平衡數(shù)值流場，從中提取了激波位置、駐點溫度和組分分布等信息，該算例成為計算高超聲速熱化學非平衡數(shù)值流場的常用對比算例；J.Muylaert 等[5]研究了完全催化邊界條件和完全非催化邊界條件對高超聲速飛行器氣動熱計算結(jié)果的影響，并與風洞試驗和飛行試驗結(jié)果進行對比分析，該算例成為氣動熱計算的典型對比算例。除此之外，高超聲速流動數(shù)值模擬的廣度不斷擴大，開始探索氣動熱力學環(huán)境多物理場耦合數(shù)值模擬。例如，B.Hassan等[6]在氣/固界面處構(gòu)造質(zhì)量、能量平衡的耦合交互面，完成了IRV-2 頭部球錐段熱化學非平衡流動、非定常材料熱響應和燒蝕過程的耦合計算研究，給出了多個彈道點上的表面熱流和溫度分布。隨著火星著落器、星際探測器的興起，飛行器的速度進一步提高，高溫氣體內(nèi)各種非線性物理機制對高超聲速飛行器氣動熱力學環(huán)境的影響已無法忽略。描述這些復雜物理機制的新的數(shù)學模型也不斷被應用于高超聲速流場數(shù)值計算。例如，Marschall等[7]總結(jié)了歐盟航空航天領域使用的多種有限速率催化和化學反應模型，并提出了有限速率表面化學反應的一般框架，用于描述氣/固界面熱量、能量輸運過程。

熱化學非平衡Navier-Stokes 方程組數(shù)值求解技術(shù)也是國內(nèi)的研究熱點。董維中等[8]采用多個振動溫度氣體模型完成了半球模型和球錐模型繞流的數(shù)值模擬，細致研究了氣體模型對激波脫落距離、壁面熱流、溫度分布的影響；高鐵鎖等[9]通過數(shù)值求解熱化學非平衡Navier-Stokes方程組獲得高溫空氣組分和溫度分布，然后基于輻射傳輸方程，利用高溫氣體組分的主要輻射機制，計算分析了輻射加熱對返回艙熱環(huán)境的影響；張子健等[10]提出了一種理論分析和數(shù)值模擬相結(jié)合的兩步漸進法，研究了振動激發(fā)過程對二維斜劈氣動力/氣動熱特性的影響規(guī)律。董維中等[11]開發(fā)了三維熱化學非平衡Navier-Stokes 方程組求解軟件AEROPH Flow，可以考慮熱力學非平衡效應、表面催化效應、燒蝕等因素的影響。鄒學峰等[12]對高超聲速飛行器熱/力/振動/噪聲多學科耦合技術(shù)的研究進展進行了綜述。

高超聲速技術(shù)的蓬勃發(fā)展，促使新型熱防護技術(shù)、主動流動控制技術(shù)以及新型控制舵面等新技術(shù)快速應用于工程實踐。因此，高超聲速飛行器的幾何外形越來越復雜，結(jié)構(gòu)形式越來越多樣化，基于結(jié)構(gòu)網(wǎng)格的數(shù)值模擬技術(shù)已不能很好地滿足工程實際需求。為此，美國國家航空航天局（NASA）啟動了Fundamental Aeronautics Program 項目[13]，發(fā)展非結(jié)構(gòu)網(wǎng)格計算技術(shù)和基于非結(jié)構(gòu)網(wǎng)格的網(wǎng)格自適應技術(shù)即為其中的重要研究內(nèi)容，并用于發(fā)展高可靠性、高效的高超聲速熱化學非平衡流動數(shù)值模擬工具。這個研究領域也引起國內(nèi)研究者的重視。李鵬等[14]在我國自主開源CFD 軟件“風雷”中也開發(fā)了非結(jié)構(gòu)網(wǎng)格計算功能，但沒有看到詳細的計算結(jié)果。

本文采用非結(jié)構(gòu)混合網(wǎng)格技術(shù)求解三維熱化學非平衡Navier-Stokes方程組，與結(jié)構(gòu)網(wǎng)格計算結(jié)果進行對比分析，驗證了數(shù)值求解工具的有效性；同時，針對邊界層網(wǎng)格分布和小組分氣體計算問題進行了計算研究，并給出了實施建議。

1 熱化學非平衡Navier-Stokes方程組

熱化學非平衡Navier-Stokes 方程組由經(jīng)典Navier-Stokes 方程組和描述物理、化學非平衡效應及其耦合效應的模型方程組構(gòu)成。本文數(shù)值計算程序可以計及振動能激發(fā)態(tài)和離解及置換化學反應，尚無法考慮電子能激發(fā)態(tài)、光輻射和電離反應。因此，熱力學非平衡效應采用二溫度模型；化學非平衡效應采用有限速度化學反應模型；二者的耦合效應僅考慮振動溫度與離解化學反應之間的相互影響。本文數(shù)值求解的熱化學非平衡Navier-Stokes方程組為

式中

其中，前ns個方程為各個組分的質(zhì)量守恒方程，其后為三個方向的動量守恒方程，然后是能量守恒方程，最后是振動能方程。U為包含ns個組分的守恒變量:ρi為組分i的密度，ρu,ρE,ρev分別為單位體積混合氣體總動量矢量、總能量和總振動能，ρ為混合氣體總密度，u,E,ev分別為混合氣體質(zhì)量平均速度矢量和單位質(zhì)量的總能量、振動能。Fc為對流通量項，I為單位矩陣，p為混合氣體壓強。Fv為黏性通量項，Ji是組分i的質(zhì)量擴散通量，遵循Fick擴散定律；τ是黏性應力張量；q是熱傳導通量，遵循Fourier 熱傳導定律，由平動—轉(zhuǎn)動溫度和振動溫度兩部分構(gòu)成；qv是僅由振動溫度確定的熱傳導通量項。Q是源項:ωi是化學反應過程中組分i的生成速率，ωvt是平動—轉(zhuǎn)動態(tài)和振動態(tài)能量轉(zhuǎn)換過程中振動能生成項，可由Landau-Teller 理論建立不同能量態(tài)轉(zhuǎn)換松弛方程是單位質(zhì)量組分i的振動能量，ωi是由化學反應導致組分i質(zhì)量變化，進而引起的振動能改變量。這些方程組是不封閉的，需補充以下數(shù)理方程。

1.1 熱力學方程

熱力學方程包括氣體狀態(tài)方程和各能量態(tài)內(nèi)能計算公式。本文采用二溫度模型，即平動—轉(zhuǎn)動溫度Ttr和振動溫度Tv，混合氣體的總能量

式中，ei是組分i平動和轉(zhuǎn)動能的和。對分子組分

原子組分沒有轉(zhuǎn)動和振動能，因此

各組分的分壓力僅與平動溫度有關(guān)，進而混合氣體的狀態(tài)方程為

式中，R為普適氣體常數(shù)；Mi、分別為組分i的摩爾質(zhì)量和振動特征溫度。

平動—轉(zhuǎn)動模式與振動模式之間存在能量交換。由Landao-Teller理論可知振動能的時間變化率為

1.2 化學反應模型

對于有nr個化學反應的模型，組分i的生成率為

Gupta 給出了高溫空氣化學反應平衡系數(shù)Kr,eq的經(jīng)驗擬合公式，進而可得逆反應的化學反應速率系數(shù)為

本文采用5 組分、17 個化學反應的空氣化學模型，式（10）、式（11）中用到的常數(shù)和經(jīng)驗擬合公式可見參考文獻[1]。

1.3 熱力學—化學反應非平衡耦合模型

利用化學反應控制溫度T建立化學反應非平衡與熱力學非平衡之間的耦合關(guān)系。5組分空氣模型包含置換和離解兩類化學反應。采用Park 優(yōu)先離解模型，置換化學反應的控制溫度與平動—轉(zhuǎn)動溫度一致，即T=Ttr；離解反應的控制溫度由平動—轉(zhuǎn)動溫度和振動溫度確定，計算公式為T=。

1.4 輸運系數(shù)

方程組（1）中需要的輸運系數(shù)有:氣體黏度μi、平動—轉(zhuǎn)動導熱系數(shù)、振動導熱系數(shù)λvi和質(zhì)量擴散系數(shù)Di。單個組分的黏度和導熱系數(shù)用Blotter高溫氣體黏度擬合公式和Eucken經(jīng)驗關(guān)系式計算；混合氣體黏度和導熱系數(shù)通過Wilke半經(jīng)驗公式，根據(jù)各個組分進行加權(quán)平均。擴散系數(shù)采用雙組元模型計算。具體計算公式可參見參考文獻[1]。

式（1）、式（5）、式（8）～式（10）和輸運系數(shù)計算公式構(gòu)成封閉的熱化學非平衡Navier-Stokes 方程組。引入適當邊界條件后構(gòu)成適定的偏微分方程組初邊值問題，即可進行數(shù)值求解。本文采用的邊界條件包括超聲速入口和出口邊界條件、對稱邊界條件、等溫壁面邊界條件、完全非催化和完全催化邊界條件。

2 數(shù)值求解方法

對方程組（1）進行空間積分，并采用格點有限元方法離散，在單元控制體上可得如下半離散方程組

式中，Ωi為網(wǎng)格點i對應的控制體積；N(i)是與點i相鄰的網(wǎng)格點集；j是與點i相鄰的點，ΔSij是網(wǎng)格邊ij對應的對偶網(wǎng)格面積。對流通量項采用AUSM格式計算；黏性通量項首先在各個控制體內(nèi)采用最小二乘法計算每個網(wǎng)格點上的空間導數(shù)值，然后對相鄰兩點進行平均獲得對偶網(wǎng)格面上流動變量、輸運系數(shù)和相關(guān)空間導數(shù)值，再計算黏性通量項；源項采用分片常數(shù)近似直接在控制內(nèi)積分。

非平衡效應和各個源項使得離散后線性方程組矩陣條件數(shù)很大，數(shù)值剛性問題嚴重。因此在時間離散時，對流通量項和源項進行隱式處理

完成空間、時間離散后，采用GMRES算法求解所得線性方程組以更新下一時間步流場信息。

2.1 非結(jié)構(gòu)網(wǎng)格布置

相對于傳統(tǒng)亞、跨、超聲速流場，高超聲速流場具有多組分、高熱焓、強激波、強熱化學非平衡的特征?；谇罢吡鲌鎏卣鞯姆墙Y(jié)構(gòu)混合網(wǎng)格布置規(guī)則，能否適用于后者，需進行專門研究。本文以公認結(jié)構(gòu)網(wǎng)格計算結(jié)果為基準，對非結(jié)構(gòu)混合網(wǎng)格y+、邊界層網(wǎng)格規(guī)模、增長率等因素進行細致研究，為高超聲速非結(jié)構(gòu)混合網(wǎng)格技術(shù)的工程應用奠定基礎。

2.2 小組分氣體計算

由于高超聲速流場熱化學非平衡效應，式（13）、式（14）具有很強的數(shù)值剛性。不斷積累的數(shù)值殘差對小組分氣體計算產(chǎn)生嚴重不利影響。本文在計算化學反應生成項時，引入元素摩爾數(shù)守恒關(guān)系式，對一部分組分氣體采用式（10）直接計算，其余組分氣體由元素摩爾數(shù)守恒關(guān)系式計算。對5組分空氣模型而言，根據(jù)元素摩爾數(shù)守恒，有

計算過程中，對分子組分N2，O2，NO 采用式（10）直接計算，而原子組分采用式（16）計算，保證各元素摩爾數(shù)守恒。與公認的結(jié)構(gòu)網(wǎng)格計算結(jié)果對比表明，這種計算方法顯著改善了小組分氣體的計算準確度。

3 數(shù)值結(jié)果與分析

本文采用基于非結(jié)構(gòu)混合網(wǎng)格的三維熱化學非平衡Navier-Stokes 方程組求解程序計算了半球模型、RAMC II模型高超聲速繞流場，并與公認的、采用結(jié)構(gòu)網(wǎng)格的計算結(jié)果進行對比分析，還討論了非結(jié)構(gòu)混合網(wǎng)格邊界層網(wǎng)格點布置和小組分氣體準確計算問題。

3.1 半球模型

半球模型的半徑r= 6.35mm；自由來流速度V∞=5280m/s，靜溫T∞= 293K，靜壓p∞= 664Pa；基于球半徑的雷諾數(shù)Re=14600。壁面采用完全非催化等溫邊界條件，壁溫Tw= 2000K。實際計算中僅取1/4 球面，并采用對稱邊界條件。

采用4套非結(jié)構(gòu)混合網(wǎng)格，著重考察邊界層網(wǎng)格y+、網(wǎng)格增長率γ對數(shù)值結(jié)果的影響。邊界層內(nèi)棱柱層數(shù)目L、邊網(wǎng)層網(wǎng)格點數(shù)nb和整體計算域網(wǎng)格點數(shù)ng等網(wǎng)格參數(shù)見表1；其中y+值是駐點附近區(qū)域的數(shù)值。

表1 4套非結(jié)構(gòu)混合網(wǎng)格參數(shù)Table 1 Main parameters of four unstructured mesh used in present computation

圖1給出了物面熱流計算結(jié)果與公認結(jié)構(gòu)網(wǎng)格計算數(shù)據(jù)[8]的對比圖。由于本文采用全三維計算，圖中熱流值是每隔15°抽取一條經(jīng)線，共計13 條經(jīng)線上熱流的平均值。橫坐標是經(jīng)線弧長與球半徑的比值?？梢钥闯?，第四套網(wǎng)格所得物面熱流與公認解基本吻合。由于公認解采用二維軸對稱熱化學非平衡Navier-Stokes 方程組，無法計及繞緯線方向的流動的影響，因此二者在駐點下游區(qū)域存在相對較大的差異。

圖1 物面熱流(w/m2)曲線與公認解[8]對比圖Fig.1 Sphere model heat flux(w/m2)of present comparison with reference[8]

圖2給出了駐點線上平動—轉(zhuǎn)動溫度和振動溫度分布與公認結(jié)構(gòu)網(wǎng)格計算數(shù)據(jù)的對比圖。在圖2（a）中，第三、第四套網(wǎng)格計算結(jié)果與公認結(jié)果吻合較好，激波位置也吻合良好。在圖2（b）中，非結(jié)構(gòu)網(wǎng)格計算結(jié)果與結(jié)構(gòu)網(wǎng)格有較明顯的差異:本文最高振動溫度出現(xiàn)在激波位置稍下游，與激波比較接近，而結(jié)構(gòu)網(wǎng)格所得最高振動溫度與激波有較遠的距離。造成差異的原因可能有:（1）振動能控制方程與化學反應、內(nèi)能模式轉(zhuǎn)化緊密相關(guān)，本文采用5組分空氣模型，而參考文獻[8]采用了7組分或11組分化學模型；（2）本文是全三維計算，而參考文獻[8]采用二維軸對稱計算，造成流場分布存在差異。

圖2 駐點線上平動—轉(zhuǎn)動溫度和振動溫度分布曲線與公認解[8]對比圖Fig.2 Translational-rotational temperature and vibrational temperature in stagnation line of present computation comparison with reference[8]

以參考文獻[8]為參考解，可以看出熱流對網(wǎng)格分布的要求比流場變量更為苛刻；為了獲得準確的熱流分布結(jié)果，邊界層網(wǎng)格y+值應不超過1，邊界層增長率應介于1.1～1.2，盡可能使邊界層網(wǎng)格點數(shù)不低于總網(wǎng)格點數(shù)的60%。

3.2 RAMC II模型

RAMC Ⅱ是一個球錐模型，頭部球半徑為0.1524m，球錐半頂角為9°，軸向長度為1.295m。計算工況為:飛行高度61km；自由來流靜壓19.2Pa，靜溫244.3K；飛行馬赫數(shù)23.9，基于球半徑雷諾數(shù)19866。物面采用完全非催化等溫邊界條件，壁溫Tw= 1500K。參考文獻[15]給出了采用結(jié)構(gòu)網(wǎng)格的詳細計算結(jié)果。通過與該公認解的對比分析，討論采用非結(jié)構(gòu)混合網(wǎng)格計算時，兩種化學反應源項計算方法對小組分氣體數(shù)值結(jié)果的影響:第一種方法是直接采用式（10）計算各個組分的化學反應生成項；第二種方法是由式（10）和元素摩爾守恒關(guān)系式（16）計算各個組分化學反應生成項。

圖3給出了駐點線上各組分質(zhì)量分數(shù)分布曲線與公認解的對比結(jié)果。其中圖3（a）、圖3（b）分別給出了直接計算法（方法一）和采用組分摩爾數(shù)守恒關(guān)系計算法（方法二）所得結(jié)果與公認解的對比結(jié)果；圖3（c）是兩種計算方法結(jié)果的對比圖?？梢钥闯?，兩種計算方法所得組分質(zhì)量分數(shù)結(jié)果與公認解均有一定差異。參考文獻[15]采用7 組分空氣模型進行二維軸對稱計算，因此與本文結(jié)果有較明顯差異。相對于公認解，采用第二種方法差異更小。由圖3（c）可見，兩種計算方法對NO組分的影響最為顯著。原因是在邊界層內(nèi)混合氣體的溫度已超過8000K，氧分子已完全離解為氧原子；此時，相對于氮氣、氮原子和氧原子，NO組分為小組分氣體，對積累的數(shù)值誤差最為敏感。因此，為保證對小組分氣體的準確計算，應采用元素摩爾數(shù)守恒關(guān)系式計算方法。

圖3 駐點線上各組分質(zhì)量分數(shù)分布曲線與公認解[15]對比圖Fig.3 Mass fraction in stagnation line of present computation comparison with reference[15]

4 結(jié)論

采用非結(jié)構(gòu)混合網(wǎng)格技術(shù)求解三維熱化學非平衡Navier-Stokes 方程組，通過與典型算例公認結(jié)構(gòu)網(wǎng)格數(shù)值結(jié)果進行對比驗證，檢驗了本文數(shù)值方法的正確性和有效性。應用本文計算程序，討論了在采用非結(jié)構(gòu)混合網(wǎng)格進行高超聲速流場數(shù)值模擬時，邊界層網(wǎng)格布置和小組分氣體準確計算問題，得到如下結(jié)論:

（1）盡管非結(jié)構(gòu)混合網(wǎng)格具有顯著的不規(guī)則性，且控制方程組離散后系統(tǒng)數(shù)值剛性很大，但采用該技術(shù)也可以獲得與結(jié)構(gòu)網(wǎng)格較一致的流場和熱流計算結(jié)果。

（2）邊界層網(wǎng)格布置對物面熱流計算結(jié)果影響顯著；對二階格式而言，為保證對物面熱流的正確計算，邊界層網(wǎng)格y+應小于1（0.5左右），網(wǎng)格增長率應不超過1.4（1.1～1.2）。

（3）小組分氣體對數(shù)值誤差積累較為敏感，為準確計算小組分氣體，在化學反應源項計算過程中，應采用元素摩爾數(shù)守恒關(guān)系式。