高焓风洞一体化数值模拟及其对气动特性影响

傅杨奥骁, 董维中, 丁明松, 刘庆宗, 江 涛

(中国空气动力研究与发展中心 计算空气动力研究所, 绵阳 621000)

0 引 言

高超声速飞行器在再入和滑翔过程中，如果飞行速度很高(一般认为马赫数10以上)，会出现高温气体非平衡效应[1]，这种效应会对飞行器的气动特性产生很大影响。目前对于高温气体非平衡效应的研究大多结合试验进行，由于飞行试验成本较高，目前试验研究大多是在地面高焓试验设备中进行的。近几十年来高焓试验设备的研制获得了高度的重视，世界各地发展出了各种不同的高焓试验设备[2-4]。在众多种类的设备中，高焓激波风洞可以模拟的总温总压高，运行成本低，是目前国际高焓流动试验研究应用的主力试验手段[5]。

高焓激波风洞通过强激波压缩得到高温高压状态的试验气体，然后通过喷管的膨胀加速得到超高速试验气流。形成强激波一般有多种方式，常见的有加热轻气体、燃烧加热、变截面、爆轰驱动等[6]。20世纪90年代，俞鸿儒[7]等提出了通过爆轰驱动获得强激波和高焓试验气流的理论，在此基础上设计的JF-10氢氧爆轰驱动激波风洞取得了巨大成功。风洞利用驱动段爆轰燃烧产生的高温高压气体作为驱动气体，为国际首创，具有有效试验时间长、耗费低、运行平稳、复现性好的特点。风洞可以获得超高速高焓试验气流，具有模拟高温气体非平衡效应的能力。

驻室中的试验气体可认为是处于热化学平衡态的多组元混合物。但随着气体在喷管中的不断加速，特别是经过喉道后的迅速膨胀，流动来不及达到当地温度、压力下的热化学平衡，喉道下游还会出现明显的热化学冻结现象[8-9]。在这种情况下，高焓风洞试验数据的分析和使用，与常规风洞试验相比，要困难和复杂得多，一般需要借助数值计算的方法进行研究。

在现有的高焓风洞数值研究中，常常将喷管流动和模型试验解耦处理，即先通过试验测量或者数值模拟等方式得到风洞喷管出口的参数(或者略作处理，如平均[10-11]、外延[12]等)，然后将其作为试验模型的均匀来流开展研究。通过这种解耦处理，可以相对独立地进行喷管流动与模型试验研究，因而应用十分广泛。曾明等在2006年[10-11]对JF-10激波风洞流场的参数测量进行了数值重建；2013年，Clemente等[12]对意大利SCIROCCO等离子风洞中的EXPERT太空舱试验进行了数值模拟；Ishihara等在2013年[13]和2014年[14]对日本HIEST自由活塞激波风洞中的钝头体试验进行了数值模拟；2015年，Holden等[15]综述了在LENS系列激波风洞中进行的试验及其数值模拟；2016年，Vasilevskii等[16]对俄罗斯VAT-104风洞的钝锥体试验进行了数值模拟。尽管解耦方法应用较多，但它在试验模型流动研究中忽略了喷管出口和试验段流场核心区参数的非均匀性，不能较为真实地反应高焓风洞试验的过程。

本文以更加真实地反映风洞试验物理过程为目的，基于一体化数值模拟思路，考虑喷管出口和试验段流场的非均匀性，开展高焓激波风洞JF-10典型运行状态下喷管/试验段/模型流场的一体化数值模拟，分析喷管出口和试验段流场非均匀性的产生原因及其对试验模型流场特性的影响，并与两种喷管/试验段模型流场解耦求解的方法进行对比，分析多个模型位置、攻角试验条件下不同计算方法对流场参数分布、试验模型表面热流分布、模型气动力系数的影响。

1 计算方法及物理化学模型

1.1 控制方程及求解方法

风洞流场包含明显的热化学冻结现象，因此流场控制方程采用包含热化学非平衡源项的三维Navier-Stokes方程，其无量纲形式为：

