任 翔，袁軍婭，蔡國飆

（北京航空航天大學，北京，100191）

0 引 言

高速流動廣泛存在于戰(zhàn)略戰(zhàn)術打擊、導彈系統(tǒng)防御、航天器發(fā)射與再入、地外天體探測等方面。高速流動一般指流速Ma＞5的流動。流動中往往存在薄激波層、壁面附近強的熵梯度和粘性耗散，并產生強烈的高溫現象。高溫條件下氣體的振動能、電子能被激發(fā)，且當分子的振動能級和電子能級達到一定程度時發(fā)生分解、電離等一系列化學反應，這稱為高溫真實氣體效應。同時高速飛行往往出現在上層大氣或者低密度環(huán)境中，流動往往伴隨著稀薄效應。此時稀薄效應和高溫真實氣體效應共同造成強烈的熱化學非平衡，這為高速流動和飛行器的氣動特性的準確預測提出挑戰(zhàn)。

考慮到高速流動的跨流域特性，CFD和DSMC方法廣泛用于高速流動的機理研究和工程應用，尤其針對熱化學非平衡過程，CFD和DSMC方法分別從宏觀和微觀角度發(fā)展了多種模型。CFD方法自20世紀60年代提出Lighthill離解模型、Landau-Teller振動松弛模型，到20世紀80年代的Park雙溫模型。而DSMC方法，主要使用L-B（Larsen-Borgnakke）或者量子L-B模型來模擬分子轉動、振動態(tài)的激發(fā)與松弛，先后發(fā)展了總碰撞能（Total Collision Energy，TCE）、振動態(tài)激發(fā)解離（Vibrationally Favored Dissociation，VFD）、廣義碰撞能（Generalized Collision Energy，GCE）和Q-K（Quantum Kinetic）等模型來模擬化學反應過程，并逐步考慮振動-離解耦合效應。

迄今為止，中國少有同時采用CFD和DSMC方法模擬高速流動中的熱化學非平衡效應，并進行對比分析。另一方面，考慮到高速流動的地面試驗困難，綜合使用CFD和DSMC方法恰好相互驗證。本文以不同來流速度的高速圓柱繞流為研究對象，分別采用基于Park雙溫模型的CFD方法以及基于L-B內能松弛模型和Q-K化學反應模型的DSMC方法進行仿真，獲得熱化學非平衡效應對流場和氣動熱的影響，為高速氣動熱的預測提供參考。

1 模型和計算方法

1.1 物理模型和網格

采用來流工質為N2的二維高速圓柱繞流[1,2]，圓柱直徑為0.3048 m，壁溫為500 K。來流克努森數為Kn=0.01，來流速度取為Ma=5～45，其他來流參數保持恒定，具體見表1。

表1 來流條件 Tab.1 Conditions of Freestream Flow

因為圓柱后側為亞聲速流動，這將導致模擬收斂速度緩慢，而本文主要關注駐點處的壓強和熱流，所以僅對圓柱前半段進行仿真，采用貼體四邊形網格，計算域和網格結構示意如圖1所示。對于不同的來流速度工況，依據文獻[1]中的 Recell,w=1準則設定靠近壁面附近的法向網格尺寸。最大網格長寬比不超過100，具體的網格數量和網格尺寸見表2。

圖1 計算域及網格示意 Fig.1 Computational Domain and Grid Structure

表2 不同來流速度下的網格尺寸設置 Tab.2 Grid Size Setting for Different Freestream Flow Velocities

續(xù)表2

1.2 數值求解方法

CFD方法是在宏觀上通過數值求解N-S方程，而DSMC方法是在微觀上基于概率方法直接模擬分子的運動和碰撞。

開源程序OpenFOAM是一個高度可擴展的CFD軟件開發(fā)包，并包含了Lagrange粒子庫可實現DSMC方法。為了模擬高速流動，基于OpenFOAM框架，Casseau等[3,4]開發(fā)了適用于高速的雙溫CFD求解器hy2Foam；White等[5]發(fā)展了基于L-B過程的內能松弛過程和基于Q-K的化學反應過程dsmcFoamPlus求解器。本文同時使用這兩個求解器進行對此圓柱繞流開展數值模擬研究。

