基于高分辨率格式的高超声速气动热数值模拟研究
作者： 孙振生 王江彬 胡宇 毛凯
摘 要：高超声速气动加热问题一直是高超声速飞行器设计中的关键，是气动布局设计、结构设计以及热防护选材等的重要依据。然而气动热的数值模拟受到网格尺度、数值格式等多种因素的影响，是目前计算流体力学的难点之一。本文选择典型的三维钝锥模型，主要围绕网格尺度、数值格式精度以及耗散对高超声速气动热问题展开相关研究。结果表明，随着网格雷诺数的不断增加，热流计算值呈现不断减小的趋势，在气动加热比较严重的部位，热流计算值对于网格尺度的敏感性更强;网格尺度较小时精度越高，热流计算越准确，而网格尺度较大时，具有高分辨特性的MDADF-HY格式更为占优;同时，适当的耗散也是数值格式能够准确计算热流的前提。
关键词：高超声速;气动加热;网格尺度;格式精度;格式耗散;飞行器
中图分类号：TJ760; V211.3
文献标识码：A
文章编号：1673-5048（2022）04-0064-06
DOI：10.12132/ISSN.1673-5048.2021.0237
0 引 言
当飞行器在高空以超声速或高超声速飞行时，会产生严重的气动加热现象。为解决高超声速飞行器研发设计中热防护等问题，迫切需要发展能够准确预测气动热环境的工具。近年来，针对高超声速气动热环境预测的CFD数值模拟方法得到了迅速的发展，但也存在很多尚未解决的问题。
对于气动力中升力的计算，目前的大部分CFD方法均可以得到相当精确的结果，但是针对气动热的计算，不同的CFD方法得到的结果可能会有较大的偏差。网格的划分、计算格式的精度及耗散、湍流模型的选择等诸多因素都会对气动热的数值模拟结果造成影响，甚至会有量级上的差别。目前为止，用数值方法求解气动热问题一直是CFD研究中非常困难的一个问题。随着国内航空航天领域技术的发展，针对高超声速气动热的数值模拟方法必然会得到更为广泛的应用，解决数值计算中的气动热问题对于高超声速飞行器的研发和设计具有很强的现实意义。
20世纪初，普朗特边界层理论的提出为气动热的数值计算提供了非常可靠的理论工具。文献[1]通过求解边界层方程得到了攻角在25°～40°之间的热流分布。文献[2]采用粘性激波层数值计算方法计算了STS-2E航天飞机在某些飞行状态下的气动热流分布，与试验数据吻合效果良好。随着一系列差分格式的提出，针对气动加热问题的数值计算方法研究得到迅速的发展。
国内对于高超声速气动加热数值计算的研究起步相对较晚，主要是通过求解N-S方程来求解高超声速气动加热问题，有力地推动了高超声速数值计算方法的发展。文献[3]基于黎曼问题解析解构造了低耗散NNLD格式，通过计算气动力热问题验证了该格式的可靠性。文献[4-5]采用有限体积法，无粘通量采用三阶MUSCL格式以及AUSMPW+格式进行计算，时间离散采用LU-SGS方法，在湍流计算中添加了先进壁面函数边界条件，该壁面条件考虑了热传导效应，最后对高超声速飞行器进行了数值模拟，得到较为可靠的气动热数据。文献[6]在WCNS格式的基础上，发展了高精度算法WCNS-E-5，并对钝锥以及双椭球进行了热流计算，通过与二阶MUSCL格式对比，证明了WCNS-E-5具有更高的流线分辨率，得到的热流密度更加准确。也有一些学者将有限元积分方法应用到高超声速流场问题中，文献[7]给出一种有限元-有限差分混合算法，针对高超声速钝头体求解N-S方程，得到比较满意的结果。
为了抑制在间断附近由于梯度过大引起的非物理振荡，许多格式会使用限制器以限制插值梯度。文献[8]采用minmod， Van leer和Osher-C三种不同的限制器对典型的高超声速飞行器外型双锥模型进行了气动热数值模拟，得出minmod限制器对于热流的计算在三种限制器中性能最优。文献[9]分析对比了传统的MUSCL限制器以及多维限制器的计算性能，发现多维限制器的鲁棒性和准确性要更优。文献[10]提出了一种新型三阶TVD限制器，并将其应用于复杂外形气动热计算，表现出良好的计算精度以及间断分辨率。
除数值格式以外，高超声速气动热数值模拟若想得到精确度较高的热流结果，拓扑结构设计合理以及间距尺度控制良好是计算网格必不可少的条件。关于网格尺度对气动热数值模拟的影响，国内外学者开展了广泛研究。文献[11]在计算二维圆柱问题时发现，根据网格雷诺数原则确定的第一层网格尺寸会偏大，需要减小一个数量级才能够满足计算精度要求。文献[12]对二维圆柱模型进行了气动热数值模拟，研究了网格雷诺数对热流沿壁面分布的影响，认为热流的准确计算需要满足网格雷诺数小于10，而过小并不会明显改善计算结果。
在计算网格生成方面，许多学者提出了值得借鉴的准则。文献[13]将壁面热流与气体分子平均自由程联系起来，针对近壁面网格提出了基于分子平均自由程的尺度准则，网格的生成只依赖于壁面局部参数，通过计算钝锥算例验证了该准则的可靠性。然而其局限性在于壁面局部参数是未知的，需采用粗网格试算获得。针对该问题，文献[14]提出了壁面参数预估方法，无需为获得壁面参数而进行提前试算，提高了该准则的实用性。同时计算了完全气体以及真实气体下的高超声速流动问题，验证了该方法的可靠性。文献[15]认为壁面法向网格尺度应该与当地温度梯度有关，因此可以预先应用工程算法估算出壁面温度梯度，之后根据估算值确定壁面法向网格尺度，从而生成计算网格。此种方法不仅能够提高计算精度，而且具有加速收敛的效果。
目前对于气动热的数值模拟，采用不同数值方法的精度、耗散特性等对热流的影响机理并不十分清楚，尤其对于格式的耗散特性在气动热方面的研究非常少。本文采用最新提出的高分辨率格式MDADF-HY对热流的网格尺度效应进行研究，并结合传统的数值格式分析了格式精度对热流的影响，最后采用一类在不影响色散情况下单独控制耗散的格式来研究耗散对热流的影响。
1 数值方法与计算模型
1.1 量热完全气体控制方程
本文的热流是通过求解可压Navier-Stokes方程来进行计算。在不计质量力和源项的情况下，Navier-Stokes方程在直角坐标系中可以写为下列向量形式：
具体参数表示可以参考文献[16]，壁面热流q为
数值离散方法采用有限体积方法，对于粘性项的离散采用中心格式;对于无粘项的离散，采用4种不同的数值格式，分别是二阶MUSCL格式、五阶WENO格式、MDCD格式以及MDADF-HY格式[17-19]。在获得单元界面两侧变量的基础上，通过Roe格式近似求解黎曼问题得到单元界面处变量值。文献[20]指出该算例在此流动条件下为层流，因此无需考虑湍流模型对计算带来的影响。时间推进采用三级三阶Runge-Kutta方法。
1.2 模型选取及计算条件
三维球头钝锥是十分典型的高超声速飞行器再入体模型，本文的三维钝锥算例选自NASA TN D-5450报告[21]，采用热电偶对不同钝度的钝锥模型测试了多种来流雷诺数及攻角情况下的热流值，试验数据丰富。其中驻点热流与理论相比测量误差小于±5%;对于0°攻角，锥面上的三条母线上的热流测量重复性误差小于±6%;由于精度是相对于热流测量水平来说的，因此对于低热流的测量误差估计在±20%以内。本文选择其中一种几何外形，如图1所示，钝锥半锥角θc=15°，球头曲率半径R=0.009 5 m，锥尖到锥底的总长度L=0.568 7 m。
计算状态：来流马赫数为10.6，来流雷诺数Re∞=3.937×106，温度T∞=47.34 K，壁面温度Tw=294 K。选取攻角为10°以对比分析迎风面以及背风面流场信息。由于热流稳定要比压力缓慢得多，因此按照文献[22]给出的热流稳定步数在压力稳定步数的10倍以上作为收敛准则。
2 网格雷诺数对热流的影响分析
网格雷诺数的定义式为Recell=ρ∞u∞Δnwμ∞，其中Δnw为壁面法向第一层网格高度，其余物理量为来流参数。图2所示为钝锥整体计算网格。
采用MDADF-HY格式计算无粘通量，壁面热流计算结果按无量纲形式Q/Q0给出，其中Q为计算得到的热流值，Q0为文献[21]给出的驻点热流值，对Q统一作归一化处理。横坐标为至尖锥（非钝体头部）距离与钝锥总长度之比。取0.1，1，10，100，1 000五种不同网格雷诺数计算结果进行分析。图3分别给出了攻角为0°以及10°时五种不同网格雷诺数情况下钝锥表面热流沿轴向的分布，并且与试验值做出了对比。从图中可以得出，五种网格雷诺数计算得到的壁面热流大致有相同的趋势，头部驻点区由于受到气流滞止的作用以及正激波的影响，气动加热效果非常严重，导致热流密度非常高。离开驻点区域后，热流密度迅速下降，在相对位置大概为0.2之后，热流密度下降的趋势逐渐变缓，到达钝锥尾部时基本保持稳定。
随着网格雷诺数的不断增加，热流计算值呈现不断减小的趋势。对于10°攻角背风面，气动加热现象不明显，五种网格雷诺数对应的网格得到的热流计算值均能与试验值较好的吻合。但对于迎风面，网格雷诺数为1 000时得到的热流计算值与试验值差距比较大，尤其是在头部附近，最大误差已经超过了50%。可见网格尺度过大时，计算的热流值难以升上去，在钝锥大部分区域维持在较低的范围，远远达不到对气动热数值模拟的要求。
对比0°攻角以及10°攻角下的壁面热流值，在一定的攻角下，钝锥迎风面的热流密度更高，并且要远远大于背风面热流值，因此在热防护设计中，飞行器迎风面是需要重点考虑的部位。
3 格式精度对热流的影响分析
CFD计算中，数值格式是关系到数值模拟结果准确与否最主要的因素。对于高超声速气动热来说，热流值是与温度梯度相关的物理量，不同精度插值用到的模板点数量也不同，这就会导致计算得到的温度值及温度梯度值的精确度不同。本文选择三种不同精度的格式：二阶MUSCL格式、四阶MDADF-HY格式及五阶WENO格式，来研究格式精度对壁面热流计算的影响。
为对比分析三种格式精度对热流的影响以及对网格雷诺数的敏感性，针对0°攻角工况，三种格式的驻点热流计算值与Fay-Riddell公式计算值[23]相对误差随网格雷诺数的变化曲线如图4所示。
可以看出，当网格雷诺数足够低时，在0.1～10范围内，三种格式得到的壁面热流相对误差都能够收敛到某个确定值，精度越低得到的壁面热流相对误差越大，并且高阶格式计算出的热流值要低于低阶格式。比较四阶MDADF-HY格式和五阶WENO格式，当网格雷诺数较低时，WENO格式得到的热流值要比MDADF-HY格式得到的热流值更接近试验值。但是随着网格雷诺数的逐渐增大，到10的量级之后，MDADF-HY格式的计算值反而更接近试验值，分析认为MDADF-HY格式的分辨率比WENO格式更好，当网格比较稀时，MDADF-HY格式对于粘性边界层的模拟以及激波等间断结构的捕捉效果要更好。
同时，随着网格雷诺数由高到低，相比于四阶MDADF-HY格式，五阶WENO格式虽然收敛到最小相对误差时所需的网格雷诺数更小，但是达到收敛时计算结果更准确。分析认为，高阶插值求解壁面热流，在计算无粘通量时需要在控制体单元两侧采用更多的模板点，利用了更多的流场信息来插值计算界面处通量，因此得到的结果更加准确，对于壁面热流的模拟更逼近试验值。除此之外，二阶MUSCL格式由于精度比较低，对于热流计算的相对误差一直大于另外两种格式，并且MUSCL格式采用的单元两侧的模板点比较少，导致其对于网格雷诺数更加敏感，需要更密的网格才能收敛到相对误差最小值。