基于NACA0012翼型的超燃冲压高超声速飞行器的气动热预示

汪 钰,李 芃

（铜陵学院，安徽 铜陵 244000）

0 引言

高速气流在穿过激波时会被压缩，由此导致减速和温度急剧升高，产生的高温环境会引起飞行器烧蚀、干扰内部设备的正常工作。因此，对飞行器气动热的准确预示至关重要。风洞试验、工程估算和计算流体力学（CFD）是气动热预示的常用方法。第一种方法耗时长、成本高，第二种方法的估算精度有限。随着计算机技术的迅速发展，CFD已成为气动热预示研究的主要方法［1-3］。NA‐CA0012是研究飞行器外流场特性的典型翼型，它是由美国国家航空咨询委员会（National Advisory Committee for Aeronautics简称NACA）设计的。目前基于NACA0012的气动热研究主要集中在两个方面：（1）翼型的几何变化对温度或温度系数的影响［4-5］；（2）温度对翼型表面结冰特性的影响［6-7］。后缘形状是影响预示精度的重要因素［8］，根据建模方法的不同，NACA0012包括尖后缘和钝后缘两种。本研究设计了两种翼型外流场，通过与参考数据的对比，分析了不同参数对仿真预示的影响，给出了高超声速飞行条件下的最佳参数配置，并进行了相关的气动热预示，进一步丰富和完善了有关NA‐CA0012的研究成果。

1 仿真基础

1.1 超燃冲压高超声速飞行器

超燃冲压高超声速飞行器是指以超燃冲压发动机为主要动力的巡航飞行器，此类发动机只有在飞行速度达到一定值时才能启动，因此通常采用组合推进方案。如以火箭为基础的联合循环（RBCC）、以涡轮为基础的联合循环（TBCC）和固体火箭双模推进。其中，RBCC充分发挥了火箭发动机大推重比和冲压发动机高比冲的优势，是未来高超声速飞行器和单级再利用飞行器的理想动力系统。

为了正确地预示RBCC高超声速飞行器的气动热，需要分析飞行器的飞行路径并确定自由来流的初始条件。如图1所示，RBCC飞行器的飞行路径由爬升段、巡航段和攻击段组成。飞行器首先由载机运送到指定地点，启动固体火箭发动机将飞行器从低速模态带入亚燃冲压模态。该阶段推进效率最高，高度和速度的增量最大。最后进入超燃冲压模态，当飞行器平稳爬升到目标高度和马赫数后，进入保持等高等速飞行的巡航段，其飞行距离约占总航程的90%，是飞行路径的主要阶段。到达目标空域后进入攻击段，对目标实施精准打击。本文提出了适用于巡航段的气动热预示方案，巡航段的飞行高度约为30km，自由来流的初始条件如下：马赫数M为6.5，大气压力Pt为1179Pa，大气温度T t为226.65K，大气密度ρ为0.0181 kg/m3，大气动力黏度μ为1.475e-05 kg/（m·s）。

图1 RBCC超燃冲压高超声速飞行器的飞行轨迹

1.2 计算域外流场

如图2所示，基于NACA4 airfoil generator和air‐foil tools设计了尖后缘和钝后缘两种NACA0012翼型并绘制了计算域外流场，其中翼型的特征长度L=1m。对外流场的数值模拟，一般建议远场距离为翼型特征长度的12-20倍。建立远场距离为20m的C型网格拓扑以最小化近壁网格的偏度。尖后缘采用ICEM CFD进行网格划分时可直接折叠后缘处的块，如图2（a）所示。钝后缘保留后缘处的块，如图2（b）所示。

图2 （a）尖后缘计算域外流场；（b）钝后缘计算域外流场

1.3 湍流模型

对于涉及可压缩和热传递的流动，需要求解质量、动量和能量守恒方程。此外，还要求解湍流带来的输运方程。考虑到湍流方程无法通过直接数值模拟（DNS）求解，采用了目前应用最为广泛的雷诺平均可压缩NS方程（RANS）。K-epsilon、SST komega和Spalart-Allmaras是被广泛应用且相对准确的RANS湍流模型［9-11］。K-epsilon模型主要用于工业CFD仿真，该模型存在对边界层分离和逆压梯度不敏感的问题，所以在翼型外流场的应用并不广泛。Spalart-Allmaras是为航空应用设计的单方程模型，能够为受逆压梯度影响的边界层仿真提供良好的结果。SST k-omega是基于标准k-omega模型设计的双方程模型，其对自由来流的初始条件不敏感，且经过校准，能够准确计算光滑表面的流动分离。因此采用双方程的SST k-omega湍流模型。