(1)

Q=(ρi,ρEV,ρ,ρu,ρv,ρw,ρE)T

W=(wi,wV,0,0,0,0,0)T

式中Q是守恒变量，Re为雷诺数，F、G、H和FV、GV、HV分别对应x、y、z三个方向的对流项和黏性项，W为热化学非平衡源项。计算采用结构网格，控制方程离散采用有限差分方法，对流项通过AUSMPW+(Advection Upstream Splitting Method by Pressure-based Weight functions)格式离散，黏性项通过中心差分格式离散，时间项采用LU-SGS(Lower-Upper Symmetric Gauss Seidel)格式。非平衡源项、对流项和黏性项均采用全隐式处理。具体计算公式和方法详见文献[17-18]。

1.2 热化学反应模型

由于风洞试验气体为空气，因此本文选用7组分(O2、N2、NO、O、N、NO+、e)Dunn-Kang化学反应模型[18]；热力学模型采用Park两温度模型[19]，考虑振动非平衡效应。具体计算公式和方法详见文献[18]。

1.3 表面催化特性模型

2 计算方法校验

针对JF-10激波风洞，开展了喷管/试验段一体化数值模拟研究，并与试验及文献结果进行对比分析。JF-10激波风洞的喷管为锥形喷管，半锥角为7.1°，出口直径为500 mm，喉道直径为11 mm，总长为2 m。试验段长2 m，内径1200 mm，计算几何示意图如图1所示。驻室总压为19.6 MPa，总温为7920 K，试验气体为空气。壁面条件采用等温壁Tw=300 K条件、完全催化条件(FCW)，入口条件采用亚声速入流边界条件。

图1 计算几何示意图

图2(a)和图2(b)给出了沿喷管轴线的温度和部分组分质量分数分布与文献[10]结果的对比，横坐标为0处为喷管喉道，文献采用的是二维轴对称计算，而本文为三维计算。

(a)温度分布

从图2中可以看出，本文计算得到的平动温度、振动温度、氧气、氮气等典型参数分布与文献[10]符合良好。同时可以看出，在喷管出口处，温度(平动-转动温度)约450 K，振动温度在3500 K以上，存在热力学振动温度的冻结现象；温度450 K的平衡态空气中不会有O、N、NO等组分[10,18]，而在喷管出口处，气流中仍存在一定浓度的O、N、NO等组分，可见，由于流速较高，气流出现化学反应的冻结现象。

图3给出了激波风洞中皮托管的数值重建结果与试验测量结果对比。皮托压力传感器安装在半径Rn=0.0095 m的球头圆柱体顶点，通过计算风洞中三维球头绕流流场，开展皮托压力测量的数值重建。以喷管/试验段一体化数值模拟中得到的皮托管安置处的自由来流参数作为来流条件，计算球头热化学非平衡流场，采用等温壁、完全催化壁(FCW)条件。图3(a)给出了d=100 mm处数值重建的压力分布云图，图3(b)给出了风洞轴线上皮托压力测量的数值重建结果与风洞试验结果的对比，其中横坐标d表示皮托管距喷管出口的距离。可以看出，数值重建结果与试验结果符合良好。

(a)球头压力分布云图

总的来说，计算结果与试验及文献结果吻合良好，计算方法和程序可信。

3 计算结果及分析

针对JF-10高焓激波风洞，开展了喷管/试验段/试验模型一体化数值模拟。计算域及计算网格如图4所示，试验模型为钝锥体，模型总长为224.4 mm，底部半径为75 mm，头部半径为50 mm，半锥角为6.9°，图中d表示试验模型距喷管出口距离。驻室总压19.6 MPa，总温7920 K，试验气体为空气，试验模型壁面温度采用等温壁Tw=300 K条件，表面催化特性采用完全催化(FCW)和非催化(NCW)条件。

(a)计算域示意图

