高超聲速三維熱化學非平衡流場的數值計算對比研究

張敏捷，向樹紅

（1. 北京衛星環境工程研究所 可靠性與環境工程技術重點實驗室；2. 北京衛星環境工程研究所：北京 100094）

?

高超聲速三維熱化學非平衡流場的數值計算對比研究

張敏捷1, 2，向樹紅1, 2

（1. 北京衛星環境工程研究所 可靠性與環境工程技術重點實驗室；2. 北京衛星環境工程研究所：北京 100094）

摘要：文章選取國內外典型的高超聲速繞流實驗，對比研究了不同化學反應模型（Dunn-Kang模型、Park85模型、Park93模型、Park2001模型、Gupta模型）對高超流動計算結果的影響。重點采用Park雙溫度模型開展熱化學非平衡效應的模擬研究，并同時與5組元單溫度模型的計算結果進行了對比。研究表明：Park85-7組元雙溫度模型與實驗結果吻合最好，結果最為可靠；而Park85-5組元單溫度模型更適用于工程計算。

關鍵詞：高超聲速繞流；熱化學非平衡；計算流體力學

隨著飛行器速度和高度的增加，飛行器周圍空氣經歷了從量熱完全氣體、熱完全氣體到高溫氣體的變化，其熱化學狀態由最初的常溫無反應氣體逐步過渡到單溫度描述的化學平衡態、化學非平衡態以及最終的多溫度熱化學非平衡態，與化學反應有關的氣體組分也由低速時期的2組元（O2、N2）逐步變化為高速時的5組元（O2、N2、O、N、NO）、

0 引言

自20世紀80—90年代起，各航天大國開始積極研制下一代低成本、可重復使用的天地往返運輸系統，掀起了新一輪的高超聲速非平衡流研究熱潮。

高超聲速與高溫流動現象一般發生于再入飛行器重返大氣層期間，包括：彈道式導彈、飛船返回艙、航天飛機、返回衛星、空天飛機等，處于高超聲速飛行階段的各類驗證機、導彈等，以及行星7組元（O2、N2、O、N、NO、NO+、e-）甚至11組元（O2、N2、O、N、NO、O2+、N2+、O+、N+、NO+、e-）。飛行器在發射和再入過程中，都不同程度的經歷了各種熱化學態區域。研究不同熱化學性質的空氣繞流流場特性及其相對應的準確、高效的數值模擬方法，在航天飛行器設計和研制中是十分重要的。

目前一般通過求解耦合了空氣熱化學特性的Navier-Stokes方程來計算模擬高超聲速流動。針對熱化學非平衡流動的數值模擬，國內外研究人員提出了多種化學反應動力學模型，有些研究人員在不同時期也提出了不同的化學反應速率模型，其中具有代表性的有Dunn-Kang模型[1]、Park85模型[2]、Park93模型[3]、Park2001模型[4]、Gupta模型[5]等。每個模型包含多種反應式集合及相應的反應速率常數，不同的模型之間化學反應速率量值有時相差數個量級，但其優劣卻難以判定。由于高溫流動的復雜性，目前還沒有明確和統一的化學反應模型選取準則。

本文選取國內外典型的高超聲速試驗進行數值模擬，對比采用不同熱化學模型的計算精度和所消耗的計算資源，系統研究了適用于工程計算的高超聲速流動計算方法。

1 理論基礎

1.1高超聲速黏性流動控制方程

空氣在高超聲速高溫氣體流動中處于熱化學非平衡狀態，氣體內能模式由平動溫度（代表平動能）、轉動溫度（代表轉動能）、每種雙原子組分的振動溫度（代表振動能）和電子溫度（代表重粒子的激發能和自由電子的平動能）這4個不同的溫度來描述。對熱力非平衡流動，本文采用Park提出的雙溫度模型[6]：認為如分子或原子的重粒子平動溫度等于其轉動溫度，用單一溫度T描述；而其振動溫度等于電子的平動溫度，由單一溫度TV來描述。在非平衡條件下，內能交換是通過粒子碰撞或各種類型如原子復合、分子離解等化學方式進行的。通常，自由電子平動模式和分子振動模式的能量交換非常迅速，很快達到平衡；對于重粒子，其平動能和轉動能之間的能量交換也很快，可迅速達到平衡。由此，在直角坐標系下，三維守恒形式的熱化學非平衡Navier-Stokes方程為[7]

