风洞试验是认识高超声速流动机理和预测飞行器性能不可或缺的手段。激波风洞是支撑高马赫数飞行地面气动试验研究的主力设备之一。它主要利用运动激波对流体的气动加速和气动压缩作用产生试验所需的高焓气流，因此，其工作效能极大地依赖于喷管上游激波管内以运动激波、稀疏波等为特征的强非定常流动。对激波风洞整体流动的充分理解和认知是此类风洞有效服务于气动试验的重要前提。

为了达到模拟高马赫数飞行环境所需的总温、总压条件，激波风洞运行时，驱动段以及激波反射后的驻室区域的气体往往处于高温、高压的状态。在极端的温度和压强条件下，气体的物理化学属性与状态显著偏离理想气体的描述，从而在流动中产生所谓的真实气体效应。这里的真实气体效应实际包括2类。一是高温真实气体效应。由于气体温度的提升，气体分子原本大量处于基态的振动能和电子能被显著激发，由此产生分子多种内能模态间（转动能、平动能、振动能和电子能）有限速率的相互转化；同时，热运动加剧使得分子内部化学键断裂乃至电子脱离原子核束缚，气体发生离解和电离。这种高温气体中发生的热弛豫和化学反应，导致气体宏观可感内能的改变以及气体物理化学属性的变化^[1]。二是高压真实气体效应。在气体温度维持低值时，高压意味着分子数密度增加；当气体状态逼近液化临界或超临界态，气体分子体积及分子间作用力逐渐变得不可忽略，经典气动理论使用的理想气体状态假设将不再适用，气体的热容、内能、熵等热力学参数同时发生一定程度的偏移^[2]。高状态运行的激波风洞流动中常常兼具上述2类真实气体效应，基于经典气体动力学理论的计算已不足以准确评估风洞的运行状态。

对激波风洞流动的测试和研究是高超声速试验技术的重要组成部分，贯穿于风洞的设计、运行和对试验结果的分析环节，如：Jacobs^[3]发展了准一维数值模拟程序L1D，用于脉冲风洞内流场模拟，并成功服务于昆士兰大学T4风洞运行状态的预测评估；Boyce等^[4]对T4激波风洞不同工况下的试验气流参数以及有效试验时间进行了分析；Hannemann^[5]和Johnston等^[6]以试验结合数值模拟考察了HEG高焓激波风洞复现X-38等飞行器再入大气条件的流动状态；Holden等^[7]在LENS-I和LENS-II激波风洞中获得了高总温和高总压的试验气流，并评估了真实气体对超高速再入飞行器性能的影响；卢洪波等^[8]在FD-21高焓激波风洞中开展高马赫数吸气推进试验技术探索。近20年来，人们对于高焓脉冲风洞的流动模拟和计算普遍注意将高温真实气体效应纳入考虑。Llobet等^[9]考虑振动非平衡和化学平衡对激波风洞的运行过程进行计算，获得试验气流参数和有效试验时间随初始压比的变化规律。2022年，Lynch等^[10]使用L1D4程序（耦合NASA-CEA^[11]计算热化学平衡）为自由活塞驱动反射型激波风洞设计了总压20 MPa、总温3700 K的工况，复现了来流马赫数为8～9的飞行环境；Luo等^[12]采用准一维数值模拟方法，考虑化学平衡模型，对JFX爆轰驱动激波风洞内的流动进行精确快速的模拟和研究，掌握了激波风洞运行过程中的波系运动情况。但学者们对高焓脉冲风洞中高压真实气体效应的研究则较少。MacLean等^[13-14]在对LENS-I激波风洞和LENS-XX膨胀管风洞的分析过程中使用了考虑体积修正的状态方程和热平衡气体，发现考虑体积修正后的自由流单位雷诺数明显增大；Candler^[15]在使用克劳修斯状态方程和热化学非平衡的基础上，对AEDC激波风洞的喷管流动进行计算分析，并指出克劳修斯状态方程在低温高密度状态下无法完全体现高压气体效应。由于考虑完全的立方型状态方程后，流场温度、压强和密度的非线性关系显著增加了数理建模与流动解算的复杂性，Sirignano^[16]提出对成熟的真实气体状态方程进行线性化处理，并用线性化后的状态方程分析了等熵流等基本流动过程，但是这些简化处理仅在有限的条件范围内适用，当流场状态跨度较大或者测试条件范围较宽时，简化方程偏差显著。因此，仍有必要采用完备的状态方程，系统地分析激波风洞流动中高压真实气体的影响效应，综合评估考虑高压气体效应前后激波风洞波系结构和试验气流参数的差异，揭示高压真实气体效应的作用机制。

本文中，以FD-14A激波风洞^[17]为基本构型，基于自主发展的准一维热化学非平衡数值模拟平台，引入完全的立方型气体状态方程描述高压真实气体状态，对该型风洞在试验气流总温2000～6000 K、总压20～200 MPa典型工况下的流动过程进行计算和分析；考虑到激波风洞运行过程中高温真实气体效应与高压真实气体效应是同时存在的，在考虑热化学非平衡基础上对比采用真实气体状态方程和理想气体状态方程的计算结果，着重关注高压真实气体效应导致的激波风洞运行全场流动时空特性和试验气流参数的差异以及产生差异的物理机制。

1.   风洞构型与数值模拟方法

图1为FD-14A激波风洞^[17]内流道的基本构型。该风洞激波管部分管道的内径均为0.15 m；当前研究设置驱动段长9 m，被驱动段长18 m。为考虑反射端的质量通量，设置喷管喉道的直径为23.9 mm。确定喉道和收缩段后，喷管其余部分不会对实际激波管流动产生显著影响，为了保留激波风洞的全部结构特征，计算过程中保留了整个喷管。