1.2.1 CFD方法

CFD 方法采用有限速率化學反應的Park T-Tv 雙溫松弛模型[6]。雙溫模型中平動-振動的能量傳輸是由Landau-Teller方程決定的[7]，而它的松弛時間采用 Millikan-White-Park模型[8]?；瘜W反應采用基于Arrhenius有限速率模型，且其速率系數取自DSMC的Q-K模型[9]。本文只考慮氮氣的離解過程 （N2+N2—2N+N2、N2+N—2N+N），未考慮電離過程。

組分的粘性系數采用冪律模型，與DSMC的變徑硬球分子的碰撞過程保持一致。熱擴散系數與粘性系數間符合Eucken關系。混合氣體的粘性系數和熱擴散系數采用Wilke混合規(guī)則。質量擴散采用Fick定律，且其二元擴散系數則遵循Gupta的曲線擬合[10]。使用基于第一激發(fā)能級的簡單諧振子模型和多電子激發(fā)態(tài)模型共同來確定氣體的熱力學參數。

不考慮壁面上的速度滑移，且認為是非催化壁。由于來流Re較小，流動屬于層流，不用考慮湍流模型。

1.2.2 DSMC方法

DSMC方法使用L-B模型模擬各內能態(tài)的激發(fā)與松弛，與CFD的雙溫模型中假設分子的平動和轉動能間是平衡相比，DSMC將平動能、轉動能和振動能都獨立分開。N2離解反應過程采用Q-K模型，設定VHS分子碰撞模型與CFD方法保持一致，并使用NTC方法模擬分子的運動和碰撞過程。除采用與CFD相同的網格外，DSMC算例還保證最少網格粒子數不少于7，時間步長小于最小平均碰撞時間。為了減小DSMC的統(tǒng)計誤差，采樣時間取為105倍的平均碰撞時間。

本文除了采用上述的可模擬熱化學非平衡過程的DSMC方法外，還采用忽略離解反應的DSMC方法進行對比。結果中使用的簡寫記號的含義見表3，其中Tt、Tr、Tv分別表示在DSMC中獨立考慮的氣體分子平動、轉動和振動的溫度。而CFD采用Park T-Tv雙溫松弛假設，其認為平動態(tài)和轉動態(tài)、振動態(tài)和電子態(tài)分別是平衡的，即采用2個溫度Ttr和Tve表征。

表3 數值方法設置 Tab.3 Numerical Method Settings

1.3 計算方法驗證

為了驗證本文的CFD和DSMC程序的有效性，使用了美國Holden等[11]做的“CUBRC”第7次試車實驗數據進行方法驗證?！癈UBRC”實驗測量稀薄超聲速來流對雙圓錐體的氣動力和氣動熱值，是一個典型的稀薄氣體動力學實驗。

圖2為CFD和DSMC驗證計算中雙圓錐體表面壓強、熱流分布與實驗對比。雙錐體表面壓強和熱流密度計算值與實驗值相符，最大值和最小值的位置一致，準確再現了雙錐流動中的分離渦和激波-邊界層干擾。

圖2 雙錐的氣動參數分布 Fig.2 Distribution of Aerodynamic Parameters on the Double Cone

2 計算結果分析

圖3為來流速度Ma=25且考慮化學反應的DSMC流場結果。結果表明，高超聲速來流在圓柱頭部形成弓形激波，經過激波后速度快速下降，而壓強和溫度快速上升，其中壓強在壁面的駐點處達到最大值。由于壁面設定為恒溫壁面，平動、轉動和振動溫度是先增加后減小，并考慮到內能松弛效應，平動、轉動和振動溫度存在明顯差別。

圖3 Ma=25時的DSMC模擬的流場 Fig.3 Flow Field from DSMC Simulation at Ma=25