其中：Q為守恒變量；E、F、G為對流項；ED、FD、GD為擴散項；S為化學源項。它們表示為

式中：u、v、w分別為x、y、z方向的速度；p、ρ、e、eV分別為混合氣體的壓強、密度、單位質量氣體總能量和單位質量氣體振動?電子能量；k為單位質量氣體流體湍流動能；ρs、fs、hs、Ds分別為組分s的密度、質量分數、單位質量靜焓及質量擴散系數；qx為x方向的熱通量；qV, x為x方向的振動?電子熱通量；τxj為剪切應力張量，其中j代表x、y、z方向；sω.為各組元的化學反應質量生成率；Vω.為振動能量生成率。F、G的表達式與E相似；FD、GD的表達式與ED相似。

1.2化學反應機理

高溫氣體中N2和O2之間將出現化學反應，其中包括碰撞離解反應、置換反應、締合電離反應、碰撞電離反應、附著反應等。考慮不同程度的化學反應，空氣模型分為2組元、5組元、7組元、8組元、9組元、11組元、15組元等。在高超聲速流動研究中較為常用的空氣模型為5組元、7組元和11組元模型。

對化學反應氣體，有限速率反應方程可寫為

其中：v′sr和v″sr分別為反應物和生成物的化學當量系數；Ms為參與反應的分子式；Ns為化學反應的方程式個數；Nr為單個化學反應方程式內反應項個數。化學源項可表示為

式中：Kfr為正向化學反應速率；Kbr為逆向化學反應速率；Mw,s為反應物的相對分子質量。根據Arrhenius公式，反應速率是溫度的函數，則有

式中，C1、η和(-ε0/k)為通過實驗數據擬合出來的系數。針對高溫氣體化學反應，研究人員提出了不同的化學反應模型。Dunn-Kang在1973年提出了一種化學反應模型，總計26種反應式，用于計算化學反應速率，完整的模型數據見文獻[1]。Dang-Kang模型雖然是當前常用的，但并不是很準確的模型[8]。Chul Park于1985年提出了一套化學反應式（總計21種反應式）以及相應的正向反應速率計算常數，完整數據見文獻[2]。Park在1993年對模型增加了幾個反應式，并對計算常數做了修改，完整數據見文獻[3]。Park于2001再對化學模型進行了總結修改，并考慮了燒蝕壁面的反應物，完整數據見文獻[4]。Gupta于1990年總結了Bortner[9]和Dunn-Kang的反應模型，得到的新反應模型總計20個反應式。模型中前7個反應式取自于Bortner模型，其余數據取自于Dunn-Kang模型，完整的Gupta模型數據見文獻[5]。不同化學反應模型系數不同，有的甚至相差數個量級。

逆向反應速率可用2種方法求出，第1種是通過平衡常數求出，即

其中平衡常數Keq根據溫度擬合曲線求出。第2種方法是根據Arrhenius公式求逆向反應速率得到：

其中，Cb,1、ηb和(-εb,0/k)為通過實驗數據擬合出來的系數。

平動和振動模式之間的能量交換是重要的源項，采用標準的Landau-Teller形式表達有

其中，e*V,s、eV,s、τV,s分別為組分s的平衡振動能、非平衡振動能和相應的振動松弛時間。

2 RAM-C飛行試驗

RAM-C（Radio Attenuation Measurements）是NASA蘭利研究中心進行的一項再入飛行實驗項目[10-11]。項目一共進行了3次飛行實驗，實驗飛行器為頭部半徑0.1524m、半錐角9°、長1.3m的球錐。飛行器有效飛行高度為53.3～85.3km，最大飛行速度7620m/s。在此范圍內測得了飛行器軸向電子數密度分布。RAM已經成為驗證高超聲速流動模型的經典工具，大量的研究通過對比飛行實驗測得的電子數密度數據，來驗證各種熱化學模型。本文選取了61km和71km高度的飛行數據進行數值計算，采用不同的化學模型（Park85、Park93、Park2001、Gupta、Dunn-Kang）以及雙溫度模型考（慮熱力非平衡效應）對比驗證不同化學模型的計算精度。另外，采用Park85-5組元（O2、N2、O、N、NO）單溫度模型進行了計算。試驗結果數據見表1，其中電子數密度單位為個/cm3。