图5给出了未放置模型时风洞流场典型参数分布，图5(a)、(b)分别为沿流向温度和速度、马赫数分布，其中横坐标x<2 m区域为喷管，x>2 m区域为试验段。由图可以看出，在试验段，高超声速气流沿流向存在一定的非均匀性：沿流动方向，平动温度有一定下降，马赫数缓慢升高。事实上，这种沿流向的非均匀性在试验中也得到了验证，从图3(b)的皮托压力试验测量可以看出，皮托压力从喷管出口(d=0 mm)处的21kPa一直下降到距喷管出口d=600 mm处的12.4 kPa。

图5(c)给出了本文计算得到的喷管马赫数分布云图，由图可以看出，在喷管出口附近，存在y方向的非均匀性，喷管轴线附近区域的马赫数略高于远离轴线区域。图5(d)为试验段中距喷管出口d=200 mm处沿径向静压分布，其中纵坐标y表示距风洞轴线的距离。由图可以看出，在试验段核心区域，气流也存在y方向的非均匀性，中心轴线区域(y=0 m)静压较低，在两端界面附近(|y|=0.22 m)出现局部峰值，最大相差12.6%。

试验段流向(x方向)和径向(y方向)的非均匀性产生原因，可结合图5和图6从两个方面进行分析：(1)由于喷管扩张段存在一定的扩张角，气流膨胀扩张，喷管出口处气流存在非x方向速度，从图6(a)给出的喷管出口处无量纲径向速度(通过轴向速度Vx无量纲化)沿径向分布以及图6(b)给出的喷管径向速度分布云图可以看出，x-y平面，在y=0.2 m附近，气流的y方向速度可达到500 m/s左右，约为轴向速度的10%，这种非x方向速度，会使气流在试验段呈扩张趋势；(2)喷管壁面不可避免的存在黏性效应，在喷管扩张段，沿x正方向，马赫数逐渐增大，黏性效应逐渐明显，壁面附近温度由于摩擦显著升高，边界层不断增厚，形成高温区域，使y方向的压力、温度、马赫数等参数分布不均匀，这种非均匀性会随着流动向下游扩展，影响区域逐渐增大。

(a)温度与速度分布

图6 喷管中垂直速度分量分布

由于这种流动的非均匀性，试验模型处于试验段的不同区域时，其来流条件和流场流动特性将会发生变化，会对其气动力/热特性的预测结果产生一定影响。而传统的解耦方法，忽略了流动的非均匀性，有可能带来预测误差。下面将针对这一现象，具体展开分析。

3.1 不同模型位置下的差别

在实际试验中，由于试验模型尺寸或者测量点设置等因素，模型在试验段中的安放位置有多种可能。这里我们将试验模型距喷管出口距离d分别取为0 mm、100 mm、200 mm和300 mm。考虑试验段流动的非均匀性，将喷管/试验段/试验模型一体化数值模拟，记为Coupled。传统的解耦方法，取喷管出口中心处参数，作为试验模型的均匀来流参数进行数值模拟，记为Case1，该方法忽略了试验段喷流流场流向和径向的非均匀性。另一种解耦方法，取试验模型所处位置的波前对称轴线处参数，作为试验模型均匀来流参数进行数值模拟，记为Case2，这种解耦方法忽略了试验段喷流流场径向的非均匀性。

为了排除网格差异的影响，Case1、Case2方法的计算网格与一体化模拟Coupled方法中,试验模型部分的网格保持一致。同时还开展了两套网格的计算研究：采用试验模型部分流向、法向、周向网格点数分别为141×61×31和141×81×41，最小网格间距分别为5×10-6m和2×10-6m的网格进行一体化计算，记为Grid1和Grid2，其热流分布如图7所示，可以看出，两套网格的热流值变化不超过2%。计算结果受网格影响较小，因此本文计算中如无特殊声明，均采用Grid2网格。

图7 不同计算网格的模型表面热流分布对比

图8给出了模型距喷管距离d=200 mm时放置模型前后的马赫数分布云图。可以看出，置入试验模型，在试验模型上游区域，流场变化幅度较小；而在试验模型下游，弓形激波与试验段超声速气流外边界存在强烈的相互干扰，气流膨胀扩张，流动结构十分复杂。

