基于高分辨率格式的高超声速气动热数值模拟研究

孙振生，王江彬, 胡 宇，毛 凯

(火箭军工程大学，西安 710025)

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方程在直角坐标系中可以写为下列向量形式：

(1)

具体参数表示可以参考文献[16]，壁面热流为

(2)

数值离散方法采用有限体积方法，对于粘性项的离散采用中心格式； 对于无粘项的离散，采用4种不同的数值格式，分别是二阶MUSCL格式、 五阶WENO格式、 MDCD格式以及MDADF-HY格式。在获得单元界面两侧变量的基础上，通过Roe格式近似求解黎曼问题得到单元界面处变量值。文献[20]指出该算例在此流动条件下为层流，因此无需考虑湍流模型对计算带来的影响。时间推进采用三级三阶Runge-Kutta方法。

作为中国生产整体硬质合金刀具的专业制造厂家，喜威一（北京）刀具有限公司拥有自己的品牌“CVE”，产品有硬质合金标准刀具（铣刀、钻头及铰刀）和非标准刀具（钻铰刀、内冷钻、阶梯钻、阶梯铰、复合刀具及成型刀具），并提供修磨、涂层服务；而刀具产品原材料均采用欧洲进口合金棒料，拥有世界一级的生产设备。另外，该公司还是意大利UFS丝锥的中国总代理。

1.2 模型选取及计算条件

三维球头钝锥是十分典型的高超声速飞行器再入体模型，本文的三维钝锥算例选自NASA TN D-5450报告，采用热电偶对不同钝度的钝锥模型测试了多种来流雷诺数及攻角情况下的热流值，试验数据丰富。其中驻点热流与理论相比测量误差小于±5%； 对于0°攻角，锥面上的三条母线上的热流测量重复性误差小于±6%； 由于精度是相对于热流测量水平来说的，因此对于低热流的测量误差估计在±20%以内。本文选择其中一种几何外形，如图1所示，钝锥半锥角=15°，球头曲率半径=0.009 5 m，锥尖到锥底的总长度=0.568 7 m。

图1 钝锥几何外形

计算状态： 来流马赫数为10.6，来流雷诺数=3.937×10，温度=47.34 K，壁面温度=294 K。选取攻角为10°以对比分析迎风面以及背风面流场信息。由于热流稳定要比压力缓慢得多，因此按照文献[22]给出的热流稳定步数在压力稳定步数的10倍以上作为收敛准则。

2 网格雷诺数对热流的影响分析

图2 整体计算网格

采用MDADF-HY格式计算无粘通量，壁面热流计算结果按无量纲形式给出，其中为计算得到的热流值，为文献[21]给出的驻点热流值，对统一作归一化处理。横坐标为至尖锥(非钝体头部)距离与钝锥总长度之比。取0.1，1，10，100，1 000五种不同网格雷诺数计算结果进行分析。图3分别给出了攻角为0°以及10°时五种不同网格雷诺数情况下钝锥表面热流沿轴向的分布，并且与试验值做出了对比。从图中可以得出，五种网格雷诺数计算得到的壁面热流大致有相同的趋势，头部驻点区由于受到气流滞止的作用以及正激波的影响，气动加热效果非常严重，导致热流密度非常高。离开驻点区域后，热流密度迅速下降，在相对位置大概为0.2之后，热流密度下降的趋势逐渐变缓，到达钝锥尾部时基本保持稳定。

图3 钝锥轴向热流分布曲线

随着网格雷诺数的不断增加，热流计算值呈现不断减小的趋势。对于10°攻角背风面，气动加热现象不明显，五种网格雷诺数对应的网格得到的热流计算值均能与试验值较好的吻合。但对于迎风面，网格雷诺数为1 000时得到的热流计算值与试验值差距比较大，尤其是在头部附近，最大误差已经超过了50%。可见网格尺度过大时，计算的热流值难以升上去，在钝锥大部分区域维持在较低的范围，远远达不到对气动热数值模拟的要求。

对比0°攻角以及10°攻角下的壁面热流值，在一定的攻角下，钝锥迎风面的热流密度更高，并且要远远大于背风面热流值，因此在热防护设计中，飞行器迎风面是需要重点考虑的部位。

3 格式精度对热流的影响分析

CFD计算中，数值格式是关系到数值模拟结果准确与否最主要的因素。对于高超声速气动热来说，热流值是与温度梯度相关的物理量，不同精度插值用到的模板点数量也不同，这就会导致计算得到的温度值及温度梯度值的精确度不同。本文选择三种不同精度的格式： 二阶MUSCL格式、 四阶MDADF-HY格式及五阶WENO格式，来研究格式精度对壁面热流计算的影响。

图4 不同精度格式驻点热流相对误差随网格雷诺数变化曲线

可以看出，当网格雷诺数足够低时，在0.1～10范围内，三种格式得到的壁面热流相对误差都能够收敛到某个确定值，精度越低得到的壁面热流相对误差越大，并且高阶格式计算出的热流值要低于低阶格式。比较四阶MDADF-HY格式和五阶WENO格式，当网格雷诺数较低时，WENO格式得到的热流值要比MDADF-HY格式得到的热流值更接近试验值。但是随着网格雷诺数的逐渐增大，到10的量级之后，MDADF-HY格式的计算值反而更接近试验值，分析认为MDADF-HY格式的分辨率比WENO格式更好，当网格比较稀时，MDADF-HY格式对于粘性边界层的模拟以及激波等间断结构的捕捉效果要更好。