图  1  FD-14A激波风洞^[17]内流道结构示意图

Figure  1.  Configuration of the inner flow channel of the FD-14A shock tunnel^[17]

下载: 全尺寸图片 幻灯片

该风洞驱动段气体为高压氢气（H₂），被驱动段气体为空气或氮气（N₂）。本研究中，为重点考察高压真实气体效应，被驱动段统一采用氮气。

1.1   准一维流动控制方程

对于以大长径比管道流动为特征的高焓脉冲风洞设备，准一维流动数值模拟一直是一种行之有效且快捷的计算分析手段^{[3, 12, 14]}。通过合理的假设，准一维数值模拟可将管道壁面的摩擦与传热、高温气体热化学非平衡等各种物理化学效应同时纳入考量。本文中所使用的考虑双温度振动非平衡、化学反应和输运效应的准一维流动控制方程如下：

式中：U为守恒项；F为对流项；P为变截面管道壁面的流向作用力；A为流道面积函数，描述管道截面空间变化；ρ、u、p、e、e_v、Y_i分别为气体的密度、速度、压强、质量总能量、振动能和组分$i$的质量分数。非平衡气体的质量总能量有：

式中：T为平动温度，T_v为振动温度；h_i⁰和$c_{V,i}^{\rm{rt}}$分别为组分i的生成焓和冻结比定容热容（温度无关）；振动能e_v = Σe_v,i为振动温度的函数。对于组分i，有：

式中：R_u为普适气体常数，w_i为组分i的分子摩尔质量，a_1,i、a_2,i、a_3,i、a_4,i和a_5,i为该组分的热力学系数。Δe_p为高压真实气体效应的能量修正项，与真实气体状态方程相关，后续章节将给出推导形式。

方程(1)中J_in和J_bd分别为内流和壁面输运源项。其中，J_bd用于描述管道壁面摩擦与传热对平均流动的影响。本文中主要采用Groth等^[18]的壁面效应模型，组合Jacobs的参考温度形式，并进行了适配二温度（平动和振动温度）非平衡模型的改良。S为热化学非平衡源项，包括化学反应组分生成源项和振动能弛豫源项。以上3个源项的具体形式为：

式中：λ_i为由组分扩散导致的第i组分含量变化，τ为由黏性引起的流体内部应力，q_tr和q_v分别表示平动能以及振动能的热传导，M表示壁面黏性导致的动量传递，Q_tr和Q_v分别表示壁面传热导致的平动能和振动能传递，${\dot \omega _i}$表示由化学反应导致的第i组分生成，s_v表示由振动弛豫导致的能量转化。对于计算过程中需使用的气体输运系数（主要是黏性系数和热传导系数），单一组分的输运系数以Sutherland公式计算，混合气体输运系数由单一组分系数按Wilke^[19]的方法混合得到。

壁面动量损失源项M和壁面能量损失源项Q的表达式分别为：

式中：D为管道水力学直径；f为Darcy-Weisbach摩擦因数^[18]；c_p为流体的比定压热容；Pr为普朗特数；T_w为当地壁温；T_aw对平动能为绝热壁温，取值为经过可压缩修正之后的气流平均温度，对振动能为当地平均振动温度。

本文中振动弛豫采用Park^[20]的双温度模型，特征弛豫时间由Millikan等^[21]的半经验表达式计算得到，并参考Park^[20]的方法对极高温度下的弛豫时间进行修正。空气在高温下的化学反应采用Gupta等^[22]的基元反应动力学机理（本文中仅涉及氢气和氮气，因而仅取其中氮气离解反应）。具体模型和数据可参考文献[22-23]。

为了在统一模型架构下求解激波风洞的高温高压真实气体流动，在考虑高温气体效应的基础上，采用以下立方型气体状态方程：

封闭方程组(1)。式(6)中：v=1/ρ为气体比容，R=R_u/w为气体常数，w为分子摩尔质量；a、b、k₁和k₂为模化分子间相互作用以及分子体积的参数或函数。当a和b取值为零时，方程(6)回归为理想气体状态方程。

1.2   数值计算方法

准一维数值模拟采用基于非结构网格的有限体积法离散并解算方程(1)，考虑到源项S中化学反应的加入会使得方程具有一定的刚性，因此，采用时间分裂步法将流动与热化学源项解耦求解，首先计算热化学冻结条件下的流动过程，再在固定网格内求解定容绝热热化学弛豫的刚性常微分方程组。纯流动过程采用二阶时空精度的MUSCL-Hancock格式求解，数值通量通过HLLC^[24]格式计算得到，时间步长由CFL (Courant-Friedrichs-Lewy)条件控制；热化学非平衡源项部分采用成熟的刚性常微分方程求解器DVODE进行计算。

1.3   数值方法验证

本研究基于自主开发的准一维热化学非平衡数值模拟程序，所采用的基本模型和方法已在1.1节和1.2节中描述。该程序的基本算法和输运模型已在多个研究中获得充分验证^[25-27]，并在实际脉冲风洞结构和工况设计中获得应用^[28]。

首先，采用激波管试验对上述壁面摩擦和传热模型进行校验。将含壁面输运模型的准一维数值程序计算得到的激波管中入射激波马赫数Ma_s沿程衰减情况与试验结果进行对比，如图2所示。

图  2  准一维数值模拟的激波沿程衰减情况与试验结果的对比

Figure  2.  Comparison of shock attenuation between quasi-one-dimensional numerical simulations and tests

下载: 全尺寸图片 幻灯片