图8 马赫数分布云图(d=200 mm)

图9和图10分别给出了d=200 mm时模型驻点线上的温度和部分组分质量分数分布。可以看出，在大体上，一体化计算(Coupled)得到的流场参数分布，与解耦方法计算(Case1、Case2)结果基本一致，但在局部细节上存在差异：解耦方法得到的波后温度、分子组分质量分数较一体化方法稍低，而原子组分则稍高；解耦方法Case2较Case1，计算结果更接近一体化计算；解耦方法Case2和Case1计算结果差异较小，可以看出，这种条件下，均匀来流的绝对值细微区别(即流向非均匀性)对于驻点线参数分布影响不大。

图9 模型驻点线上的温度分布(d=200 mm, FCW)

图10 模型驻点线上的部分质量组分分布(d=200 mm, FCW)

图11给出了不同模型位置头部驻点热流的对比。由图可以看出：采用一体化计算(Coupled)得到的头部驻点热流随着距离d的增加，逐渐下降，在完全催化条件下(FCW)，其下降幅度更显著；与一体化计算相比较，采用传统的解耦方法Case1计算得到头部驻点热流随着距离d的增加，偏差逐渐增大，在d=300 mm时偏差可达10%左右；采用解耦方法Case2，其计算结果也与一体化计算也存在一定差异，但同Case1方法相比，其偏差幅度整体上小一些；解耦方法Case2和一体化计算Coupled差别小于解耦方法Case1与Case2差别，这说明对于头部驻点区域而言，喷流流向非均匀影响大于径向非均匀性影响。

图11 不同位置下的驻点热流对比

图12进一步给出了d=200 mm时模型表面热流分布。可以看出：沿流动方向，一体化计算与解耦计算的热流差异呈增大趋势，在身部区域热流值最大可相差20%；Case2较Case1更接近Coupled的结果；在靠近底部的锥身区域，解耦方法Case2和一体化Coupled的计算差别，大于解耦方法Case1与Case2的差别，造成这种现象的原因，是由于在身部区域，解耦方法Case1与Case2的差别主要由流向非均匀性引起。而一体化计算Coupled与解耦方法Case2的差异，既包含了径向非均匀性的影响，同时也包含了由于模型身部在流向位置上的差异引起的影响。

图12 模型表面热流分布(d=200 mm)

图13给出了不同模型位置时试验模型的阻力系数CD。可以看出：一体化计算结果(Coupled)与解耦计算(Case1、Case2)存在一定差异，解耦计算得到的阻力系数偏大，最大可达5%左右；随着距离d增大，一体化计算Coupled得到的阻力系数呈增大趋势，解耦方法Case1与一体化计算结果差异减小，而解耦方法Case2与一体化的差异基本保持不变；解耦方法Case2和一体化Coupled之间的计算差别，大于解耦方法Case1与Case2之间的差别，造成这种现象的原因，是由于流向非均匀性导致的Case1与Case2的差异，只是这种来流参数绝对值的差异对于无量纲的气动系数影响较小，而一体化计算中考虑了流动的扩张，扩张导致的径向非均匀性影响了流场结构、压力分布，因此造成了更大的气动力系数差异。

图13 不同位置下的阻力系数对比

模型表面压力的差别，是气动力系数计算差异的主要来源之一。图14(a)给出了本文计算的模型表面的压力分布(d=200 mm)。图14(b)给出了文献[12]的计算结果(文献的计算条件和外形，与本文不一致)，图中红线为解耦计算结果，蓝线考虑了试验段流场的非均匀性。可以看出，考虑试验段流场的非均匀性时，计算结果与试验结果符合更好。对比图14(a)和图14(b)，可以看出，尽管本文计算条件和外形与文献不一致，但其规律变化相同：采用解耦方法，计算得到的模型表面压力偏高，沿流动方向，其差距呈增大趋势。由此可见，要准确模拟高焓风洞试验段的非平衡流动，得到较为精确的气动力/热特性，必须考虑试验段流场的非均匀性，采用一体化或类似的数值模拟方法。