2 验证测试

2.1 网格策略

因为巡航段的马赫数M大于5，此时的计算域外流场处于高超声速阶段，所以采用参考文献［12］中的数据来验证仿真参数。计算域外流场如图3（a）所示，其中粗黑点为坐标原点，a、b和c是输入边界INPUT boundary，e、f和g是壁面边界WALL boundary，d是输出边界OUTPUT boundary。湍流流动受翼型壁面的影响，因此，近壁面区域的网格划分十分重要，需要保证边界层的第一层网格处于sub viscosity layer，关键在于计算第一层网格高度y H。相关计算公式如式（1）：

其中R e是雷诺数，ρ、U、μ和L分别是自由来流的密度、速度、动力黏度和翼型的特征长度。

其中y P是第一层网格中间位置到壁面的距离，y H=2y P。气动热预示对壁面边界层的精度要求较高，一般建议近壁面网格y+≤1。验证测试的参数取值如下：R e=10e+06，L=1m，μ=4.46082e-5kg/（m·s），ρ=0.0247kg/m3，Pt=576Pa，Tt=81.2K，U=1806m/s，近壁面网格加密区如图3（b）所示。为了保证数值仿真结果的可靠性进行了网格无关性测试，如表1所示，测试了420 000、590 000和830 000三种网格数，P和T是翼型表面不同位置的压力和温度，830 000网格T/Tt的平均误差率优于590 000网格，但是后者的P/Pt的平均误差率数值和总平均误差率均优于830 000网格，因此采用590 000网格。

图3 （a）计算域外流场；（b）近壁面网格加密区

单元雷诺数［13］和近激波面的aspect ratio［14］是影响仿真准确度的两个重要网格参数。文献［13］建议对二维圆柱体，其单元雷诺数应不小于8。针对NACA0012翼型，对单元雷诺数分别取值4、8和16来寻找理想的数值，对应的y H和y+分别为（3.5e-6m 0.0825）（8e-6m 0.175）和（1e-5m 0.35）。基于网格总数和单元雷诺数不变的前提，再比较不同as‐pect ratio数值下的误差率，来确认最合适的网格参数，见表1。

表1 网格无关性

2.2 数值方法

表2给出了在高超声速阶段采用的数值仿真参数。求解器选用专为高速可压缩流设计的密度基，通量采用ROE和AUSM类型，其仿真结果与试验数据较为接近。SST k-omega湍流模型应用intensity and viscosity ratio turbulence method，其中的湍流强度（%）和湍流黏度比均为1。计算域的边界条件如下：输入边界INPUT boundary类型为pressure far field，velocity specification method选 择magnitude and direction，攻角为0°，自由来流X方向分量为cos（0°）=1、Y分量为sin（0°）=0。输出边界OUTPUT boundary类型为pressure outlet，壁面边界WALL boundary类型为no-slip，等温壁，初始值311K。梯度计算的常用方法有cell-based和node-based，第二种方法的准确度更高，尤其是针对具有倾斜度的非结构化网格，但是其对于硬件资源的要求更高。考虑到生成的结构化网格和有限计算资源，采用least-squares cell-based方案计算梯度，该方法可以提供与node-based相当的精度且计算开销较低。对于压力，采用基于中心差分的second-order方案。

表2 数值仿真参数

2.3 数值结果分析

基于一种湍流模型、两个通量类型和三个单元雷诺数，对每一种后缘翼型进行如下6组参数配置的 仿 真：SST+AUSM+yplus0.35、SST+ROE+yplus0.35、SST+AUSM+yplus0.175、SST+ROE+yplus0.175、SST+AUSM+yplus0.0825、SST+ROE+yplus0.0825。

对尖后缘，上述参数配置下的P/Pt分布如图4（a）所示，其中X轴为x/L，Y轴为P/Pt的数值。T/Tt的分布如图4（b）所示。P/Pt和T/Tt的平均误差率分别为（4.01%4.17%）（19.91%14.30%）（5.84%3.86%）（11.64%9.21%）（4.09%4.10%）和（12.96%12.88%）。如图4（c）所示，SST+ROE+yplus0.35配置下的总平均仿真误差率最低。接着基于该配置，对近激波面处网格的aspect ratio分别取值30、60、120、180、240和300再执行6组仿真，如图4（d）所示，对应的P/Pt和T/Tt的平均误差率分别为（5.06%10.05%）（4.01%4.17%）（4.74%3.98%）（5.47%3.10%）（4.71%3.14%）和（10.48%7.06%）。当as‐pect ratio取值较小时，压力和温度的平均误差率较大。之后随着aspect ratio的增大，平均误差率开始下降。当aspect ratio等于240时可获得最佳结果，总的平均误差率约3.93%。

