对于超高速流动的研究手段主要有试验方法及数值方法。膨胀管/风洞是目前世界上能够产生超高速试验流场的少数地面试验设备之一，其结构与反射型激波风洞类似，通过去除后者的驻室，在激波管下游串联一个等截面的加速段，通过主激波波后气体的非定常膨胀使试验气流进一步加速，得到更高的气流速度和总焓。膨胀管的概念是由Resler和Bloxsom[4]首次提出的，Trimpi[5]经过详细的理论研究，分析了膨胀管在产生超高速气流方面的优势，指出只要改变膨胀管各段的充气压力，就可以得到大范围的试验气流。自从20世纪60年代膨胀管研究的进入第1个高潮期，但一直没有得到稳定的试验气流。直到20世纪80年代末，Paull等[6]将自由活塞驱动与膨胀管相结合，改善了膨胀管的性能，获得了稳定的可用试验气流。随着宇宙探测、星际旅行对超高速模拟试验的需求，20世纪90年代又出现了膨胀管研制的一个新高潮。世界范围内的膨胀管/风洞设备主要集中在澳大利亚及美国，澳大利亚Queensland大学发展了自由活塞驱动X-系列膨胀管，最大速度可达13 km/s[7]，但是试验时间很短。美国CURBC研究中心在验证型膨胀风洞LENS-X的基础上得到了尺寸放大版膨胀风洞LENS-XX，总长超过70 m，试验时间可达毫秒量级，是世界上最大的膨胀管/风洞设备[8]。国内目前已知的膨胀管设备，是中国科学院力学研究所高温气体动力学国家重点实验室(LHD)于2008年建立的爆轰驱动膨胀管JF-16，可以实现速度超过10 km/s的超高速试验气流，试验时间为100 μs，为超高速地面试验提供了基础支撑条件[9]。
然而，由于试验设备能力及测量手段的限制，地面模拟难以获得大量可靠有效的试验数据，以膨胀管为例，通常情况其试验时间仅为微秒量级，传统的测量手段很难在如此短的时间内响应。针对超高速流动热化学非平衡效应的数值方法近年来发展迅速，通过数值模拟结合有效试验数据的一体化研究技术在超高速流动研究领域发挥着越来越巨大的作用[10-11]。
超高速流动中存在着复杂的流体物理与热力学过程，强激波的压缩(飞行器头部驻点区域)和剧烈膨胀(头部以下流场)，使得超高速流动成为一种处于非平衡热化学状态介质的非定常流动，高温气体宏观热力学性质的确定非常复杂和困难[12]。****们在理论分析中提出了一系列物理模型来模化上述复杂的物理、化学过程[13-14]。然而，这些模型都不同程度地带有各自的局限性、经验性及不确定性，同一物理问题条件下，不同物理模型计算结果差异显著，各自的物理模型只是依赖于各自的试验数据[15]。因此，对于物理模型的适用性研究工作在不断进行，许多模型也在发展中得到改进[16-17]。
超高速流动数值模拟时，选取不同的热化学反应模型，试验气流特性及试验模型的气体动力学特性都会受到影响。本文研究内容主要是对超高速来流条件下的典型试验模型，比较分析不同热化学反应模型对流场结构及模型壁面参数的影响机制，进一步评估各热化学反应模型的适用范围，从而保证超高速流动模拟的准确性及可靠性。同时，也为实验室即将开展的膨胀管超高速流动显示工作提供计算流体力学(CFD)数据支持。
1 JF-16膨胀管机理及气流参数 JF-16膨胀管主要有3个部分组成，分别是爆轰段(detonation tube)、激波管(shock tube)及加速段(acceleration tube)，其运行过程如图 1所示。当爆轰波冲破爆轰驱动段/激波管之间的主膜后，驱动气体压缩激波管中的试验气体，在激波管中形成一道激波，称为主激波(PSW)，同时中心膨胀波向驱动段的高压气体中传播，主激波和膨胀波之间由驱动/试验(2/3区)气体界面(PCS)分开。当主激波运动到第2道膜并将其冲开后，在加速段形成第2道入射激波(SSW)，同时形成中心膨胀波向上游的2区试验气体中传播，由于2区气体处于超声速状态，因此该膨胀波在实验室坐标系下是向下游运动，并形成试验/加速(5/6区)气体界面(SCS)。2区试验气体经该非定常膨胀波后，进一步加速并降低静温到达最终试验状态，即5区高焓、高速的试验气流。
 |
