0 引言
太空探测器返回舱再入大气层过程中，速度可达第二宇宙速度11.2km/s，返回舱周围流场最高温度超过10000K，高温气体辐射加热在总气动加热中的比重将明显增大。此时，除了考虑边界层向物面的对流加热(包括热传导和扩散传热)外，需要考虑高温流场气体辐射加热对返回舱热环境的影响[1, 2]。高温气体辐射加热与热传导等其它传热方式在机理上有显著不同，前者以电磁波或光子作为载体，后者则是以分子作为载体。热传导与温度梯度近似成正比关系，而热辐射则与绝对温度直接相关。
高温流场气体辐射加热机理复杂，涉及高温气体动力学、原子分子光谱学、传热学等诸多学科的交叉与融合。对辐射加热的可靠预测不但与计算方法和计算模型本身有关，还与方法和模型中用到的大量基础数据如高温气体组分光谱数据和化学反应动力学数据密切相关，目前对辐射加热计算的精准度和通用性还不能完全满足工程要求，有关物理化学模型和计算方法及相关基础数据还在发展和完善之中[3, 4, 5, 6, 7, 8]。高温流场气体辐射加热计算方法通常有两种，即耦合求解方法和非耦合求解方法[9, 10, 11, 12, 13, 14, 15, 16, 17]。耦合求解方法考虑流场和辐射场之间的耦合干扰效应，把辐射传输方程和包含辐射源项的非平衡N-S方程耦合起来进行求解。非耦合方法则是通过数值求解不含辐射源项N-S方程对流场进行数值模拟，再基于数值模拟获得的流场参数分布，通过求解辐射传输方程并对物面某点立体角范围的所有流体微元辐射进行积分获得该点的辐射加热。耦合求解方法计算精度高，但耦合求解方法计算量很大，对计算资源的要求很高，目前只用于简单流动气体的辐射加热计算。在非耦合求解方法中，没有考虑高温气体辐射对流场产生的冷却效应，故采用此方法获得的物面辐射热流比耦合求解方法偏高，但此方法计算量相对较小，相比一般工程方法的通用性强，作为一种保守的预测方法在辐射加热工程分析中获得了很好应用[2, 11, 12, 13, 14]。
针对高温流场气体辐射加热问题，国外已经发展形成相对成熟的高温非平衡/平衡流场与流场气体辐射耦合/非耦合数值分析软件，并在航天工程中获得广泛应用。如美国NASA研制的LAURA/LORAN与DPLR/NEQAIR流场辐射耦合计算代码[9, 11, 12, 13, 14, 15, 16]、俄罗斯科学院研制的NERAT+ASTEROID计算代码[7]、欧洲有关研究机构研制的URANUS/PARADE计算代码[17]。国内在此方面的发展相对滞后，目前采用的预测方法主要是一些近似的工程方法，模拟辐射气体动力学耦合效应数值分析代码还不够成熟。
本文发展了一套复杂外形高温绕流流场与气体辐射非耦合气动加热计算分析代码。采用此计算代码，首先对轴对称外形驻点辐射加热及对流加热效应进行研究，比较分析模型与算法的可行性，然后针对典型再入条件，研究返回舱周围高温流场气体辐射加热对返回舱热环境的影响。 1 计算方法 1.1 流场数值方法
一般曲线坐标系下，雷诺平均的无量纲化三维N-S方程为：
其中：守恒变量Q=J－1(ρi,ρ,ρu,ρv,ρw,ρE)T，W=J－1(wi,0,0,0,0,0)T，J－1是坐标变换系数，ρi是组分i的密度，u、v、w对应迪卡尔坐标下三个方向的速度，E为总能，Re是雷诺数，F、G、H与FV、GV、HV分别对应不同方向的无粘与粘性通量，为化学反应非平衡源项。
采用LU-SGS全隐式有限差分方法对方程(1)进行离散求解，其中无粘项采用AUSMPW+格式离散，粘性项采用中心差分格式离散，化学反应模型采用11组分的Dunn-Kang模型(O2,N2,NO,O,N,NO+,O2+,N2+,O+,N+,e)，详见文献[18, 19]。一般对流加热方式下物面热流表达式为：
其中：k为热传导系数，Di，hi和ci为组分i的扩散系数、焓与质量分数，T为温度，n为物面法向坐标。 1.2 辐射计算方法
热力学非平衡条件下辐射传输方程：
式(3)中Jν和kν为沿光学路径l的气体光谱发射和吸收系数。单位面积单位时间的物面法向辐射加热通量为：
其中Ω、ν分别为物面法向立体角和光的频率，设球坐标系下光学路径方向为球坐标矢径r的方向，它与物面法线方向夹角为θ，在光学薄的假设下，由(3)与式(4)得：
对方程(5)进行数值积分可得：
式中积分上下限νi、νf为光谱频率上、下限值，φi、φf为方位角上、下限值，θi、θf为极距角上、下限值，L为矢径r的积分长度，光谱发射系数Jν采用考虑高温空气非平衡效应的线-线光谱辐射模型，本文在计算中主要考虑如下高温空气辐射机制[16, 20]：
1) NO的β带系(B2Π→Χ2Π)；
2) NO的γ带系(A2Π→Χ2Π)；
3) O2的S－R带系(B3∑－u→X3∑－g)；
4) N2的第一正带系(B3Πg→A3∑+u)；
5) N2的第二正带系(C3Πu→B3Πg)；
6) N2+第一负带系(B2Π+u→X3∑+g)；
7) O和N原子的线状谱(Bound-Bound跃迁)；
8) O和N原子的连续谱(Bound-Free跃迁)；
9) O+和N+离子的连续谱(Free-Free跃迁)。
通过式(6)计算高温气体流场对物面的辐射加热时，需要以物面点为坐标原点建立当地球坐标系下的光学计算网格，并把直角坐标系下的流场网格结点的组分、温度等参数插值到光学计算网格结点上。其次，针对高温气气体流场各种化学组分的不同辐射机制，计算光学网格点的气体的发射系数，把每个组分不同辐射机制下的发射系数进行积分获得每个网格结点的气体总的发射系数。 2 计算结果分析
(1) 数值和工程计算结果的比较
Tauber与Sutton等人[21]通过研究地球大气层再入和火星大气层再入条件，提出了再入飞行器驻点辐射加热工程计算公式。公式(7)适用于速度范围：V∞≤9km/s，Johnson等人[22]在此工作基础上对式(7)进行拓展，形成了更高速度范围的辐射热流计算公式(8)，其适用范围9km≤V∞≤16km/s，这两个公式结合在一起，在实际工程设计中得到很好应用[1, 2, 10]。
式中：ρ0为参考密度，Rn为飞行器头部半径，f(V∞)、k1、k2和k3均为来流速度V∞的函数，c与b为常数。
算例1：外形为半径为1.22m的半球，计算高度58km，来流速度在6km/s~15.24km/s之间。流场计算采用11组分的化学模型。图 1给出非催化壁、速度15.24km/s条件下半球绕流温度和原子及其离子参数分布，可见在此条件下波后空气发生大量离解和电离，生成大量的O与N原子及O+和N+离子。图 2给出用Gupta[23]和Martin[24]以及Tauber等人[21]工程估算公式获得的驻点辐射热流与本文计算值的对比，可见Gupta工程计算结果明显偏高，而Martin计算结果在高速条件下偏低，在较低速度条件下与Tauber等人计算结果趋于一致，总体上看，本文数值计算结果与公式(7)和式(8)的工程预测结果较为一致。
 |