图4 （a）尖后缘P/Pt的分布图；（b）尖后缘T/Tt的分布图；（c）尖后缘在不同参数配置下的平均误差率；（d）尖后缘在不同aspect ratios下的平均误差率

对钝后缘，上述参数配置下的P/Pt分布如图5（a）所示，其中X轴为x/L，Y轴为P/Pt的数值。类似的，T/Tt的分布如图5（b）所示。P/Pt和T/Tt的平均误差率分别为（5.67%3.35%）（4.22%4.14%）（7.94%3.08%）（8.91%6.70%）（5.27%3.36%）和（9.47%6.02%）。如图5（c）所示，SST+AUSM+yplus0.35配置下的平均误差率最低。接着基于该配置，对近激波面处网格的aspect ratio分别取值40、80、160、240、320和400再执行6组仿真，如图5（d）所示，对应的P/Pt和T/Tt的平均误差率分别为（19.64%7.75%）（4.22%4.14%）（3.67%3.28%）（5.91%10.03%）（3.06%2.83%）和（2.74%9.58%）。当aspect ratio取值较小时，压力和温度的平均误差率较大。之后随着aspect ratio的增大，平均误差率开始下降。当as‐pect ratio等于320时可获得最佳结果，总的平均误差率约2.95%。

因此，对NACA0012翼型的后缘形状，钝后缘的性能要优于尖后缘，前者较后者带来的仿真精度提升了约24.97%，后者较前者导致仿真精度下降了约33.28%，后缘翼型选择不当会导致较大的仿真误差。所以在高超声速阶段，建议采用的参数配置如下：钝后缘翼型+y+取值0.35+单元雷诺数为16+近激波面aspect ratio为320+密度基求解器应用AUSM类型的通量。如果计算域外流场选用了尖后缘翼型，那么建议采用如下参数配置：y+取值0.35+单元雷诺数为16+近激波面aspect ratio为240+密度基求解器应用ROE类型的通量。

图5（a）钝后缘P/Pt的分布图；（b）钝后缘T/Tt的分布图；（c）钝后缘在不同参数配置下的平均误差率；（d）钝后缘在不同aspect ratios下的平均误差率

3 气动热预示

基于NACA0012钝后缘翼型建立了远场距离为20m的C型计算域外流场，网格y+取值0.35（对应y H=1.6e-5m），单元雷诺数为16，近激波面aspect ra‐tio取值320。高超声速飞行器巡航段的气动热预示方案仿真参数如表3所示，在巡航段飞行器保持6.5马赫的恒定速度和30km的恒定高度，自由来流的初始压力和温度分别为1179Pa和226.65K。初始化选用基于拉普拉斯方程的hybrid非均匀方法，稳态分析，计算迭代5000次。气动热温度最终达到了2158K（1885℃），预示结果表明RBCC高超声速飞行器在飞行过程中将承受极端的气动热问题。

表3 数值仿真参数

4 结论

基于NACA0012翼型，采用文献［12］的初始条件和数据进行了验证测试，分析了不同后缘形状和网格参数（单元雷诺数和近激波面aspect ratio）在高超声速条件下对数值仿真精度的影响，给出了最优参数配置并对RBCC高超声速飞行器在巡航段的气动热进行了预示，结论如下。

（1）钝后缘的仿真性能优于尖后缘。

（2）不同于二维圆柱体，对NACA0012翼型，单元雷诺数的建议取值应当不小于16，近激波面网格的aspect ratio数值取决于选用的后缘形状。对尖后缘翼型，aspect ratio建议取值为240。对钝后缘翼型，aspect ratio的建议取值为320。

（3）在20米远场距离条件下采用钝后缘+SST k-omega+单元雷诺数16+近激波面aspect ratio320+AUSM通量类型的密度基求解器可获得最佳仿真效果，P/Pt和T/Tt的平均误差率分别为3.06%和2.83%。RBCC高超声速飞行器在飞行过程中的气动热可达到1885℃，热障问题十分突出。