图2显示，准一维数值模拟可很好地刻画广泛条件下运动激波由于壁面摩擦和传热而发生衰减的现象。

接下来，对引入立方型状态方程后的程序进行补充验证。选取表1所示的典型工况，计算激波风洞驻室压强随时间演化情况，并与试验测量压强数据进行对比，如图3所示。

表  1  激波风洞运行典型工况

Table  1.  Typical operating conditions of the shock tunnel

管段	压强/Pa	温度/K	初始组分	长度/m
驱动段	5×107	300	氢气	9
被驱动段	1×105		空气	18
喷管	1		空气	7.5

下载: 导出CSV 

| 显示表格

图  3  不同状态方程下激波风洞驻室压强数值模拟结果与试验结果的对比

Figure  3.  Comparison of numerically-simulated pressures in the stagnation zone of a shock tunnel based on different equations of state with the experimental results

下载: 全尺寸图片 幻灯片

由图3可知，数值模拟计算得到激波到达时刻以及驻室区压强变化情况与试验数据基本吻合。考虑立方型气体状态方程后的数值模拟更准确地捕捉到了稀疏波系的运动情况，计算结果与试验结果的吻合度显著提升。以上结果表明，本文中采用的准一维数值模拟方法和程序可以准确预测激波风洞的整体运行状态。

2.   结果与讨论

2.1   立方型气体状态方程和热力学参数

状态方程建立气体压强、温度和密度（或比容、数密度等）之间的关系。理想气体（ideal gas, IG）状态方程将气体分子视为质点，而不考虑气体分子体积以及分子间作用力，因此在低温高压条件下，它无法准确描述气体状态。van der Waals^[2]建立了第一个适用于真实气体的立方型状态方程，即van der Waals（VDW）方程，该方程将分子体积和分子间作用力的影响分别以斥力项p_r和引力项p_a的形式体现在压强表达式中^[29]：

式中：参数b表征分子体积，参数a度量分子间吸引力；考虑到分子体积会对分子间作用力产生影响，因此引力项中还包含关于体积（比容v）的函数g(v)，VDW方程中取g(v)=v²。参数a、b与气体临界参数（临界温度T_cr、临界压强p_cr）相关，其中，a₀和b₀为状态方程设定的常系数：

后续研究者在VDW方程基础上发展出数百种立方型状态方程^[30]，其中，Soave-Redlich-Kwong (SRK)方程^[31]和Peng-Robinson (PR)方程^[29]是2种代表性的立方型气体状态方程。他们通过采用形式更复杂的g(v)并引入温度修正系数α来更好地描述气体在高压条件下的状态。

SRK方程^[31]的具体形式为：

PR方程^[29]的具体形式为：

温度修正系数α有：

式中：f_ω为分子偏心因子ω的函数。

SRK方程^[31]和PR方程^[29]适用于单一组分，对于混合气体，由于所含各组分的参数a、b不同，需要在考虑各组分摩尔分数的基础上进行平均化处理。常见平均化处理方式有以下2种。一种是基于组分二元交互作用的经典混合规则^[32]：

式中：x_i、a_i和b_i分别为第i组分的摩尔分数以及对应的参数，δ_ij为二元作用参数（这里参考Poling等^[33]提出的最简单的处理方法，即取δ_ij=0）。Peng等^[29]在提出PR方程时对于混合气体参数采用了上述规则^[32]。另一种混合气体平均化处理方法是Kay^[34]提出的准临界参数混合规则，该规则首先对混合气体的临界温度和临界压强进行平均处理，再由式(8)计算得到相应的平均参数a和b。混合气体的临界温度和临界压强的计算公式分别为：

式中：T_cr,i和p_cr,i分别为i组分的临界温度和临界压强。

将2种混合规则计算得到的70%氢气和30%氮气混合气体的可压缩性因子（z=pv/RT）随气体压强p的变化与NIST试验数据^[35]进行对比，如图4所示。可以看到，使用经典混合规则计算方法获得的可压缩性因子大幅偏离试验值，而使用准临界参数混合规则可以较好地描述混合气体状态。因此，后续计算将使用式(14)进行平均处理。

图  4  不同混合规则下300 K氢-氮混合气体可压缩性随压强的变化

Figure  4.  Variations of compressibility factors with pressure for hydrogen-nitrogen mixture at 300 K under different mixing rules

下载: 全尺寸图片 幻灯片

气体的热力学参数（内能、焓、熵和比热等）的具体形式与气体状态方程密切相关，气体状态方程变化后，需要对上述参数进行重新推导。针对立方型方程一般形式(6)，由热力学基本方程推导获得真实气体状态方程与理想气体状态方程内能之差（式(2)中Δe_p项），表达如下：

式中：Δe_0,ref=Δe₀(v=v_ref, T=T_ref)，是与满足理想气体假设的参考体积v_ref和参考温度T_ref有关的常数。获得内能增量后，可进一步计算得到比定容热容c_V、比定压热容c_p以及声速c：

式中：c_V,IG为理想气体状态下气体的比定容热容，γ为比热比。上述方程中保留部分导数以保持方程形式简洁，这些导数均可由状态方程本身或$\alpha$表达式直接求导得出，这里不再展开。

气体状态方程的选取，特别是温度修正系数的具体形式会对气体的热力学参数产生不可忽视的影响，选取合适的状态方程是分析激波风洞时空特性和试验气流参数的前提。依据以上各式计算不同状态方程（VDW方程、SRK方程和PR方程）下氢气的可压缩性因子和声速随压强的变化，并与NIST试验数据^[35]进行对比，如图5～6所示。