| 图 1 温度和组分质量分数云图(V∞=15.24km/s)Fig. 1 Contours for temperature and species mass fraction (V∞=15.24km/s) |
 |
| 图 2 数值和工程计算的驻点辐射热流比较Fig. 2 Radiative heat flux at stagnation-point |
算例2：类联盟号飞船返回舱外形。返回舱大钝头球冠半径取1.25m，其参考再入弹道类似于BLAST返回舱再入弹道[2]，如图 3所示。为了和上述Tauber等人的工程估算结果进行对比，假定飞行迎角为0°。图 4为返回舱驻点对流和辐射热流随再入时间的变化情况，图中曲线为工程估算方法得到的结果，离散点为第一次再入峰值热流附近点H=66km、63km的数值计算结果，对应飞行速度分别为10.3km/s、9.8km/s。从图 4可见，在接近辐射热流峰值的66km高度，本文辐射热流数值计算值相对工程估算值偏低，而在63km二者符合较好。在63km高度，同时给出完全催化壁(FCW)和完全非催化壁(NCW)条件下的对流与辐射热流计算结果，可见不同催化条件下的辐射热流相差很小，但对流热流值相 差较大，表明物面催化条件对辐射加热影响很小，这与文献[13, 15]关于FIRE II辐射加热的数值分析结果一致。从理论上分析，由于物面催化效应主要影响物面附近的原子和离子的复合过程，而这些化学过程对流场温度分布总体上影响不大，而流场气体温度对原子和分子电子能级的激发/反激发过程起主导作用，导致物面催化效应对辐射加热的影响不大。从图 4还可以看出，在63km热流峰值高度，采用Sutton-Grave、Detra-Kempt-Riddell 和Scott等工程预测方法[2, 25]获得的对流热流值介于完全催化壁(FCW)和完全非催化壁(NCW)条件下的数值模拟值之间，符合物面催化效应对气动加热的影响规律。总体上看，本文数值计算的物面对流与辐射热流与工程估算结果具有较好的一致性。
 |