同时，随着网格雷诺数由高到低，相比于四阶MDADF-HY格式，五阶WENO格式虽然收敛到最小相对误差时所需的网格雷诺数更小，但是达到收敛时计算结果更准确。分析认为，高阶插值求解壁面热流，在计算无粘通量时需要在控制体单元两侧采用更多的模板点，利用了更多的流场信息来插值计算界面处通量，因此得到的结果更加准确，对于壁面热流的模拟更逼近试验值。除此之外，二阶MUSCL格式由于精度比较低，对于热流计算的相对误差一直大于另外两种格式，并且MUSCL格式采用的单元两侧的模板点比较少，导致其对于网格雷诺数更加敏感，需要更密的网格才能收敛到相对误差最小值。

4 格式耗散对热流的影响分析

格式耗散的大小，可以通过分析数值解与解析解的振幅之比来定性说明。以前文提到的一维线性波动方程为例，线性波动方程的解析解与时间无关，不会随着时间的变化而变化，而数值解的幅值则会随时间变化。当数值解与解析解的振幅之比小于1时，数值解的振幅是随时间不断衰减的，则认为格式具有正的耗散，稳定性比较好； 当振幅之比为1时，数值解的幅值保持不变，则认为格式具有零耗散； 当振幅之比大于1时，数值解的振幅随时间不断增大，则认为格式具有负的耗散，稳定性比较差。

而数值格式的色散特性，则表征数值解与解析解之间相位的差别。对于色散误差的优化，需要用到目标函数，当目标函数取最小值时格式的色散特性达到最优。MDCD格式的优势在于调节耗散参数不会影响已经优化的色散特性，因此采用MDCD格式，通过选取不同的耗散参数得到一组精度、 色散相同，而耗散不同的通量计算方法，针对本文的高超声速气动热问题来研究耗散对热流计算的影响。文献[17]对色散误差进行了优化，色散参数与低波数误差占总误差的比重相关，本文取色散参数=0.04，对于耗散参数的选取，文献[17]提出为保证线性权重，需要满足≤的条件，因此取0.004, 0.012, 0.020, 0.028, 0.036五种不同工况进行数值计算。考虑到计算效率以及热流对于网格雷诺数和法向网格增长比的要求，计算网格的网格雷诺数取为10，法向网格增长比取为1.1。

图5为选取耗散参数为0.004, 0.012, 0.020, 0.028, 0.036五种情况下钝锥壁面热流分布曲线。对于攻角为0°以及攻角为10°背风面，不同耗散参数计算所得壁面热流结果大致相当，均能与试验值较好地吻合。在驻点附近，热流有一个急剧下降的过程，从驻点值过渡到下游的稳定值这一过程中，可以看出，耗散参数越大，壁面流计算值就会越大。

对于10°攻角下钝锥迎风面的数值模拟，随着耗散参数的增加，热流计算值大体也呈现出增加的趋势，尤其在过渡段更为明显。在钝锥后半部分，除=0.004，其余四种参数计算的壁面热流值几乎重合，与试验值的误差非常小。对于=0.004的情况，热流分布曲线比较曲折，热流值偏差也比较大，分析认为是迎风面气动加热现象比较严重，尤其是驻点附近温度梯度变化剧烈，过低的耗散会导致计算的稳定性变差，难以收敛到准确值。

在气动热计算过程中，数值格式的耗散也是十分重要的因素。对于离散求解N-S方程，过低的耗散会导致计算稳定性变差，使得计算难以收敛到稳定值。由于气动热的计算比气动力的计算更加困难，耗散比较低的情况下很多时候即使气动力已经收敛，气动热还未达到稳定值。当耗散过高时，热流计算值与试验值的误差也会随之增加，虽然此时的计算稳定性比较好，但是对于流场中微小尺度结构的分辨率会变得很差，同时会存在“抹平”间断的现象，这对于准确捕捉高超声速激波结构非常不利。所以合适的数值耗散对于热流密度的准确模拟是十分必要的。

图5 不同耗散参数情况下钝锥热流分布曲线

5 结 论

本文选取典型的高超声速飞行器再入体三维钝锥模型，按照不同网格雷诺数生成多组计算网格，计算了0°攻角以及10°攻角两种工况，通过分析迎风面、 背风面的壁面热流，研究了网格尺度、 格式精度、 格式耗散三种因素对热流计算的影响。

分析了网格雷诺数对热流计算的影响规律，由于10°攻角迎风面气动加热现象更为严重，因此五组网格雷诺数在此工况下对热流的计算值相差比较大，对于较大的网格雷诺数，热流计算误差更明显。从上文可以看出，当网格雷诺数小于等于1后，钝锥身部计算的壁面热流差别很小。当网格雷诺数小于等于10后，驻点计算的壁面热流差别很小。综合来看，当网格雷诺数小于等于1后，计算的热流才能达到网格无关性。

比较了MUSCL, MDADF-HY, WENO三种格式计算得到的热流计算结果，在网格雷诺数较低时，MUSCL误差大于10%，MDADF-HY与WENO维持在较低的误差水平，同时WENO计算更精确。当网格雷诺数增加到一定值时，MDADF-HY的热流计算反而比WENO更精确。

采用耗散可控格式进行了格式耗散对热流计算的研究，0°攻角以及10°攻角背风面的热流计算值相差比较小，与试验值均能达到较好的吻合。而对于10°攻角的迎风面，当耗散比较小时，头部区域热流值出现了波动，说明对于气动加热现象严重的部位，需要数值格式具有一定的耗散以维持计算的稳定，减少非物理振荡。