图  5  不同状态方程300 K氢气可压缩性因子随压强的变化

Figure  5.  Variation of compressibility factor with pressure for hydrogen at 300 K with different equations of state

下载: 全尺寸图片 幻灯片

图  6  不同状态方程300 K氢气声速随压强的变化

Figure  6.  Variation of sound speed with pressure for hydrogen at 300 K with different equations of state

下载: 全尺寸图片 幻灯片

3种立方型方程均能很好地呈现气体可压缩性和声速随压强升高逐步升高的特征。其中，VDW方程在高压下显著高估可压缩性和声速，而SRK和PR方程两者的结果比较接近，且与NIST试验数据较接近。对比后两者，SRK方程在可压缩性预测上更有优势，而PR方程对声速的预测适用范围更宽。对于压强150 MPa的氢气，PR方程计算得到的声速为2 600 m/s，相比试验值2 480 m/s高约5%；而此时，理想气体方程计算得到的声速为1317 m/s，相比试验数据低约47%。综合考虑，PR方程更适用于分析高压氢气驱动的激波风洞流动问题，后续数值模拟采用PR方程。

2.2   激波风洞流场结构和试验气流参数

运用前述准一维数值模拟方法对激波风洞的全场内流进行计算，通过对比采用立方型气体状态方程和理想气体状态方程的算例，考察高压真实气体效应导致的流场结构和流动参数的差异，考虑到试验时间最大化的需求，本文全部工况中均在理想气体状态方程下通过调整驱动段和被驱动段气体状态，使激波风洞处于缝合接触面运行状态。同一工况、不同状态方程的计算，激波管初始条件相同。

为避免引入多余扰动，计算过程考虑一个理想的瞬态、无损耗的破膜过程。是否考虑非理想破膜对于本文中关注的高压真实气体效应作用规律和机理没有本质影响。

2.2.1   激波风洞全场流动结构

以表2所示的典型工况为例考察激波风洞全场准一维流动的时空结构，该工况以150 MPa氢气驱动110 kPa氮气。在使用理想气体状态方程并考虑管壁沿程损耗情况下，该工况中经过激波压缩后的试验气流总参数为：总压100 MPa、总温6000 K。

表  2  激波风洞运行工况

Table  2.  Operating conditions of shock tunnel

管段	压强/Pa	温度/K	初始组分	长度/m
驱动段	15×107	300	氢气	9.0
被驱动段	11×104		氮气	18.0
喷管	10		氮气	7.5

下载: 导出CSV 

| 显示表格

数值模拟获得理想气体状态方程下激波风洞全场压强、温度x-t云图分别如图7～8所示。由图7～8可以看到：主膜片破裂产生一道右行入射激波、一道接触面和一簇左行中心稀疏波。入射激波到达喷管入口后，喉道前方膜片破裂，气流进入喷管并在喉道发生壅塞，从而反射一道激波向上游运动。喉道前方气流在反射激波作用后成为高温高压的驻室气流。

图  7  激波风洞全场压强x-t云图（理想气体状态方程）

Figure  7.  x-t diagram of pressure contour in the shock tunnel with ideal gas state equation

下载: 全尺寸图片 幻灯片

图  8  激波风洞全场温度x-t云图（理想气体状态方程）

Figure  8.  x-t diagram of temperature contour in the shock tunnel with ideal gas state equation

下载: 全尺寸图片 幻灯片

驻室气流的参数及其稳定性直接决定了喷管出口试验气流的状态及有效试验时间。如图7所示，在缝合接触面运行状态下，对驻室气流产生影响的波系主要为反射稀疏波头和接触面。定义入射激波抵达驻室时间为t₁，反射稀疏波头抵达驻室的时间为t₂，接触面抵达驻室的时间为t₃，则驻室气流参数相对稳定的有效试验时间t_eff可以定义为：

采用与表2中相同的初始条件，使用PR状态方程进行考虑高压真实气体效应的模拟计算，得到的激波风洞全场的压强、温度x-t云图分别如图9～10所示。

图  9  激波风洞全场压强x-t云图（PR状态方程）

Figure  9.  x-t diagram of pressure contour in the shock tunnel with PR state equation

下载: 全尺寸图片 幻灯片

图  10  激波风洞全场温度x-t云图（PR状态方程）

Figure  10.  x-t diagram of temperature contour in the shock tunnel with PR state equation

下载: 全尺寸图片 幻灯片

对比2种状态方程下的整体流场结构。首先，两者波系的时空结构特征相似，并且，入射激波以及接触面的运动轨迹在激波反射前没有明显区别，入射激波均在7 ms左右到达喷管入口并发生反射，说明入射激波强度相当；其次，激波反射之后，反射激波以及接触面的后续运动状态两者存在可见差异，图7中反射激波的运动速度明显慢于图9中的，接触面在与反射激波作用后，速度减小，缓慢向喷管移动，其中，图7中接触面在约23 ms进入喷管，而图9中接触面在25 ms内未能进入喷管；其三，稀疏波运动轨迹存在显著差异，PR方程下稀疏波和反射稀疏波速度显著高于理想气体，图7中稀疏波头在约7.0 ms时刻到达风洞左壁面，图9中稀疏波头到达左壁面时间提前至约3.5 ms，反射稀疏波到达驻室的时间也大幅提前。