表1　RAM-C飛行實驗數據：電子數密度峰值分布Table 1 RAM-C flight data: peak electron number density distribution

2.1建模

應用CAD軟件CATIA建模工具對RAM飛行器進行實體建模。用商業網格劃分軟件ICEM對實體模型進行網格劃分。劃分網格時，在壁面及激波處需加密處理。垂直壁面方向布置80個網格點，壁面第一層網格高度為2×10-5。網格量為117078。表2列出了飛行實驗自由來流條件。RAM飛行器模型及計算網格如圖1所示。

表2　詳細計算條件Table 2 Details for computation

圖1　RAM-C計算流體力學網格Fig. 1 RAM-C grid for CFD

2.2計算結果及分析

本文的計算均在聯想24核心D30工作站上進行，收斂殘差準則為10-3，同時檢測給定點的各項物理參數，保證計算收斂。

表3列出了本文計算結果與其他文獻計算結果。

表3　計算結果對比Table 3 Comparison of various computation results

對比雙溫度模型和其他文獻計算結果，在關鍵性的激波位置指標上，在61km處，文獻計算結果范圍為6.75%～9.18%，Park85化學反應模型的計算結果為7.667%；在71km處，文獻計算范圍為8.47%～12.46%，Park85為8.365%。對峰值溫度，在61km處，文獻計算范圍為18400～22500K，Park85為20940K；在71km處，文獻計算結果為19900～25000K，Park85為21920K。由對比看出，Park85模型計算結果最為準確。

對比Park85-7組元雙溫度模型和Park85-5組元單溫度模型計算結果，在61km處，激波位置分別為7.667%、7.666%；峰值溫度分別為20940K、19058K。在71km處，激波位置分別為8.365%、8.159%；峰值溫度分別為21920K、19802K。比較二者差別不大。

Park85-7組元雙溫度模型計算出的駐點線上各組分質量分數分布如圖2所示。

圖2　沿駐點線質量分數分布Fig. 2 Mass fraction distribution along stagnation line

從圖2可看出：在61km處，由于激波的存在，各組分質量分數在激波處變化非常劇烈。在71km處，空氣更加稀薄，化學非平衡效應更加明顯，激波厚度更大，激波處各組分質量分數變化較為緩慢。

工程上更為關心飛行器壁面的熱流密度分布。表4列出了71km處本文熱流密度峰值計算結果與其他文獻的計算結果。總體而言，文獻的計算結果范圍為0.58～1.02MW/m2，且集中在1MW/m2處，與Park85模型計算結果吻合較好。

表4　駐點熱流密度Table 4 Stagnation heat flux

圖3為Park85-7組元雙溫度計算結果與NASA實驗數據的對比，可看出二者吻合的較好。圖4為不同化學模型計算出的駐點線溫度分布。

圖3　電子數峰值密度Fig. 3 Peak electron number density compared with experiment data

圖4　駐點線溫度分布Fig. 4 Temperature distribution along the stagnation line

由以上可以看出，Park85-7組元雙溫度化學模型能較好地捕捉流動特征，計算結果較為準確。另一方面，5組元化學模型所包含的化學反應方程式相對較少（總共17個化學反應方程式），且只有1個特征溫度，計算量相比雙溫度模型要小得多。圖4給出了61km和71km高度下Park85-5組元單溫度和7組元雙溫度模型的計算對比結果。由圖可看出，在61km處，對比單溫度（熱力平衡）模型和雙溫度（熱力非平衡）模型，計算出的平動溫度T重合度較高；對雙溫度模型，計算出的氣體平動溫度和振動溫度分離相差不遠，化學非平衡現象不明顯。而在71km處，2種計算模型計算所得平動溫度相差較大；對雙溫度模型，計算出的氣體平動溫度與振動溫度亦相差較多，化學非平衡現象較為明顯。這與實際物理規律是相符合的：高度越高，氣體越稀薄，化學非平衡現象越明顯。