| 图 3 参考再入弹道Fig. 3 Reference reentry trajectory |
 |
| 图 4 驻点热流沿弹道变化Fig. 4 Stagnation-point heat flux along reference reentry trajectory |
表 1给出了流场气体组分不同辐射机制对驻点辐射加热热流的贡献，由于在此高度条件下，飞行速度很高，返回舱处于第一次再入的热环境严酷区，此时流场温度很高，大量氧气和氮气分子离解成浓度很高的O和N，它们在高温条件下产生自由-自由与束缚-自由跃迁而成为高温流场气体的重要辐射机制。另外，N2第一正带系与N2+的第一负带系对总辐射热流的贡献也比较大。
| 组分辐射机制 | Qr(66km) | Qr(63km) |
| N(B-B、B-F) | 248.182 | 198.442 |
| N+(F-F) | 1.947 | 1.061 |
| O(B-B、B-F) | 97.710 | 86.510 |
| O+(F-F) | 3.648 | 2.606 |
| N2+(1-) | 11.863 | 12.075 |
| N2(1+) | 45.038 | 43.759 |
| N2(2+) | 0.111 | 0.255 |
| NO(β) | 1.354 | 1.558 |
| NO(γ) | 0.030 | 0.036 |
| O2(S-R) | 0.005 | 0.014 |
| Total | 409.888 | 346.316 |
(2) 飞船返回舱辐射加热计算分析
下面重点分析类联盟号飞船返回舱在第一次再入峰值加热弹道点附近H=66km和63km的辐射加热效应，球冠半径取1.25m，飞行迎角-23°。
图 5给出H=63km时返回舱绕流温度与组分的质量分数分布。由于此时流场温度很高，最高温度超过15000K，流场中O2和大底附近流场区域中的N2分子几乎完全离解成O和N原子，它们在流场中的质量分数很高。
 |
| 图 5 返回舱绕流温度与组分质量分数分布(H=63km)Fig. 5 Contours of temperature and species mass fraction over a reentry capsule (H=63km) |
图 6和图 7给出了H=63km和H=66km返回舱对称面上迎风与背风线辐射加热与对流热流分布。从图 6和图 7可以看出，迎风面的对流与辐射热流明显高于背风面的对流与辐射热流。迎风面辐射热流峰值出现在迎风面大底驻点附近，这与对流加热情况不同，对流热流峰值则出现在迎风面肩部靠前区域，辐射热流可达总热流的10%左右(完全催化壁)，在H=66km时最大超过10%，在H=63km时则略低一些；对于非催化物面，这个比例会稍高一些。由于实际物面是有限催化条件，此时对流加热热流介于完全催化与非催化壁之间，而物面催化效应对辐射热流影响很小，因此可以预测，实际情况下辐射加热的热流占总热流的比例将超过10%。
 |