(a)d=200 mm

3.2 不同模型攻角下的差别

保持试验模型位置不变，模型距喷管出口距离d=200 mm，试验模型飞行攻角分别取为0°、5°、10°、15°。

图15给出了不同飞行攻角时模型头部驻点热流值和迎风面肩部顶点的热流值对比，图中Coupled、Case1和Case2含义与3.1节相同。可以看出，随攻角增大，一体化计算(Coupled)的头部驻点热流小幅下降，在15°攻角时下降了5%左右，这是由于径向的非均匀性造成的，而解耦计算(Case1、Case2)的热流值则基本保持不变；三者的迎风面肩部热流呈上升趋势；解耦方法得到的热流要高于一体化计算结果，二者的差距随着攻角增大而增大，在15°攻角时Case1与Coupled差别可达23%；Case2热流值较Case1更接近于Coupled。

(a)驻点热流对比

图16进一步给出了FCW条件、10°攻角下的试验模型表面热流分布，图中Windward表示迎风面，Leeward表示背风面。可以看出，模型攻角的变化使模型表面热流分布发生了明显变化，随着攻角增大，迎风一侧的热流升高；在头部区域，一体化计算热流值在Case1与Case2之间，更接近Case2，而Case1与Case2差异明显，这说明该区域流向非均匀性影响较大；而在靠近底部的锥身区域，解耦计算热流结果要高于一体化计算结果，Case1与Case2的差异小于Coupled与Case2的差异，这与图11和图12现象一致。

图16 模型表面热流分布(α=10°)

图17给出了不同攻角下得到的气动力系数对比，其中CL为升力系数，CMZ为俯仰力矩，其中力矩参考点为头部驻点。可以看出：解耦计算方法和一体化计算方法结果存在一定差异，随攻角增大，差异呈增大趋势；解耦方法Case1与Case2差别相对较小，这说明对于该条件下的气动力特性而言，喷流径向非均匀影响大于流向非均匀性影响，这与图13结论一致。

图17 不同攻角下的气动力系数对比

4 结 论

本文考虑风洞试验中的高温空气化学反应，振动激发及非平衡效应，通过数值求解三维热化学非平衡Navier-Stokes方程，完善了高焓激波风洞流场一体化数值模拟计算方法和计算程序。考核校验表明：数值计算结果与文献及试验结果符合良好，数值模拟结果可信度高。

1)通过计算发现，试验段流场在流向和径向存在一定的非均匀性。试验模型处于试验段的不同区域时，其流场流动特性将会发生一定变化，对试验模型气动力特性、气动热环境的预测结果产生影响。在进行数值模拟和试验方案设计时有必要考虑试验段流场的非均匀性。

2)通过不同解耦方法(Case1、Case2)与一体化数值模拟结果对比可以看出：流场非均匀性造成的气动力热差异，随着模型与喷管出口距离d增大，呈增大的趋势；对试验模型不同区域的影响存在差别，在驻点区域，流向非均匀性影响较大，在身部区域，径向的非均匀性影响较大。

3)解耦计算(Case1、Case2)与一体化计算的差异在不同攻角下的差别，主要体现了流场径向非均匀性的影响。随飞行攻角增大，差异整体上呈增大趋势；三者在身部区域差别较大，驻点附近区域相对较小，在迎风面的肩部顶点，最大差别可达20%以上。

从本文的结果来看，高焓风洞试验段流场的非均匀性对于模型气动力/热的影响，与试验模型的安装位置、攻角、外形相关，为了得到更为准确的试验数据，需要结合具体试验的一体化数值模拟的结果，对试验方案以及试验数据进行修正。基于一体化数值模拟的思想，数值计算可以更加深入地理解高焓风洞模拟特点，改善风洞试验数值模拟精度。

本文的数值模拟方法还可以应用于JF-12激波风洞[20]等其他高焓设备的模拟研究。