吕水燕，张传侠，叶 坤，徐 健

(1.中国兵器工业试验测试研究院， 陕西 华阴 714200； 2.西北工业大学航空学院， 西安 710072)

随着高超声速飞行器的进一步发展，严重的气动热是外形设计中首要解决的问题之一。关于气动热问题的研究方法主要有气动热的地面试验研究、数值计算研究和工程估算研究[1]。相对于地面风洞试验和实飞试验，数值计算运行成本低、设计周期短，可以解决风洞试验无法避免的壁面干扰和气体物性问题。数值计算主要是求解N-S方程各种近似形式即可得到热流分布，这种方法较为精确，在CFD计算热流的应用中最为广泛。美国NASA Langley研究中心开发了许多该类相关程序，最为成熟的有LAURA(Langlay Aerothermo-dynamic Upwind Relaxation Algorithm Code)、GASP(General Aerodynamic Simulation Program)、DPLR(Data-Parallel Line Relaxation Method for the Navier-Stokes Equations)等，这些软件涵盖了稀薄气体效应和高温气体效应、湍流等[2]。美国的这些软件已经发展应用得相当成熟，且目前气动热的数值计算已经开始向非结构网格技术发展[3-4]。我国气动热的研究起步较晚，数值计算主要通过求解NS方程预测高超声速气动加热问题，国内研究以NND格式系列化研究为代表，发展了一系列提高计算精度的方法[5]。关于气动热的数值计算，由于对近壁的第一层高度要求很高，因此，在采用的网格模式上，一般采用结构化网格，这就在一定程度上使得外形复杂或波系干扰强烈的问题很难划分高精度的结构网格进行数值计算求解[6-7]。闫超、潘沙等[8-9]在气动热的网格相关性上进行了相应的研究，而对非结构网格在气动热的分析上， 国内目前少有研究。

本文采用气动热的经典双椭球模型为研究对象，参考文献[10]中的试验条件，对研究对象在8.04Ma和10.02Ma的工况下进行气动热数值模拟研究。计算分别采用结构网格和非结构网格，对模型表面的Y+值进行评估，计算获得了双椭球体表面的气动压力分布和热流分布，并将两种网格模式下的模型中心子午线上的热流数据进行了对比分析。

1 计算模型

本文计算了文献[10]中所述的双椭球体标模的气动热特性，并与文献中的试验数据进行了对比分析。试验是在Ma=8.04，Re=1.13×107/m和Ma=10.02，Re=2.20×106m两种流场条件下进行的。本文计算针对0°攻角的计算结果进行对比分析，双椭球模型如图1所示。

图1 双椭球模型几何图

双椭球体给出几何公式如下，以保证高精度模型的生成。

试验按照以下两种状态进行，如表1所示。

表1 试验状态

2 边界条件推算

根据表1中的试验状态，推算试验仿真边界条件，包括：静压PS、静温TS和来流速度V。边界条件推算采用如下公式[11]：

总温计算公式：

马赫数计算公式：

萨特兰公式计算动力黏度：

萨特兰公式计算密度：

雷诺数计算公式：

利用以上公式推导出了计算所需的参数，获得边界条件汇总如表2所示。

表2 根据试验状态推导仿真边界

3 数值方法

结构网格划分采用ICEM 进行，计算采用两套网格进行计算，第1套网格单元数量为375万，近壁层网格高度为10-5m，层高按照1∶1.1的比例向外发散，第2套网格近壁网格高度值调整为10-6m。结构网格模型如图2(a)所示。

非结构网格同样采用ICEM进行划分两套网格，近壁采用棱柱层网格划分附面层，非结构网格的近壁网格高度同样划分为10-5m和10-6m，1∶1.1等比例向外发散，网格单元数554万。非结构网格模型如图2(b)所示。

图2 计算模型网格

针对复杂结构模型，采用结构化网格难度大，网格生成周期也很长，难以在工程上进行应用，这就需要研究结构网格和非结构网格计算结果的差异。本文同时采用非结构网格对双椭球体标模进行计算，与结构网格的计算结果进行对比分析，以期获得结构网格和非结构网格计算结果差异的规律。

SST湍流模型通过混合函数将标准k-ω模型和k-ε模型结合到一起，核心思想是在近壁面利用k-ω模型的鲁棒性，以捕捉到黏性底层的流动，而在主流区域则利用k-ε模型又可以避免k-ω模型对于入口湍动参数过于敏感的劣势，SST模型在工程上应用很广泛，本文计算采用SST湍流模型进行计算。

CFD计算格式的上风格式主要分为两大类：以Van Leer格式为代表的FVS(Flux Vector Splitting)格式和以Roe为代表的FDS(Flux Difference Splitting)格式。针对这两种格式各自的优缺点，又发展出AUSM[12-13](Advection Upstream Splitting Method)的混合格式。这种格式构造简单，无矩阵运算，激波分辨率高而且稳定性好，既具有FDS格式在边界层中解的精确性，也具有FVS格式在捕捉强间断时的健壮性，具有非常优秀的流动计算性能。经过10多年的发展，在高超声速的各种流动及气动热的数值计算研究中，获得了广泛的发展应用。一些商业软件如Fluent、Overflow等在高超声速的计算中也添加了该类格式。本文计算采用AUSM的数值离散格式。

4 不同网格模式对气动热的影响分析

4.1 Y+分析