| 图 1 JF-16波动过程图 Fig. 1 Sketch of wave process for JF-16 |
| 图选项 |
JF-16膨胀管可以产生超过10 km/s的超高速气流，其典型的试验气流和加速气流参数如表 1所示，第3节中根据问题的研究特性分别选取二者作为计算的来流条件。
表 1 试验气流和加速气流参数 Table 1 Flow parameters of test and acceleration gas
| 气流参数 | u/(km·s-1) | T/K | Ma | p/Pa | ρ/(kg·m-3) | |
| 试验气流 | 8.0 | 2 390 | 7.6 | 9 524.550 | 0.011 6 | |
| 加速气流 | 8.8 | 9 335 | 3.2 | 9 220.575 | 0.001 7 | |
| 注：u—x方向速度；T—温度；Ma—马赫数；p—压力；ρ—密度。 | ||||||
表选项
2 数值方法及热化学反应模型 高超声速飞行时，波后气体受到激波压缩同时升温，高温空气会因飞行速度不同发生各种化学反应，使得空气特性偏离量热完全气体或热完全气体[18]。随着飞行器轨道和速度的逐步提升，驻点区域空气会经历量热完全气体、热完全气体到高温气体的变化，热化学状态也由常温无反应气体过渡到单温度(热平衡过程)描述的化学平衡态、化学非平衡态到最终的多温度(热非平衡过程)热化学非平衡态；与此同时，与化学反应相关的气体组分也由2组分(N2、O2)发展到5组分(N2、O2、NO、N、O)甚至11组分(N2、O2、NO、N、O、N2+、O2+、NO+、N+、O+、e-)，需要根据流动的具体特征选择最适合的热化学反应模型[19]。
本文主要关注超高速流动模拟时不同热化学反应模型对流场特性的影响机制，对数值方法进行了简化，忽略黏性和热传导，采用了多组分Euler方程，其二维的控制方程如下：
 | (1) |
其中：
 | (2) |
 | (3) |
 | (4) |
 | (5) |
式中：U为未知变量；F和G分别为x、y方向的对流通量；Sc为化学反应源项；
本文数值模拟时，分别采用了3种热化学反应模型。Model 1和2为热平衡化学非平衡模型，热力学参数则由温度相关的多项式直接给出，与时间无关，化学反应模型选择Dunn和Kang[20]的高温空气化学反应模型，其中Model 1为5组分17基元(1~17) 反应模型，Model 2为11组分31基元(1~31) 反应模型，如表 2所示。Model 3为热化学非平衡模型，热力学参数随时间变化，热力学模型采用Park的双温度模型(T，Tv)，流动过程及热化学反应过程由平动温度T和振动温度Tv来描述，化学反应模型为Dunn和Kang[20]的5组分17基元反应模型。对流项的离散应用了修正后的频散可控耗散格式(DCD)[21]。关于数值方法的可行性及可靠性论证在文献[9, 21]中有详细阐述，该方法已成功应用于膨胀管超高速流动的数值模拟。
表 2 Dunn & Kang高温空气化学反应模型[20] Table 2 Dunn & Kang chemical reaction model for air at high temperature[20]
| 基元数 | 反应式 |
| 1 | O2+N=2O+N |
| 2 | O2+NO=2O+NO |
| 3 | N2+O=2N+O |
| 4 | N2+NO=2N+NO |
| 5 | N2+O2=2N+O2 |
| 6 | NO+O2=N++O+O2 |
| 7 | NO+N2=N++O+N2 |
| 8 | O+NO=N++O2 |
| 9 | O+N2=N++NO |
| 10 | N+N2=2N+N |
| 11 | O2+O=2O+O |
| 12 | O2+O2=2O+O2 |
| 13 | O2+N2=2O+N2 |
| 14 | N2+N2=2N+N2 |
| 15 | NO+O=N++2O |
| 16 | NO+N=O2+N |
| 17 | NO+NO=N++O+NO |
| 18 | O+N=NO++e- |
| 19 | O+e-=O++2e- |
| 20 | N+e-=N++2e- |
| 21 | O+O=O2++e- |
| 22 | O+O2+=O2+O+ |
| 23 | N2+N+=N++N2+ |
| 24 | N+N=N++2+e- |
| 25 | O+NO+=NO+O+ |
| 26 | N2+O+=O+N2+ |
| 27 | N+NO+=NO+N+ |
| 28 | O2+NO+=NO+O2+ |
| 29 | O+NO+=O2+N+ |
| 30 | O2+N2=NO+NO++e- |
| 31 | NO+N2=NO++N2+e- |
表选项
3 热化学反应模型对比研究 3.1 Model 1和2(试验气流条件) 采用Model 1、2的2种热化学反应模型，对处于试验气流条件下的双楔模型进行数值模拟对比研究，为后续进行双楔模型超高速流动显示试验提供CFD数据支持。双楔模型选取了2组外形，前楔角均为15°，后楔角分别为35°和55°，楔角大小的变化可以影响斜激波的强度及激波交互处的反射规律。
图 2为15°~35°双楔流场压力分布，p∞为来流压力。通过比较Model 1和2的结果可以看出，二者的流场结构近似相同，斜激波角度及流场压力分布保持一致，初步判断试验气流条件下，Model 1和2对15°~35°双楔模型计算结果的影响较小。
 |