为了定量分析特征波系的运动情况，提取入射激波、反射激波、稀疏波头和反射稀疏波头运动速度的沿程分布情况，如图11所示。图11中空心点为IG方程计算结果，实心点为PR气体方程计算结果。根据一维非定常气动理论，稀疏波相对气体以声速传播，左行稀疏波传播速度为u−c，右行稀疏波传播速度为u+c。由于驱动段气体状态均匀，因此左行中心稀疏波头大致以恒定的声速向风洞左端运动。其中，理想气体方程下稀疏波头运动速度为1260 m/s，PR气体方程下稀疏波头运动速度为2370 m/s。中心稀疏波在管壁反射后向右传播，同样，PR气体方程下的反射稀疏波头绝对运动速度大于理想气体方程下的稀疏波速。从图11还可以看到，尽管被驱动段初始流场均匀，但受管壁摩擦以及热传导的影响，入射激波沿程强度有所衰减，激波运动速度逐渐降低，使用PR气体方程计算得到的入射激波初始运动由2631 m/s逐渐减低至2 544 m/s，使用理想气体方程计算得到的入射激波初始运动由2751 m/s逐渐减低至2635 m/s。PR气体方程计算得到的入射激波速度略低于理想气体方程计算得到的入射激波初始速度。而对于反射激波，两者对比趋势相反，PR气体方程的反射激波速度大于理想气体方程的反射激波速度。

图  11  激波管内波系运动速度沿程分布

Figure  11.  Velocities of waves along the shock tube

下载: 全尺寸图片 幻灯片

关于气体的热化学非平衡效应，计算显示，在喉道上游的激波管流动中，流场温度尚不足以导致显著的氮分子离解（尽管驻室区温度可达6000 K，但高压抑制了离解的发生）；同时激波后振动弛豫的尺度远小于管径尺度，因此，整个激波管流动近似处于振动平衡和化学冻结的状态。后续讨论中将不再单独讨论振动温度和离解组分。

2.2.2   驻室状态

激波风洞内波系速度的变化会导致波系时空相干关系的改变，并共同导致驻室状态参数和有效试验时间的变化。图12为不考虑真实气体效应（即，热化学冻结（thermochemically frozen, TF），采用理想气体状态方程（IG EOS））、仅考虑高温气体效应（即，热化学非平衡（thermochemical nonequilibrium，TN，采用理想气体状态方程（IG EOS））和考虑高温、高压气体效应（即，热化学非平衡（thermochemical nonequilibrium，TN，采用PR状态方程（PR EOS））3种情况下激波风洞驻室区（激波管右端，距离喉道0.15 m处）的压强、平动温度和当地气流组分摩尔分数随时间变化情况。

图  12  激波风洞驻室参数随时间变化情况

Figure  12.  Parameters changing with time in the stagnation zone of the shock tunnel

下载: 全尺寸图片 幻灯片

如图12所示，主膜片破膜约7 ms后，入射激波首先抵达右端，当地压强、温度跃升；在经历短暂平台后，反射激波抵达，压强和温度再次跃升；在经历一段波动后，压强和温度信号趋于平稳，其中，压强信号缓慢上升，温度信号则呈缓慢降低趋势，这一区域内参数无明显波动，是激波风洞试验主要利用的区域。

对比图12中3种工况计算得到的驻室参数情况，激波管流动的高温气体效应主要体现为氮气分子振动能激发和小部分氮分子的离解降低了驻室区的温度。该效应对全场流动时空特性影响不大。为了分析高压气体效应的影响，后续不再详述高温气体效应的影响，仅在同等考虑高温气体效应情况下，针对使用理想气体状态方程和PR状态方程的流动差异进行讨论。

理想气体状态方程和PR状态方程下的驻室区压强信号分别在20.6和13.3 ms时刻开始下降。对于温度，理想气体状态方程下的驻室区温度在16.7 ms附近迅速下降；而PR状态方程下的温度曲线无明显趋势改变，仅在13.3 ms附近轻微向下折转。由图12的组分变化可知，理想气体下16.7 ms后氢气组分开始占据驻室区，表明此时接触面抵达；而PR方程下，25 ms内驻室区均以氮气为主，未见明显的氢组分，表明接触面尚未进入观测区间。

结合波系结构和组分演变可知，图12中压强的迅速下降主要对应于反射稀疏波的抵达，而温度的迅速下降主要对应于接触面（冷的驱动气体）的抵达。由于试验气流需要同时维持压强、温度和组分的稳定，有效试验时间由最先的抵达的稀疏波或接触面决定。

提取各个波系到达驻室的时刻如表3所示。受高压真实气体效应影响，入射激波抵达时间由6.8 ms轻微延迟至7.1 ms，而反射稀疏波头到达驻室的时间则由20.6 ms大幅提前至13.5 ms，接触面到达驻室时间由17.2 ms延迟为25.5 ms。这使得限制有效试验时间的因素由接触面变为反射稀疏波，实际有效试验时间由10.4 ms缩短为6.4 ms，降幅达38.5%。

表  3  各个波系到达驻室入口时刻

Table  3.  Arrival time of each wave to stagnation zone

方程	${t_1}$/ms	${t_2}$/ms	${t_3}$/ms	${t_{{\text{eff}}}}$/ms
理想气体状态方程	6.8	20.6	17.2	10.4
PR状态方程	7.1	13.5	25.5	6.4

下载: 导出CSV 

| 显示表格

考察有效试验时间内的驻室气流参数情况。首先，由于PR方程下入射激波强度较弱，因此其有效试验时间内的压强和温度均略低于理想气体计算结果。其次，由于入射激波沿程传播受管壁摩擦和传热的影响强度逐步减弱，波后气流非均匀，反射激波后的气流参数不能维持平台（试验数据同样显示上述特征，如图3所示）。从图12可以看到，有效试验时间内：在理想气体方程下，驻室压强由80 MPa逐步上升到120 MPa，平均驻室压强约为100 MPa，驻室温度由6600 K逐渐降低到6100 K，平均驻室温度约为6350 K；而在PR气体方程下，驻室压强由80 MPa逐渐上升到100 MPa，平均驻室压强约为90 MPa，驻室温度由6300 K逐渐降低到5800 K，平均驻室温度约为6050 K。