热流的仿真计算中，近壁网格密度对热流计算结果影响很大，这已经成为共识。目前大多数讨论热流计算网格的文献把近壁面的Y+或者网格雷诺数作为衡量仿真精度的一个非常重要的无量纲值[8-9]。该值越低，计算精度越高，但是网格尺寸也越小，网格规模越大。本文参考Fluent内核代码中关于Y+的定义，为网格高度的特征尺寸与黏性尺度的比值，其公式为

其中:y是表征近壁网格高度的特征尺寸，对于结构网格来说y是第一层网格高度的一半，对于非结构网格来说，y第一层网格高度的1/3；μ是空气黏性系数；ρ是空气密度，τw壁面剪切应力。

本文从Y+入手，图3给出了本计算中双椭球体上下面中心子午线上的Y+的数值曲线。图3(a)为8.04Ma工况网格子午线上的Y+值，图3(b)为10.02Ma工况网格子午线上Y+值。由图可以看到，计算获得Y+数值均已控制在10以下，且结构网格和非结构网格在同等近壁网格尺度下获得的Y+数值基本一致，在第二椭球位置之后，非结构网格的Y+数值优于结构网格。10.02Ma工况Y+数值均达到了4以下，趋势规律与8.04Ma一致。

图3 双椭球体中心子午线的网格Y+值

4.2 计算结果分析

图4以Ma=8.04攻角为0°状态下为例，给出了流场静压试验与仿真的对比。从图4可以看到，流场的激波与试验纹影图吻合很好。

采用非结构网格，由于边界层外网格的稀疏性，导致压力云图在激波边界处的捕捉清晰度远远不如结构网格，压力云图的均匀性较差，流场激波区域显得模糊。但是，流场的高压区域的位置和结构网格的计算结果还是一致的。依照文献[4-5]的气动力的计算结果看，采用非结构网格，在气动静压方面仿真与试验数值吻合较好。而在气动热方面的差异性是本论文的重点。

壁面温度给定280 K，采用冷壁热流方法计算双椭球体中心子午线上的热流结果如图5所示。从计算结果来看，无论是结构网格还是非结构网格，仿真计算结果的热流分布趋势与实验数据均有较好的吻合度。从图5(a)、图5(c)可以看到，热流顺来流方向先是大幅度下降，然后趋于平缓，气流流经第2个小椭球体时形成第2道激波，热流在两椭球相贯线前明显有一个下降过程，两种仿真边界的结果同试验结果一致，这表明该位置处边界层发生了分离。小椭球处的2次激波的波后热流迅速升高，形成了第2个热流峰值，热流达到峰值后逐渐下降。下表面中心线上的热流分布如图5(b)、图5(d)所示，热流经过前缘高热流区域后逐渐下降，至模型中部趋于平缓。

图4 流场结构试验照片对比

图5 双椭球体上、下表面中心子午线热流

由于8.04Ma流场的压力和密度比10.02Ma的流场高出许多，参见表2，所以，8.04Ma的上表面热流较10.02Ma高出一些。

在8.04Ma上、下表面子午线热流曲线上，在前缘驻点位置，非结构网格的热流数值达到了1 000 kW/m2以上，与试验数据484 kW/m2偏差超过100%，而近壁层高0.01 mm的结构网格驻点位置的热流计算结果为557 kW/m2与试验结果的偏差缩小很多。这一点在10.02Ma的计算中同样得到印证。说明非结构网格在驻点位置的热流计算失真，这是由于在同样计算资源和网格规模的情况下，非结构网格在前缘驻点位置也就是在热流梯度很大，热流变化剧烈的位置，网格密度远远达不到结构网格的精细度所导致的。

在离开驻点位置之后至第二椭球边界层分离之前，在双椭球体的上子午线试验数据和结构网格与非结构网格的偏差量一致，两种网格在该位置有同等计算精度。在两椭球相贯位置，小椭球后的二次激波使得该处热流再次激增，形成第二个热流峰值，非结构网格在此处的模拟与试验数据更加接近，优于结构网格的模拟精度，这是由于该处非结构获得的Y+数值优于结构网格。随着热流缓慢下降，非结构网格比结构网格计算的热流趋势线与试验值吻合程度更好。在下子午线，试验数据与非结构网格偏差更小，与结构网格偏差稍有增大。

在除驻点热流梯度很大的位置外，计算模型表面大范围内非结构网格计算结果均优于结构网格，究其原因为本算例中，非结构网格的划分获得了比结构网格更优的Y+数值。这说明，在非结构网格近壁获得更佳Y+数值的情况下，可以获得比结构网格更好的模拟精度。

5 结论

1) 采用非结构网格计算对激波边界的捕捉不清晰，流场激波边界较为模糊。这是由于结构网格在外场主流区域的网格比非结构网格分布更为合理细密。

2) 本文计算热流Y+结构网格控制在10以下，非结构网格控制在3以下。非结构网格在获得比结构网格Y+更优数值的前提下，非结构网格在热流曲线上比结构网格与试验值有更高的吻合度，表面大范围内可以获得比结构网格更优的计算结果。

3) 非结构网格在驻点处热流计算，误差超过100%，这是因为在驻点处热流梯度很大，非结构网格在驻点处的分布很难做到像结构网格一样细密，在同等网格规模的情况下，非结构网格对高热流梯度的位置难以准确模拟。