| 图 2 15°~35°双楔压力分布 Fig. 2 Pressure distribution of 15°~35° double-wedge |
| 图选项 |
提取压力、密度和温度3个主要流场参数沿模型壁面分布情况，如图 3所示。可以直观地看出，3个参数沿壁面分布曲线完全重合，证明试验气流条件下，Model 1和2对计算结果的影响可以忽略。通常情况认为空气发生电离的温度临界点为9 000 K左右，从壁面温度曲线可以看出，第2道斜激波波后温度瞬间峰值接近8 400 K，此时空气的电离程度微弱，Model 2虽然考虑了11组分31个基元反应，但电离反应尚未明显激发，与Model 1对比没有明显差别。
 |
| pw—壁面压力；ρw—壁面密度；ρ∞—来流密度；Tw—壁面温度；T∞—来流温度。 图 3 15°~35°双楔流动参数沿壁面分布 Fig. 3 Flow parameter distribution along wall of 15°~35° double-wedge |
| 图选项 |
提取流场中电子e-的摩尔分数分布图，如图 4所示，并将全流场中粒子的摩尔分数峰值列成表 3。图 4及表 3数据表明，流场中带电粒子摩尔分数量级为小量，再次验证了之前的分析。综上所述，在试验气流条件下，流场中电离反应很弱，带电粒子比例很小，Model 1和2对数值模拟结果的影响可以忽略。
 |
| 图 4 15°~35°双楔电子摩尔分数分布 Fig. 4 Electronic mole fraction distribution of 15°~35° double-wedge |
| 图选项 |
表 3 15°~35°双楔粒子摩尔分数峰值 Table 3 Max mole fraction of species for 15°~35° double-wedge
| 粒子 | 摩尔分数峰值 |
| N | 0.34 |
| O | 0.33 |
| NO | 0.04 |
| N+ | 6.31×10-5 |
| O+ | 5.08×10-4 |
| NO+ | 4.74×10-4 |
| e- | 1.06×10-3 |
表选项
采用相同的分析方法，计算得到15°~55°双楔模型流场压力分布如图 5所示，Model 1和2的计算结果仍然类似。第2楔角由35°增大至55°，激波交汇处的反射规律发生了变化，第2楔角为35°时，两道激波在拐角处是规则反射，随着楔角增大至55°，壁面的偏转角过大，气流经过第2道斜激波后无法立刻满足壁面条件，而是先经过一段过渡区域再转向平行于壁面，激波交汇处出现了一道新的强间断，即马赫杆，此时称之为马赫反射。
 |