| 图 6 迎风与背风物面热流分布(H=63km，FCW)Fig. 6 Heat flux along the upward and leeward wall (H=63km，FCW) |
 |
| 图 7 对称面迎风与背风物面热流分布(H=66km)Fig. 7 Heat flux along the upward and leeward wall(H=66km) |
为了分析返回舱外形尺寸对高温流场气体辐射加热的影响，按比例把返回舱的外形尺寸增大，大底球冠半径由1.25m增加到3m，飞行迎角为0°或-23°。图 8给出飞行迎角0°时物面驻点热流随再入时间的分布，其中曲线为工程预测结果，离散点为峰值加热时刻辐射热流的数值结果。可见，球冠尺寸对返回舱绕流气体辐射加热影响较大，在峰值加热时刻，球冠半径为3m的驻点辐射热流占总热流比例达30％。对于有迎角的飞行状态，图 9给出66km时不同尺寸的返回舱对称面迎风线上的热流分布，可以看出，除靠近肩部附近区域，随着球冠半径从1.25m增加到3m，辐射热流相比对流热流明显增大。在同一再入飞行条件下，物面对流热流随飞行器球头半径增加而 减小，而辐射热流随着球头半径增加而增加(参见式(7)或式(8))，由此就不难理解再入返回舱绕流高温气体辐射加热对物面热流的影响规律。
 |
| 图 8 不同球冠半径的驻点热流(α=0°)Fig. 8 Heat flux at stagnation-point for different sphere nose(α=0°) |
 |