2.3   高压真实气体对激波管流动的影响机制

上述数值结果表明，考虑高压真实气体效应后风洞流场结构、试验气流参数和有效试验时间均发生了明显变化。为了揭示高压真实气体对激波管流动的影响机制，接下来将对采用理想气体状态方程和PR状态方程下激波管内稀疏波、激波以及后续接触面的传播过程进行对比分析。

2.3.1   稀疏波传播速度

考虑高压真实气体效应后稀疏波以及反射稀疏波的传播速度明显加快。稀疏波以声速传播，因此稀疏波的传播速度加快实际上是当地声速在高压真实气体下发生增大所导致的。

根据式(19)，以PR方程为状态方程时，气体的声速可写作：

式中：$\sqrt {\gamma RT}$为理想气体状态方程下的声速。与状态方程本身相对应，PR方程下高压真实气体效应对声速的影响可分解为分子体积和分子间作用力两因素，分别由两无量纲因子η_v和η_F表达。当分子体积与分子间作用力趋于0时，η_v和η_F均趋于1，声速回归至理想气体状态方程下的声速。

为了评估上述两因素的影响程度，分别提取η_v和η_F随压强变化情况，如图13所示。可以看到，分子体积项的影响随着压强的增高而显著增强，而分子间作用力的影响随压强增高变化不大，数值维持在1～0.95之间。由此可知，分子体积对空间的占据效应对于声速的影响占据主导地位。

图  13  分子体积和分子间作用力对氢气声速的影响因子随压强的变化

Figure  13.  Influencing factors of molecular volume and intermolecular force on the sound speed of hydrogen changing with pressure

下载: 全尺寸图片 幻灯片

2.3.2   激波管入射波系和驱动能力

除改变稀疏波传播速度外，气体状态方程对激波和稀疏波的强度也存在不同程度的影响，这些影响可改变激波管和激波风洞的工作状态。这里首先从气体动力学基本方程出发分析真实气体对激波管入射波系和驱动能力的影响机制。

激波管特征波系由向低压区传播的入射激波、向高压区传播的稀疏波和中间隔离驱动气体与被驱动气体的接触面构成。三波系将激波管流动分为4个平台区，依次标记为区域①~④，入射激波在端壁发生反射将区域②气流的动能转化为气体内能并叠加到区域②气流既有的内能中，由此产生高温和高压的区域⑤。如图14所示为激波管流动波系和流场分区的示意图。

图  14  激波管流动波系和区域示意图

Figure  14.  Schematic diagram of wave system and flow regions in the shock tube

下载: 全尺寸图片 幻灯片

由于高温热力学参数和PR方程形式的复杂性，无法解析获得高温高压真实气体下的激波关系和稀疏波关系的具体形式。为此，只能由原始方程出发解算激波管的入射波系及其流场参数：(1)入射激波关系由跨间断的质量、动量、能量守恒方程外加PR方程组成，4个方程可建立波后压强p、比容v、温度T和速度u等4个未知数和入射激波强度的关系。PR方程和平衡热力学参数主要通过能量方程中的焓h发生作用。间断关系方程全部为代数方程，因而以牛顿迭代数值求解（由已知的区域①参数出发，获得区域②的参数）方程组：

式中：下标为1的表示波前气体参数，下标为2的表示波后气体参数。

(2)稀疏波关系由小扰动波的连续性关系和等熵关系外加PR方程组成，3个方程可建立压强、比容、温度和速度（p、v、T、u）4个未知数中任意两者的关系。稀疏波关系中的相容关系为常微分方程，因而以R-K方法积分求解（以区域④的条件为初始条件求解区域③的参数）如下联立方程组：

式中：s、c_V分别为气体的比熵、比定容热容，c_V,IG为理想气体状态下气体的比定压热容，c和u为当地声速和当地速度。

区域②和区域③由接触面隔开，因此区域②与区域③的速度（u）和压强（p）应当分别相等。可在p-u图上寻找激波管问题的解点。

在激波管区域④和区域①参数已知的情况下，可分别绘制出入射激波后和稀疏波后气流的p-u图，如图15所示。由图15可知，高压真实气体效应对稀疏波关系产生了较明显的影响，驱动段气体由相同压强膨胀时，稀疏波后相同速度处基于PR状态方程的流场压强更低，并且这种差异会随着稀疏波的增强而扩大。与此同时，高压真实气体效应对入射激波间断关系几乎没有影响，两者的波后p-u曲线几乎重合（图15中灰色实线与红色虚线）。这主要是因为：尽管激波后压强显著增大，但同时发生的温度提升大幅压制了气体密度或分子数密度的增长，使得高压气体效应无法显现。为了更清晰地体现这一效应，图16给出了PR方程所定义的氮气可压缩性因子随温度和压强的变化曲面（温度300～4 000 K，压强0.1～200 MPa），任意热力学过程均对应于该曲面上的一条曲线。图16同时给出了初始温度300 K、初始压强0.1 MPa情况下对应等温、等熵和激波压缩3种路径的参数变化曲线。可以看到：等温压缩路径下温度维持不变，可压缩因子具有最显著的上升趋势；等熵压缩路径下温度和可压缩因子变化趋势均居中；激波压缩路径下温度上升最快，可压缩因子则上升不到0.5%，表明该初始条件下的激波压缩过程几乎不呈现高压真实气体效应。