另一方面，相比5組元單溫度模型，7組元雙溫度模型由于增加了一個用于描述氣體振動?電子能量的溫度TV，增加了所需求解的偏微分方程的個數，并且7組元模型所包含的化學反應方程式較多（總共24個），計算量大大增加，這在工程上是難以接受的。目前關于單溫度熱力平衡流動的計算技術已經發展得非常成熟，而多溫度模型在不同能量模態間的能量交換等方面尚存在不少值得研究的地方。對此，本文給出的建議是，對于高超聲速流動的工程計算，在計算高度不是特別高，速度不是特別快時，一般采用5組元單溫度（熱力平衡化學非平衡）模型進行計算，可以兼顧計算效率和計算精度。

3 球繞流驗證實驗

中國空氣動力研究與發展中心超高速所在二級輕氣炮（彈道靶）中對處于高超聲速非平衡態下的球模型激波脫體距離進行準確測量[19]，為非平衡流的理論和計算研究提供用于計算程序驗證的、可靠的實驗數據。所用實驗模型為球模型，直徑12mm。為了進一步對比驗證Park85-5組元單溫度模型和7組元雙溫度模型的計算精度，對這一實驗進行了仿真分析。實驗條件及測量結果如表5所示。

表5　球模型脫體激波測量結果Table 5 Experimental results of shockwave distance for sphere model

球模型及計算網格如圖5所示，在垂直壁面附近及激波附近對網格進行了加密，垂直壁面方向布置120個網格點，壁面第一層網格高度為2×10-5。網格量為225703。

圖5　球模型計算網格Fig. 5 CFD grid for sphere model

根據實驗條件取等溫壁面600K，湍流模型簡化為層流模型，收斂準則取為10-4。經過約2萬次迭代后，計算結果收斂。值得注意的是：物理上，激波為極薄的一小段厚度，在這個厚度內流動參數發生跳躍式突變；而數學上，流動參數變化是連續的。本文取氣體平動溫度峰值處作為激波所在位置。表6列出了計算結果與實驗結果。可看出7組元雙溫度模型計算結果與實驗結果吻合較好，5組元單溫度模型的結果與實驗結果相差也不大，但其消耗的計算時間卻大大地減少。

表6　計算結果和實驗結果的對比Table 6 Comparison between CFD results and experiment data

4 結束語

本文針對國內外典型的高超聲速流動實驗數據，采用不同的化學反應模型（Park85、Park93、Park2001、Gupta、Dunn-Kang）開展了對比研究，結果表明：Park85-7組元雙溫度模型具有較好的計算精度，能較為精確地捕捉流動特征，所得計算結果與實驗結果吻合很好。在此基礎上，采用Park85-5組元單溫度模型進行了計算。對比分析結果表明，Park85-5組元單溫度模型計算量小，計算精度與7組元雙溫度模型相比，相差不大，適合于工程計算。對此，本文給出的建議是：對于高超聲速流動的工程計算，在計算高度不是特別高，計算的飛行速度不是特別快時，一般采用Park85-5組元單溫度模型進行計算，可以兼顧計算效率和計算精度。

參考文獻（References）

[1] Dunn M G, Kang S W. Theoretical and experimental studies of reentry plasma[R]. NASA Contractor report, 1973

[2] Park C. On convergence of computation of chemistry reacting flows[C]//AlAA 23rdAerospace Sciences Meeting. Reno, Nevada,1985

[3] Park C. Review of chemical-kinetic problems of future NASA missions I: Earth entries[J]. Journal Of Thermophysics and Heat Transfer, 1993, 7(3): 385-398

[4] Park C. Chemical-kinetic parameters of hyperbolic earth entry[J]. Journal Of Thermophysics and Heat Transfer, 2001, 15(1): 76-90

[5] Gupta R N, Jerrold M Y, Richard A T, et al. A review of reaction rates and thermodynamic and transport properties for an 11-species air model for chemical and thermal nonequilibrium calculations to 30000 K[R], NASA Reference Publication 1232, 1990

[6] Park C. Assessment of two temperature kinetic model for ionizing air[C]//AIAA 22ndThermophysics Conference. Honolulu, Hawaii, 1987