| 图 9 再入返回舱迎风物面热流(α=-23°)Fig. 9 Heat flux along upward wall of a reentry capsule(α=-23°) |
从以上计算分析可得到以下初步结论：
(1) 建立了多组分高温气体非平衡流场及其辐射加热的数值计算方法和计算代码。通过比较分析典型状态下的数值与工程计算的辐射热流，初步验证了数值算法和计算代码的可行性。
(2) 在返回舱再入热环境严酷区，高温流场中O和N原子产生的连续谱和线状谱是辐射加热重要来源，另外，电子带系N2(1+)和N2+(1-)对总辐射热流的贡献也不能忽略；物面材料催化效应严重影响高温流场气体边界层的对流热流，但对辐射热流影响很小。
(3) 计算分析了典型再入条件下返回舱绕流对物面的辐射和对流加热特性。计算分析表明，在第一次再入峰值热流时刻附近，物面辐射热流占总热流比例随着返回舱外形尺寸增大而增大，在某些情况下辐射热流占总热流的比例高达30％。返回舱辐射加热最严重的区域在返回舱大底驻点附近物面区域，而对流加热最严重的区域在肩部附近靠前区域。
由于辐射气体动力学干扰问题的复杂性，对高超声速飞行器气动热环境的准确预测还存在较大困难。辐射加热预测精准度依赖于流场和辐射场计算模型、方法以及所使用的化学动力学数据、原子、分子光谱数据等基础数据的可靠性。对于再入速度达第二宇宙速度的情况，辐射气体动力学干扰效应很强，对高温绕流气体辐射耦合效应的数值模拟将面临更大挑战。在建立了非耦合的辐射加热预测分析工具的基础上，下一步将重点发展高温流场气体辐射耦合效应的数值模拟手段。
| [1] | Anderson J D. Hypersonic and high temperature gas dynamics[M]. Second Edition, American Institute of Aeronautics and Astronautics, Inc. 2006: 773-785. |
| [2] | Clemente M D, Ferrarella D. Convective and radiative heat flux estimation on a re-entry capsule[R]. AIAA-2010-5067. |
| [3] | Shang J S, Surzhikov S T. Nonequilibrium radiative hypersonic flow simulation[J]. Progress in Aerospace Sciences, 2012, 53: 46-65. |
| [4] | Hash D, Olejniczak J, et al. FIRE II calculations for hypersonic nonequilibrium aerothermodynamics code verification: DPLR, LAURA and US3D[R]. AIAA-2007-605. |
| [5] | Surzhikov S T, Shang J S. Numerical rebuilding of fire-II flight data with the use of different physical-chemical kinetics and radiation models[R]. AIAA-2013-0190. |
| [6] | Olynick D R, Henline W D, et al. Comparisions of coupled radiative Navier-Stokes flow solutions with the project Fire II flight data[R]. AIAA-94-1955. |
| [7] | Djadkin A, Beloshitsky A, Shuvalov M, et al. Uncertainties in heating predictions of segmental-conical space vehicle resulting from data on chemicaland physical kinetics[R]. AIAA-2013-1056. |
| [8] | Brandis A M, Huo W M, Johnston C O. Validation of Hyper-Rad for earth entries[R]. AIAA-2013-2777. |
| [9] | Nelson H F. Radiative heating in scramjet combustors[R]. AIAA-1996-1843. |
| [10] | Gollan R J, Jacobs P A, et al. Numerival modelling of radiating superorbital flows[J]. ANZIAM Journal, 2004, 45(E): C248-C268. |
| [11] | Hartung L C, Mitcheltree R A, Gnoffo P A. Coupled radiation effects in thermochemical nonequilibrium shock-capturing flowfield calculation[R]. AIAA-92-2868. |
| [12] | Johnston C O. Nonequilibrium shock-layer radiative heating for Earth and Titan entry[D]. Virginia Polytechnic Institute and State University, 2007. |
| [13] | Scalabrin L C, Boyd I D. Numerical simulations of the FIRE-II convective and radiative heating rates[R]. AIAA-2007-4044. |
| [14] | Hartung L C. Predicting radiative heat transfer in thermochemical nonequilibrium flow fields[R]. Theory and user's manual for the LORAN code[R]. NASA TM-4564, 1994. |
| [15] | Palmer G E, White T, Pace A. Direct coupling of the NEQAIR radiation and DPLR CFD codes[R]. AIAA-2010-5051. |
| [16] | Park C. Nonequilibrium air radiation(NEQAIR) program: user's manual[R]. NASA TM-86707, 1985. |
| [17] | Liebhart H, Wernitz R, et al. Advances for radiation modeling for Earth Re-entry in PARADE: Application to the STARDUST atmospheric entry[R]. AIAA-2012-3196. |
| [18] | 董维中. 热化学非平衡效应对高超声速流动影响的数值计算与分析[D]. 北京: 北京航空航天大学, 1996. |
| [19] | 高铁锁, 李椿萱, 等. 高超声速电离绕流的数值模拟[J]. 空气动力学学报, 2002, 20(2): 184-189. |
| [20] | 乐嘉陵. 再入物理[M]. 北京: 国防工业出版社, 2005. |
| [21] | Tauber M E, et al. Stagnation-point radiative heating relations for Earth and Mars entries[J]. Journal of Spacecraft and Rockets, 1991, 28(1):40-42. |
| [22] | Johnson J E, Starkey R P, Lewis M I. Aerothermodynamic optimization of rentry heat shield shapes for a crew exploration vehicles[J]. Journal of Spacecraft and Rockets, 2007, 44(4): 849-859. |
| [23] | Gupta R N, et al. Stagnation-point heat-transfer rate prediction at aeroassist flight conditions[R]. NASA TP-3208, 1992. |
| [24] | Martin J J. Atmospheric Re-entry, An introduction to its science and engineering[M]. Prentice-Hall, Englewood Cliffs, 1966. |
| [25] | Sutton K, Graves R A. A general stagnation-point convective-heating eqution of arbitary gas mixture[R]. NASA TR R-376,1971. |