图  15  不同状态方程下激波管区域②和区域③气流的p-u曲线

Figure  15.  p-u diagrams of the flows in regions ② and ③ of the shock tube for different equations of state

下载: 全尺寸图片 幻灯片

图  16  不同压缩过程高压气体效应的影响程度

Figure  16.  Influences of high pressure gas effect in different compression processes

下载: 全尺寸图片 幻灯片

当激波管中激波和稀疏波后气流在接触面处匹配流场速度和压强时（即图15中激波和稀疏波曲线交点），尽管激波关系几乎不变，但由于稀疏波关系曲线的偏移，考虑高压真实气体效应获得的激波强度仍要低于理想气体情形。然而这里激波强度降低的幅度是有限的，理论计算显示，在表1条件下，激波管产生的初始激波马赫数大致由7.9降至7.6，这与准一维数值模拟结果是一致的。

为呈现广泛条件下驱动段、被驱动段间压强配比与入射激波间的关系，按照以上理论方法对2种状态方程下不同高、低压条件驱动产生的入射激波强度进行计算。图17所示为以压强200和50 MPa氢气驱动压强10～1000 kPa的氮气的入射激波马赫数，其中空心点为理想气体状态方程的计算结果，实心点为PR状态方程的计算结果。由图17可以看到：首先，考虑高压效应的入射激波马赫数要低于理想气体情形，这一规律在全范围内维持不变；其次，高压真实气体效应对入射激波马赫数的影响幅度随驱动段压强和高、低压段压比的增大而增大，但整体变化幅度并不十分显著，在所测试的最极端情形下（200 MPa氢气驱动10 kPa氮气，压比20 000），理想气体状态方程下入射激波马赫数约为11.37，PR方程下激波马赫数约为10.78，两者差异仅为0.55。

图  17  理想和PR状态方程下不同初始压强配置所产生的入射激波马赫数

Figure  17.  Mach numbers of incident shock waves produced under different initial pressure conditions with ideal and PR equations of state

下载: 全尺寸图片 幻灯片

2.3.3   激波管后续波系传播

根据理论计算，考虑高压真实气体效应后，高压驱动段的气体声速由理想气体中1 317 m/s提高至2 596 m/s，这使得以声速传播的入射稀疏波头更快抵达左端壁。为进一步分析反射波系的传播特征，将图15中理想气体状态方程和PR状态方程计算得到稀疏波后声速和激波后声速给出，如图18所示，图18中2条竖直线分别对应图15中高压真实气体和理想气体的解。

图  18  不同状态方程下激波管区域② 和区域③气流的c-u曲线

Figure  18.  c-u diagrams of the flows in regions ② and ③ of the shock tube for different equations of state

下载: 全尺寸图片 幻灯片

由图18可知，真实气体入射稀疏波后的驱动气体声速（图18中蓝色实线）在高压真实气体解线之前始终高于理想气体中相应部分的声速（图18中黑色虚线）。因此，当真实气体的左行稀疏波从端壁反射后，其传播速度仍然显著高于理想气体。同时，真实气体入射激波后声速则略低于理想气体；当右端反射激波与接触面作用后，根据气体动力学理论，由于真实气体中接触面左侧（区域③）气体有较高的声速，反射激波进入左侧气体后有更快的传播速度。

由于喉道处于壅塞状态，反射激波作用后接触面的运动速度很大程度上取决于驻室气流的状态。根据一维等熵流理论，壅塞流速大致正比于驻室区的声速（或温度0.5次方），由于理想气体中区域⑤的气流温度略高于真实气体的（图12），因此真实气体接触面的右行速度本来就略低于理想气体，但两者差异并不明显。实际流动中真正使得真实气体中接触面延迟抵达的原因是反射稀疏波提前抵达并作用于驻室区，使得当地声速急剧降低，从而极大降低了接触面的后续运动速度。

2.4   不同条件下高压效应对试验时间的影响

激波风洞试验往往需要在同一设备上获取不同总温和不同总压的试验气流。通过改变激波风洞运行条件，获得试验气流总温2000～6000 K、总压20～200 MPa范围内缝合接触面运行的不同试验工况，以此考察不同条件下高压真实气体效应对试验时间的影响。试验气流总温的调整和缝合条件的实现通过调整驱动段组分（氢气与氮气的配比）和高、低压段压比实现，试验气流总压的改变通过在维持高、低压段压比情况下调整压强实现。各段初始气体温度恒为300 K。

首先，在维持试验气流名义总温为6000 K的情况下，改变试验气流的总压，分别提取理想气体状态方程和PR状态方程下数值模拟获得的有效试验时间进行比较。不同总压条件下，入射激波、反射稀疏波头和接触面到达驻室的时刻（分别标记为t₁、t₂和t₃）如图19所示。

图  19  不同状态方程下各个波系到达时刻随总压的变化

Figure  19.  Arrival times of waves changing with total pressure under different equations of state

下载: 全尺寸图片 幻灯片

在理想气体状态方程下，由于激波管初始压比和温度维持不变，初始入射激波和稀疏波的速度均不变，而在湍流区激波沿程衰减幅度同样变化不大，因此，各波系的时空关系基本维持不变，反射激波、反射稀疏波和接触面的抵达时刻均无显著改变。此时各工况均为接触面率先抵达，有效试验由t₂−t₁决定。