| 图 5 15°~55°双楔压力分布 Fig. 5 Pressure distribution of 15°~55° double-wedge |
| 图选项 |
提取出3个流场参数沿壁面分布情况，如图 6所示，3个参数沿壁面分布曲线仍保持完全重合，Model 1和2对计算结果的影响仍然可以忽略。进一步分析温度曲线，由于第2楔角增大至55°，第2道激波强度增大，波后温度瞬间上升至峰值，接近12 000 K，随即迅速回落至8 000 K左右，这反映的是高温真实气体效应对波后气体降温的过程，波后气体振动激发、离解及电离等过程会吸收能量，有效降低波后温度。虽然波后瞬时温度峰值已达到了空气电离的临界温度值，但这一高温状态并没有持续保持，因此，波后气体的电离程度仍然较弱，Model 1和2的差异性并未得到体现。
 |
| 图 6 15°~55°双楔流动参数沿壁面分布情况 Fig. 6 Flow parameter distribution along wall of 15°~55° double-wedge |
| 图选项 |
电子e-摩尔分数分布及粒子摩尔分数峰值如图 7和表 4所示。从图 7及表 4进一步分析可以验证之前的判断，带电粒子摩尔分数量级仍旧是小量，电离反应仍未被激发至可观状态。试验气流条件下，Model 1和2对15°~55°双楔流场的计算结果仍保持一致。
 |
| 图 7 15°~55°双楔电子摩尔分数分布 Fig. 7 Electronic mole fraction distribution of 15°~55° double-wedge |
| 图选项 |
表 4 15°~55°双楔粒子摩尔分数峰值 Table 4 Max mole fraction of species for 15°~55° double-wedg
| 粒子 | 摩尔分数峰值 |
| N | 0.69 |
| O | 0.31 |
| NO | 0.03 |
| N+ | 1.73×10-2 |
| O+ | 6.37×10-2 |
| NO+ | 6.29×10-3 |
| e- | 7.66×10-2 |
表选项
通过上述2组双楔外形的计算结果，可以得到以下结论：在试验气流条件下，2种热化学反应模型计算结果保持一致，5组分11基元反应的热化学反应模型Model 1即可满足计算要求，而不必采用组分更多更复杂的Model 2，可以简化热化学反应模型以提高计算效率。
3.2 Mode l和2(加速气流条件) 为了进一步研究Model 1和2对数值模拟结果的影响，选择另一种比较极端的来流情况，即加速气流条件，加速气流拥有极高的静温，模型选择15°~35°双楔外形。计算得到压力分布如图 8所示，可以明显看出，Model 1和2的流场结构不同，斜激波的角度有明显差异。2组流场的激波结构对比如图 8(c)所示，第1道斜激波角度大小不同，第2道斜激波的角度差异更为明显，且Model 1激波角度相对较大。
 |
| 图 8 15°~35°双楔压力分布(加速气流条件) Fig. 8 Pressure distribution of 15°~35° double-wedge (acceleration gas condition) |
| 图选项 |
3个气流参数沿壁面分布如图 9所示，第2道斜激波后壁面参数分布差异性更明显。从密度分布曲线来看，Model 2波后密度相对较高，由质量守恒定律，可以判断其激波角度相对较小，验证了之前的分析。从温度分布曲线来看，Model 2波后温度相对较低，但已经超过了10 000 K，此时电离反应得到了有效激发，且电离过程会吸收波后能量，导致波后温度相比Model 1有了更大程度的回落。
 |
| 图 9 15°~35°双楔流动参数沿壁面分布(加速气流条件) Fig. 9 Flow parameter distribution along wall of 15°~35° double-wedge (acceleration gas condition) |
| 图选项 |
电子e-摩尔分数分布及各组分峰值如图 10及表 5所示，从中可以看出带电粒子的摩尔分数量级已相当可观，不可忽略。这说明空气的电离程度已经较强，可以明显影响到气流状态及流场结构。
 |