[7] Anderson J D. Hypersonic and high-temperature gas dynamics[M]. AIAA Education Series, 2006

[8] 卞蔭貴，徐立功. 氣動熱力學[M]. 合肥: 中國科學技術大學出版社, 2011: 114-116

[9] Bortner M H. Suggested standard chemical kinetics for flow field calculations: a consensus opinion[C]//AMRAC proceedings. University of Michigan, 1966-04-18: 569-581

[10] Jones W L, Cross A E. Electrostatic-probe measurements of plasma parameters for two entry flight experiments at 25000 feet per second, NASA TN D-6617[R], 1972

[11] Akey N D, Cross A E. Radio blackout alleviation and plasma diagnostic results from a 25000 foot per second blunt-body reentry, NASA TN D-5615[R], 1970

[12] Candler G V, MacCormack R W. Computation of weakly ionized hypersonic flows in thermochemical nonequilibrium[J]. Journal of Thermophysics, 1991, 5(3): 266-273

[13] Coquel F，Marmignon C. Simulation numérique d'écoulements endéséquilibre thermochimique et faiblement ionizes[R]. ONERA Rapport Technique de Synthése TRS 14/1929 AY. F, 1929

[14] Josyula E, Bailey F W. Governing equations for weakly ionized plasma flowfields of aerospace vehicles[J]. Journal of Spacecraft and Rockets, 2003, 32(6): 845-857

[15] Tchuen G, Zeitoun D E. Computation of weakly ionized air flow in thermochemical nonequilibrium over sphere-cones[J]. International Journal of Heat and Fluid Flow, 2008, 29: 1393-1401

[16] Tchuen G, Zeitoun D E. Effects of chemistry in nonequilibrium hypersonic flow around blunt bodies[J]. Journal of Thermophysics and Heat Transfer, 2009, 23(3): 433-442

[17] Walpot L M. Development and application of a hypersonic flow solver[D]. PhD Thesis, TU Delft, 2002

[18] Soubrie T, Rouzaud O, Hylkema J. Computation of weakly multi-ionized gases for atmospheric entryusing an extended roe scheme[C]. European Community on Computational Methods in Applied Sciences. Finland, Jyv?skyl?, 2004: 805-810

[19] 柳軍. 熱化學非平衡流及其輻射現象的實驗和數值計算研究[D]. 長沙: 國防科學技術大學, 2004

（編輯：肖福根）

A comparative study of the computation of 3-D hypersonic flow in thermochemical nonequilibrium state

Zhang Minjie1, 2, Xiang Shuhong1, 2
(1. Science and Technology on Reliability and Environment Engineering Laboratory, Beijing Institute of Spacecraft Environment Engineering;
2. Beijing Institute of Spacecraft Environment Engineering: Beijing 100094, China)

Abstract:The influence of the chemical kinetic model selection on the computation of the hypersonic flow in the thermochemical nonequilibrium state is examined, with consideration of several different chemical reaction models (the Dunn-Kang model, the Park85 model, the Park93 model, the Park2001 model, and the Gupta model). Typical experimental results at home and aboard are used to identify the reliable model. The thermal nonequilibrium processes in the gas are accounted in the Park’s two temperature model. Cases are computed using five species, and one temperature model is used for comparison. It is shown that the Park85-7 species-two temperature model gives the best results, as compared to experiment data. On the other hand, the Park85-5 species-one temperature model is more suitable for engineering applications.

Keywords:hypersonic flow; thermochemical nonequilibrium; CFD著陸探測器。

作者簡介：張敏捷（1991—），男，碩士研究生，研究方向為高超聲速計算流體力學、航天器噪聲分析、工程振動；E-mail: spacezmj@163.com。指導教師：向樹紅（1963—），男，研究員，博士生導師，主要從事航天器噪聲分析、控制和試驗、爆炸沖擊機理、防護和試驗、振動試驗方法、極端力熱等綜合環境分析、空間環境觀測搭載等技術的研究；E-mail: Xshxsh638@sohu.com。

收稿日期：2015-09-11；修回日期：2016-01-06

DOI:10.3969/j.issn.1673-1379.2016.01.006

中圖分類號：V434+.11

文獻標志碼：A

文章編號：1673-1379(2016)01-0035-07