在采用PR气体状态方程之后，如前面的分析：随着总压的提高（驱动压强增高），入射激波的强度基本未受影响（轻微减弱）；稀疏波传播速度加快，反射稀疏波头抵达时间提前；相应的，受提前抵达的反射稀疏波作用，驻室区气流声速降低，其经由喉道流出的速度降低，这使得接触面抵达时间延迟。从图19可以看到，在试验气流总压20 MPa运行条件下，反射稀疏波提前幅度不大，接触面仍在稀疏波前抵达端部，因此，接触面并未出现显著延迟，此时有效试验时间仍由t₂−t₁决定。而当总压提高至50 MPa，反射稀疏波越过接触面率先抵达端部，此时试验气流部分受到反射稀疏波作用，导致接触面相对理想气体情况发生滞后。这一效应随着试验气流总压的提升愈发显著。总压50 MPa以上的工况，真实气体有效试验时间转而由t₃−t₁决定。在试验气流总压100、150和200 MPa等3种工况下，高压真实气体效应导致有效试验时间相对理想气体分别减少约40%、50%和60%。

进一步对试验气流总温2 000～6 000 K、总压20～200 MPa范围的其余工况进行计算。为达到缝合接触面运行条件，不同试验气流名义总温对应激波管驱动段掺入氮气摩尔分数以及高、低压段的压比如表4所示。

表  4  激波风洞试验气流名义总温对应的初始压比和驱动气体氮气摩尔分数

Table  4.  Initial pressure ratios and molar fractions of nitrogen in driver gas corresponding to different nominal total temperatures of test flow in the shock tunnel

总温/K	初始压比	氮气摩尔分数/%
2000	176.5	15.0
3000	444.4	9.5
4000	800.0	5.0
5000	1272.7	1.8
6000	1600.0	0

下载: 导出CSV 

| 显示表格

将以上测试范围内高压真实气体效应对有效试验时间的影响情况总结为图20，其中Δ定义为2种状态方程下有效试验时间的差值相对理想气体方程有效试验时间的比：

图  20  不同总温、总压条件下高压气体效应对试验时间的影响

Figure  20.  Influences of the high-pressure gas effects on test time under different total pressure and temperature conditions

下载: 全尺寸图片 幻灯片

分析固定总压条件下，有效试验时间随总温的变化情况可知，随着总温升高，Δ逐渐减小。结合2.2节与2.3节中的分析，温度的影响主要体现在试验气流参数上，对于试验时间影响较小。Δ随总温升高逐渐减小的原因为高温对高压气体效应的影响具有抑制作用。随着温度的升高，高压气体效应引起的稀疏波运动速度变化幅度减小，在限制有效试验时间的因素为反射稀疏波到达后，理想气体状态方程和PR状态方程下有效试验时间的差异进一步减小。

如图20所示，不同驱动状态下高压气体效应的影响程度不同：(1)在试验气流总温较低时，试验时间主要受到稀疏波限制，高压真实气体效应致使有效试验时间直接呈现缩短趋势。(2)当试验气流总温提高至3 000 K以上，理想气体中接触面先于反射稀疏波抵达端部，此时如总压较低，高压效应不足以改变两者的抵达顺序，则有效试验时间几乎不变；而当总压提升，高压效应导致反射稀疏波先于接触面抵达端部，则有效试验时间随试验气流总压的提高而逐步缩短。

综上所述，对于高压氢气驱动的激波风洞内的全场流动时空结构和驻室区气流参数的分析获得如下认知：(1)高压真实气体效应对激波管流动的影响主要体现为稀疏波传播速度加快侵蚀有效试验时间长度；(2)同等压强下，温度越高则高压气体效应越弱。因此，为了减小高压真实气体效应的影响，工程中易于实现的方法有：①延长驱动段的长度，使得稀疏波需要经过更长的距离才能反射到达驻室，延长有效试验时间；②对激波管驱动段进行加热处理（同气体介质、同等压强下，高温气体驱动能力更强，气体密度更小，可削弱高压气体效应）。

3.   结　论

推导了立方型气体状态方程下气体的热力学参数（内能、比热容、比热比等）和声速，并将上述参数运用于耦合高温热化学非平衡过程的准一维数值模拟，对考虑高压真实气体效应的高焓激波风洞流动结构以及试验气流参数进行了准一维数值计算和分析，获得了如下主要结论。

(1)采用PR气体状态方程可以很好地刻画气体的高压真实气体效应，耦合该方程的准一维数值模拟可以准确复现高焓激波风洞的整体流动。

(2)高温真实气体效应和高压真实气体效应在物理上近乎处于相互排斥的地位，两者分别在激波管流动的热端（被驱动气体）和冷端（驱动气体）起主导性作用。高压真实气体效应对激波管的驱动能力以及入射激波和反射激波后的流动状态影响甚微，仅在驱动压强较高时轻微削弱所产生激波的强度；但在温度不高的驱动气体中，高压真实气体效应使气体声速增大，从而导致稀疏波和反射稀疏波传播速度加快，这一效应可以极大地改变激波管中各波系的相干时空关系。

(3)高压真实气体效应对激波风洞驻室试验气流的影响主要呈现为它对有效试验时间长度的侵蚀。这种作用主要通过高压效应加快稀疏波的传播、从而致使稀疏波提前抵达驻室区并显著降低当地气流总压达成。如高压效应不能使反射稀疏波先于接触面抵达喉道，则高压效应不对有效试验时间形成实质侵蚀。

(4)激波风洞在试验气流总温较低而总压较高的工况下，有效试验时间一般以反射稀疏波头抵达时刻为终点，高压真实气体效应更容易在此类工况下对有效试验时间产生大幅影响。对于采用冷高压气体驱动的高焓激波风洞，可在理想设计的基础上适当延长驱动段长度以抵消反射稀疏波提前抵达对试验时间的影响。