| 图 10 15°~35°双楔电子摩尔分数分布(加速气流条件) Fig. 10 Electronic mole fraction distribution of 15°~35° double-wedge (acceleration gas condition) |
| 图选项 |
表 5 15°~35°双楔粒子摩尔分数峰值(加速气流条件) Table 5 Max mole fraction of species for 15°~35° double-wedge (acceleration gas condition)
| 粒子 | 摩尔分数峰值 |
| N | 0.77 |
| O | 0.21 |
| NO | 6.88×10-5 |
| N+ | 0.11 |
| O+ | 0.05 |
| NO+ | 3.55×10-4 |
| e- | 0.13 |
表选项
通过以上对比研究可以得到以下结论，在加速气流下，Model 1和2对计算结果影响的差异性已经有明显体现，Model 1只考虑了5组分的化学反应，忽略了电离反应，对加速气流条件已不再适用，必须考虑11组分带电离反应的Model 2。因此，对于不同的来流条件，必须根据气流特性选择合适的热化学反应模型。
3.3 Model 1、2和3(试验气流条件) 3.1节和3.2节选取的皆为热平衡模型，本节采用热非平衡模型Model 3进行数值模拟对比研究，选取试验气流条件及15°~35°的双楔外形，并将计算结果与Model 1、2进行对比分析，压力分布如图 11所示。从图中可以看出，Model 3与1、2的流场结构略有不同，第1道斜激波倾角有一定差异，第2道斜激波倾角几乎一致。
 |
| 图 11 15°~35°双楔压力分布(Model 1、2、3) Fig. 11 Pressure distribution of 15°~35° double-wedge (Model 1、2、3) |
| 图选项 |
3个主要气流参数沿壁面分布如图 12所示。观察温度曲线，Model 3波后温度相对较高，且随着时间发展逐步与Model 1、2到达一致。对于Model 3，激波与气体相互作用的前期，波后气体平动温度和振动温度急剧上升且平动温度的上升速率高于振动温度，气体处于热非平衡状态；随后由于粒子间的高速碰撞作用，使得粒子的平动能量和转动能量逐步转移给振动能量，平动温度随时间发展逐步降低并与振动温度达到一致，此时为热平衡状态。而Model 1和2，振动激发不需要时间，波后气体平动温度直接降低至平衡状态温度值，相较于Model 3波后平动温度由峰值随时间逐步降低，Model 3的波后温度前期较高。
 |