圖4為DSMC模擬的駐點線上平動、轉動、振動溫度和氮原子質量分數分布，溫度使用來流靜溫無量綱化。結果表明，平動和轉動溫度間的非平衡效應主要出現在激波位置處，而振動溫度在整個波后位置均表現出與前二者的非平衡。考慮氮氣離解反應時各模式溫度的極大值明顯降低。

圖4 Ma=25時的DSMC模擬的駐點線上的溫度分布 和氮原子質量分數分布 Fig.4 Temperature and N% Distributions on the Standing Line for the DSMC Simulation at Ma =25

圖5為各模式溫度和N質量分數的最大值隨Ma的變化，為了統(tǒng)一不同來流速度下的溫度，這里的溫度使用了來流總溫無量綱化。圖5a表明，DSMC與CFD在考慮氮氣離解時，獲得的氮原子質量分數與來流速度的關系基本相同，并且當Ma＜15時，氮氣離解是可以忽略的，此時圖5b中DSMC有無離解反應的平動溫度rT以及轉動溫度tT是重合的。由于在DSMC中考慮了轉動能的松弛，圖5b中的轉動溫度tT的分布曲線總是稍低于平動溫度rT，而平動溫度rT會出現大于考慮平動-轉動平衡的理想總溫值的現象。當Ma＞15時，隨著馬赫數的增加，氮氣離解越來越多，考慮化學反應的平動溫度rT和轉動溫度tT與無化學反應得到的溫度差距越來越大，DSMC在考慮氮氣離解時所獲得的平動溫度rT和轉動溫度tT要比CFD的下降的快。圖5c中DSMC在考慮氮氣離解時所獲得的振動溫度vT也比CFD低，隨著來流速度的增加，振動溫度vT呈現先減小后增大再減小的變化規(guī)律，這是由于振動溫度vT的最大值是受駐點線上流動特征時間、振動態(tài)的弛豫特征時間以及弛豫過程中離解反應對振動能消耗的綜合影響。

圖5 各模式溫度和N質量分數的最大值隨Ma的變化 Fig.5 Maxima of Temperture and N% with Freestream Velocity

續(xù)圖5

圖6為DSMC和CFD方法獲得的駐點處壓強系數Cp和熱流系數CH。

圖6 駐點處壓強和熱流的預測 Fig.6 Prediction of Pressure and Heat Flux at the Stationary Point

分別使用下式進行無量綱化：

式中wp為壁面壓強；U∞為來流速度；wq為壁面熱流。

由圖6可知，高速條件下，駐點壓強系數一直維持在1.80～1.85之間，呈現出馬赫數無關效應，也不隨化學反應而變化。對于熱流系數，當Ma＜15時，兩種方法預測的結果基本一致，隨馬赫數的增加而增大，并且化學反應效應不明顯；在15≤Ma＜40時，化學效應顯現，隨著來流速度的增加，考慮化學反應的熱流系數先減小后增大，且CFD和DSMC方法的規(guī)律基本一致，而忽略化學反應會造成熱流預測偏高。

對比CFD和DSMC預測的表面參數，兩者的結果基本是一致的，并表現出隨來流速度相同的變化規(guī)律，但是兩者對流場的預測出現一定偏差，需要進一步細致研究兩者在物理模擬方法的差異。

3 結 論

本文采用可描述熱化學非平衡過程的CFD和DSMC方法對不同來流速度下的高速圓柱繞流進行了數值模擬，獲得了駐點線上的各模式溫度、駐點氣動參數隨來流速度的變化關系，并討論了熱化學非平衡效應對氣動熱的影響，得到了如下結論：

a）CFD和DSMC方法分別在宏觀和微觀上發(fā)展了模擬熱化學非平衡的模型，在相同的熱物性設定基礎上可以獲得基本一致的氣動參數分布；

b）在稀薄來流條件下，分子內能的激發(fā)與離解過程存在松弛效應，并且隨著來流速度增大，熱化學非平衡效應對氣動熱的準確預測影響越顯著；

c）在充分考慮熱化學非平衡過程時，熱流系數隨來流速度的增加呈現先增后減再增的變化規(guī)律，而忽略化學反應會造成熱流預測偏高。