任 翔，袁軍婭，蔡國(guó)飆

（北京航空航天大學(xué)，北京，100191）

0 引 言

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

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

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

1 模型和計(jì)算方法

1.1 物理模型和網(wǎng)格

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

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

因?yàn)閳A柱后側(cè)為亞聲速流動(dòng)，這將導(dǎo)致模擬收斂速度緩慢，而本文主要關(guān)注駐點(diǎn)處的壓強(qiáng)和熱流，所以僅對(duì)圓柱前半段進(jìn)行仿真，采用貼體四邊形網(wǎng)格，計(jì)算域和網(wǎng)格結(jié)構(gòu)示意如圖1所示。對(duì)于不同的來流速度工況，依據(jù)文獻(xiàn)[1]中的 Recell,w=1準(zhǔn)則設(shè)定靠近壁面附近的法向網(wǎng)格尺寸。最大網(wǎng)格長(zhǎng)寬比不超過100，具體的網(wǎng)格數(shù)量和網(wǎng)格尺寸見表2。

圖1 計(jì)算域及網(wǎng)格示意 Fig.1 Computational Domain and Grid Structure

表2 不同來流速度下的網(wǎng)格尺寸設(shè)置 Tab.2 Grid Size Setting for Different Freestream Flow Velocities

續(xù)表2

1.2 數(shù)值求解方法

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

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

1.2.1 CFD方法

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

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

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

1.2.2 DSMC方法

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

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

表3 數(shù)值方法設(shè)置 Tab.3 Numerical Method Settings

1.3 計(jì)算方法驗(yàn)證

為了驗(yàn)證本文的CFD和DSMC程序的有效性，使用了美國(guó)Holden等[11]做的“CUBRC”第7次試車實(shí)驗(yàn)數(shù)據(jù)進(jìn)行方法驗(yàn)證?！癈UBRC”實(shí)驗(yàn)測(cè)量稀薄超聲速來流對(duì)雙圓錐體的氣動(dòng)力和氣動(dòng)熱值，是一個(gè)典型的稀薄氣體動(dòng)力學(xué)實(shí)驗(yàn)。

圖2為CFD和DSMC驗(yàn)證計(jì)算中雙圓錐體表面壓強(qiáng)、熱流分布與實(shí)驗(yàn)對(duì)比。雙錐體表面壓強(qiáng)和熱流密度計(jì)算值與實(shí)驗(yàn)值相符，最大值和最小值的位置一致，準(zhǔn)確再現(xiàn)了雙錐流動(dòng)中的分離渦和激波-邊界層干擾。

圖2 雙錐的氣動(dòng)參數(shù)分布 Fig.2 Distribution of Aerodynamic Parameters on the Double Cone

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

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

圖3 Ma=25時(shí)的DSMC模擬的流場(chǎng) Fig.3 Flow Field from DSMC Simulation at Ma=25

圖4為DSMC模擬的駐點(diǎn)線上平動(dòng)、轉(zhuǎn)動(dòng)、振動(dòng)溫度和氮原子質(zhì)量分?jǐn)?shù)分布，溫度使用來流靜溫?zé)o量綱化。結(jié)果表明，平動(dòng)和轉(zhuǎn)動(dòng)溫度間的非平衡效應(yīng)主要出現(xiàn)在激波位置處，而振動(dòng)溫度在整個(gè)波后位置均表現(xiàn)出與前二者的非平衡?？紤]氮?dú)怆x解反應(yīng)時(shí)各模式溫度的極大值明顯降低。

圖4 Ma=25時(shí)的DSMC模擬的駐點(diǎn)線上的溫度分布 和氮原子質(zhì)量分?jǐn)?shù)分布 Fig.4 Temperature and N% Distributions on the Standing Line for the DSMC Simulation at Ma =25

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

圖5 各模式溫度和N質(zhì)量分?jǐn)?shù)的最大值隨Ma的變化 Fig.5 Maxima of Temperture and N% with Freestream Velocity

續(xù)圖5

圖6為DSMC和CFD方法獲得的駐點(diǎn)處壓強(qiáng)系數(shù)Cp和熱流系數(shù)CH。

圖6 駐點(diǎn)處壓強(qiáng)和熱流的預(yù)測(cè) Fig.6 Prediction of Pressure and Heat Flux at the Stationary Point

分別使用下式進(jìn)行無量綱化：

式中wp為壁面壓強(qiáng)；U∞為來流速度；wq為壁面熱流。

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

對(duì)比CFD和DSMC預(yù)測(cè)的表面參數(shù)，兩者的結(jié)果基本是一致的，并表現(xiàn)出隨來流速度相同的變化規(guī)律，但是兩者對(duì)流場(chǎng)的預(yù)測(cè)出現(xiàn)一定偏差，需要進(jìn)一步細(xì)致研究?jī)烧咴谖锢砟M方法的差異。

3 結(jié) 論

本文采用可描述熱化學(xué)非平衡過程的CFD和DSMC方法對(duì)不同來流速度下的高速圓柱繞流進(jìn)行了數(shù)值模擬，獲得了駐點(diǎn)線上的各模式溫度、駐點(diǎn)氣動(dòng)參數(shù)隨來流速度的變化關(guān)系，并討論了熱化學(xué)非平衡效應(yīng)對(duì)氣動(dòng)熱的影響，得到了如下結(jié)論：

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

b）在稀薄來流條件下，分子內(nèi)能的激發(fā)與離解過程存在松弛效應(yīng)，并且隨著來流速度增大，熱化學(xué)非平衡效應(yīng)對(duì)氣動(dòng)熱的準(zhǔn)確預(yù)測(cè)影響越顯著；

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