| 图 12 15°~35°双楔流动参数沿壁面分布(Model 1、2、3) Fig. 12 Flow parameter distribution along wall of 15°~35° double-wedge (Model 1、2、3) |
| 图选项 |
从尺度效应来看，大尺度模型流动时间尺度相对较长，波后气体热非平衡状态最终随时间发展逐步达到热平衡状态。但是对于小尺度模型，流动特征时间与能量弛豫时间可比，热非平衡效应明显，此时考虑了热非平衡效应的Model 3更为适合。
4 结论 本文简单介绍了爆轰驱动膨胀管JF-16的机理及其获得的超高速气流参数，并主要对比研究超高速流动模拟不同热化学反应模型对流场结构及模型壁面参数的影响机制，得到了以下结论：
1) 膨胀管中速度8.0 km/s的试验气流条件下，流场中电离反应程度微弱，带电粒子比例可以忽略，考虑5组分17基元反应的Model 1已经满足了模拟要求，而不必考虑更多的组分和化学反应。
2) 膨胀管中速度8.8 km/s的加速气流条件下，由于波后温度极高，流场中电离反应得到明显激发，带电粒子比例不可忽略，此时考虑5组分的Model 1难以满足数值模拟要求，必须采用考虑电离反应的11组分的Model 2来保证模拟的准确性。
3) 对高焓流动、小尺度模型而言，流动特征时间与能量弛豫时间可比，此时流动中热非平衡效应显著，则热化学非平衡模型Model 3更符合模拟要求。
参考文献
| [1] | ANDERSON J D. Hypersonic and high temperature gas dynamics[M].New York: McGraw-Hill Book Company, 2006: 449-461. |
| [2] | NEELY A J, MORGAN R G. The superorbital expansion tube concept, experiment and analysis[J].Aeronautical Journal, 1994, 98(973): 97–105.DOI:10.1017/S0001924000050107 |
| [3] | HORNUNG H G.Experimental hypervelocity flow simulation, needs, achievements, and limitations[C]//Proceedings of 1st Pacific International Conference on Aero-space Science and Technology.Tainan:National Cheng-Kung University, 1994:1-10. |
| [4] | RESLER E L, BLOXSOM D E. Very high Mach number flows by unsteady flow principles[M].New York: Cornell University Graduate School of Aeronautical Engineering, 1952: 2-7. |
| [5] | TRIMPI R L.A preliminary theoretical study of the expansion tube:A new device for producing high-enthalpy short-duration hypersonic gas flows:NASA-TR-R133[R].Washington, D.C.:NASA, 1962.https://link.springer.com/chapter/10.1007/978-3-319-23745-9_9/fulltext.html |
| [6] | PAULL A, STALKER R J, STRINGER I.Experiments on an expansion tube with a free piston driver[C]//Proceedings of 15th Aerodynamic Testing Conference.Reston:AIAA, 1988:173-178. |
| [7] | GILDFIND D E, MORGAN R G, MCGILVRAY M, et al. Simulation of high Mach number scramjet flow conditions using the X2 expansion tube:AIAA-2012-5954[R]. Reston:AIAA, 2012. |
| [8] | DUFRENE A, MACLEAN M, PARKER R, et al.Experimental characterization of the LENS expansion tunnel facility including blunt body surface heating:AIAA-2011-262[R]. Reston:AIAA, 2011. |
| [9] | 周凯, 汪球, 胡宗民, 等. 爆轰驱动膨胀管性能研究[J].航空学报, 2016, 37(3): 810–816. ZHOU K, WANG Q, HU Z M, et al. Performance study of a detonation-driven expansion tube[J].Acta Aeronautica et Astronautica Sinica, 2016, 37(3): 810–816.(in Chinese) |
| [10] | HOLLIS B R, PRABHU D K. Assessment of laminar, convective aeroheating prediction uncertainties for Mars entry vehicles[J].Journal of Spacecraft and Rockets, 2013, 50(1): 56–68.DOI:10.2514/1.A32257 |
| [11] | HOLDEN M S.Development of experimental facilities coupled with CFD to research key aerothermal phenomena in hypervelocity flight[C]//Proceedings of AIAA Aerospace Planes Meeting.Reston:AIAA, 2011:243-253. |
| [12] | SARMA G S R. Physico-chemical modelling in hypersonic flow simulation[J].Progress in Aerospace Sciences, 2000, 36(3): 281–349. |
| [13] | GNOFFO P A, GUPTA R N, SHINN J L.Conservation equations and physical models for hypersonic air flows in thermal and chemical nonequilibrium:NASA-TP-2867[R].Washington, D.C.:NASA, 1989. |
| [14] | LOSEV S A, MAKAROV V N, NIKOLSKY V S.Thermochemical nonequilibrium kinetic models in strong shock waves on air:AIAA-1994-1990[R].Reston:AIAA, 1994. |
| [15] | BUCK M L, BENSON B R, SIERON T R, et al. Aerodynamic and performance analyses of a superorbital re-entry vehicle[J].Dynamics of Manned Lifting Planetary Entry, 1963, 15(2): 376–407. |
| [16] | PARK C. Assessment of a two-temperature kinetic model for dissociating and weakly ionizing nitrogen[J].Journal of Thermophysics and Heat Transfer, 1988, 2(1): 8–16.DOI:10.2514/3.55 |
| [17] | PARK C.The limits of two-temperature model:AIAA-2010-911[R].Reston:AIAA, 2010. |
| [18] | GUPTA R N, YOS J M, THOMPSON R A.A review of reaction rates and thermodynamic and transport properties for the 11-species air model for chemical and thermal nonequilibrium calculations to 30 000 K:NASA-TR-1232[R].Washington, D.C.:NASA, 1989. |
| [19] | 柳军. 热化学非平衡流及其辐射现象的实验和数值计算研究[D]. 北京: 国防科学技术大学, 2004. LIU J.Experimental and numerical research on thermo-chemical nonequilibrium flow with radiation phenomenon[D].Beijing:National University of Defense Technology, 2004(in Chinese). |
| [20] | DUNN M G, KANG S.Theoretical and experimental studies of reentry plasmas:NASA-CR-2232[R].Washington, D.C.:NASA, 1973. |
| [21] | HU Z M, WANG C, JIANG Z L, et al. On the numerical technique for the simulation of hypervelocity test flows[J].Computer and Fluids, 2015, 106: 12–18.DOI:10.1016/j.compfluid.2014.09.039 |
