高超声速全机外形气动加热与结构传热快速计算方法

2019-12-30李佳伟王江峰程克明伍贻兆

空气动力学学报 2019年6期

李佳伟, 王江峰, 程克明, 伍贻兆

(南京航空航天大学航空学院, 江苏南京 210016)

0 引言

高超声速飞行器在大气层中持续飞行时，飞行器表面将经受严峻的气动加热[1]。高超声速气动热会导致飞行器结构温度急剧升高，严重时局部高温还可能导致飞行器结构材料的物理属性和刚度发生改变，对飞行器的飞行安全造成极大隐患。因此快速预测气动加热与结构传热的物理过程，对高超声速飞行器设计的初步选型和热防护系统轻量化设计，具有重要作用[2]。

高超声速气动加热问题[3]的主要研究方法包括：数值模拟方法[4]、工程计算方法[5-6]、地面风洞试验[7]及飞行试验等，其中后两种方法因代价高昂等因素不适于工程设计的初期选型阶段。

数值模拟方法主要是对N-S(Navier-Stokes)方程及相关简化形式的控制方程进行求解，具有计算精度高、可处理复杂流动和全机外形等优点，但在气动热与结构传热耦合问题的求解方面对计算资源需求巨大且非常耗时。在这方面的研究，相关学者做了大量工作，主要工作集中在计算效率与精度等方面。董维中[8]等开展了高超声速飞行器气动加热与飞行器表面温度耦合的数值模拟研究，结果表明在气动热环境的预测中，应采用表面温度分布与气动热耦合计算的方法，以减小表面温度分布对气动热计算结果的影响。杨恺[9]等采用基于点隐式的二阶精度迎风TVD格式的N-S方程有限体积法，开展了高温热化学非平衡条件下球锥和RAM-CII模型的气动热数值模拟计算，分析了不同非平衡模型在热化学非平衡条件下流场的影响。计算结果表明建立的数值模拟方法具有较高的精度。夏刚[10]等依据流场特征时间步长远小于结构传热时间步长的特点，提出了一种松耦合的流-固-热耦合计算方法，并采用二维高超声速圆管绕流与结构传热的非定常算例对方法进行了验证。耿湘人[11]等基于Levelset方法和有限差分方法，将流场与固体结构统一到同一组控制方程中进行求解，计算结果与试验值吻合良好，探索了流-固-热一体化计算的可行性。张胜涛[12]等针对高超声速流场-热-结构松耦合分析策略框架，采用自适应耦合计算时间步长、混合插值策略等方法，建立了多场耦合分析平台。由于流-固-热物理问题复杂，该问题的求解在数值计算技术及计算设备硬件条件等方面要求甚高，采用全流场黏性数值模拟方法难以快速得到工程设计初期选型所需的参考数据。

气动加热工程计算方法[13]对简单外形的气动热求解具有高效、准确的特点，因此在实际工程应用中率先得到发展。但在复杂气动外形的气动热分析方面适应性较差，并且需要基于大量试验数据对计算方法与结果进行人工修正。在这方面比较著名的是NASA兰利研究中心研发的一套气动加热预测程序(MINIVER)[14]，该程序中驻点区域使用经典的Fay-Riddle公式，采用参考焓方法来计算高速流动压缩性效应，另外可对过渡流区气动热进行计算。Hamilton[15]针对三维钝头体发展了一种适用于空气平衡气体的高超声速流动热流密度计算方法，该方法可计算不同高速流动状态(层流、转捩、湍流)下的热流密度。Zoby[16]等人开发了LATCH方法，该方法基于参考焓和修正雷诺比拟计算热流密度，可用于有化学反应参与的钝头体气动热计算。乐嘉陵[17]针对高超声速流动的黏性效应，利用边界层近似方法和工程方法对高超声速飞行器的传热系数和表面摩擦阻力进行了相应计算。李建林[18]等人对气动热工程计算方法进行拓展应用，对乘波体构型高速飞行器进行气动热快速估算，得到较好结果。总体而言，工程计算方法的缺陷在于需要大量人工干涉以及庞大试验数据，在处理复杂外形飞行器方面的通用性上仍需完善。

综上，本文依据高超声速边界层[19]理论，利用数值计算方法与工程算法的各自优势，将气动热的计算简化为绕飞行器的无黏外流(边界层以外)数值解和边界层内热流求解两个部分，同时耦合防热系统结构传热计算模型，考虑了高温化学非平衡效应，发展了一种可用于快速计算三维复杂外形飞行器气动加热与结构传热特性的方法，实现了复杂飞行条件下飞行器全机表面热流密度与防热结构温度场时变特性的快速计算。

该方法的优点在于：综合了全流场数值模拟与工程算法各自的优势，可以快速、高效地对三维复杂外形高速飞行器的多种飞行状态下的气动加热与结构耦合传热特性进行计算与分析，给出热流密度与防热结构层内温度等参数分布，弥补了全流场黏性数值模拟方法计算代价高、效率低、周期长等缺陷，同时拓展了气动热工程算法的应用范围。

1 无黏外流场数值计算方法

气动加热与结构传热耦合数值计算技术考虑了瞬态气动加热、巡航状态长时传热和弹道状态长时传热三种模型。根据不同飞行状态完成无黏流场的求解，取无黏流场解结果中的物面参数(如表1所示)作为边界层外缘参数提供给工程方法中的边界层简化算法作为输入条件，之后可得到用于防热系统中结构传热计算所需的飞行器表面热流及传热系数等参数。

表1 无黏外流解物面输出参数

在无黏外流场求解方面，本文采用基于块结构网格技术及分布式并行计算技术的三维流场数值计算方法，为气动热工程估算方法提供如表1所列的物面流动参数。

2 物面热流与结构传热工程估算方法

2.1 物面热流密度计算

物面热流的计算采用工程中简化求解边界层方程的方法。将该工程方法所需的边界层外缘参数设定为无黏流场解的物面流动参数(见表1)，计算输出结果为物面热流密度。

本文研究针对的是全机外形，因此根据热流密度工程计算方法将飞行器表面热流的计算划分为驻点区和非驻点区两个区域。驻点区热流计算采用目前广泛使用的Fay-Riddell公式[20]：

(1)

式中ρw、μw、hw分别表示物面上气体的密度、黏性系数和焓，ρs和μs分别表示驻点处的气体密度及气体黏性系数，hD为平均空气电离焓。计算中取普朗特数Pr=0.71，路易斯数Le=1.0。非驻点区的热流密度计算公式参见文献[21]。

2.2 结构传热耦合计算方法

采用工程中常用的平板传热模型进行结构传热计算。传热方程为一维热传导方程，热防护结构内表面采用绝热壁边界条件，采用差分法进行求解。在传热计算模型的选取上，根据式(2)定义的防热结构材料的毕奥数Bi的取值来选取[22]。

(2)

式中α为传热系数，δ为材料结构厚度，λδ为当前材料热传导系数。

当Bi≤0.1时，采用式(3)所示的热薄壁传热模型：

(3)

初始条件为：

Tw|t=0=T0

(4)

式中ρδ、cδ和ε分别表示防热材料的密度、比热容和表面辐射系数。采用差分方法对式(3)进行时间t的推进求解，即可得到热薄壁条件下热防护结构层温度随时间的变化。

当Bi>0.1时，采用热厚壁传热模型。在具体计算时，将热厚壁由内向外分为j层进行逐层推进求解。计算方程为：

表层：

(5)

最内层：

(6)

中间层：

(7)

初始条件为：

Tj|t=0=T1|t=0=Tn|t=0=T0

(8)

式中下标m为材料类型，n表示结构材料的分层数，λm为材料m的热传导系数，该计算方法可用于计算多种材料组成的热防护结构的传热情况。将式(5)～式(7)进行差分离散并按时间步长推进求解，可以计算得出各层材料温度随时间的变化。气动加热与结构传热过程的耦合计算[23-27]比较复杂，在求解气动热参数时，需要以物面温度Tw作为物面边界输入条件，然后通过简化边界层工程方法计算得到表面热流密度qn；而在计算结构传热时，又以表面热流密度qn作为热边界输入条件，计算得到壁面温度。因此，本文采用图1所示的耦合迭代计算思路[23]。

图1 气动加热与结构传热耦合求解示意图

3 弹道状态动态插值方法

文中所述的弹道状态是指飞行器在真实飞行中随飞行时间改变飞行高度、迎角、马赫数三个参数的状态(暂不考虑侧滑角)。

根据第2小节中所述的物面热流与结构传热工程计算方法的特点，将其扩展到随时间变化非定常弹道(飞行包线)状态下的气动加热与结构传热耦合计算，此时需要足够多的无黏外流的数值解作为输入参数。为了在满足工程设计误差要求前提下提高耦合算法的计算效率，提出了一种无黏外流解动态插值方法：即通过已经预先完成的有限个无黏外流解的流场结果，插值得到当前弹道时间点上飞行条件下的流场解，以供气动加热与结构传热耦合计算使用。

由于高超声速飞行器飞行高度跨度大，通常涵盖整个大气层高度，若针对飞行高度H进行插值，必须要先构建飞行器在各个飞行高度条件下的流场数据库，其计算量非常庞大，难以达到快速计算的目的。针对该问题，参考国内外有关大气参数随大气高度的变化特征及拟合方法，采用如下方法对飞行高度参数的插值进行简化。

首先，在飞行高度H对流场参数的影响规律方面，研究发现飞行器物面上的当地流场参数与来流参数的比值在不同飞行高度下几乎保持不变[28]。根据这一规律，对同一飞行器，若已知某高度H0下的边界层外缘参数(以P0表示)以及该高度下的流动参数Pinf,0，且已知任意高度下的大气参数Pinf,x，那么就可以通过Px=Pinf,xP0/Pinf,0计算获得该飞行器在任意高度上的边界层外缘参数[24,28]。

由此一来，无黏外流解的计算可仅考虑马赫数和迎角两个参数。具体方法为：在对整个弹道状态进行计算时，无黏外流的求解只需计算设定参考飞行高度下的两个马赫数和两个迎角，共四个飞行状态，其他飞行状态可根据插值方法快速得到。这四个计算状态分别为设定飞行高度下的两个马赫数和两个迎角的组合，即：(Ma1,α1),(Ma1,α2),(Ma2,α1),(Ma2,α2)。与这四个状态所对应的飞行器表面上任一点处的流动参数的值分别记作P11、P12、P21和P22，则某个弹道时间点上飞行条件(Max,αx)状态下的飞行器表面同一位置处的流动参数Pxx可通过如下插值算法得到：

(9)

构造该插值算法的目的在于减少弹道状态无黏外流解的计算次数以提高数值仿真效率，该算法具有简单可行、易实现等特点。由式(9)表示的插值算法可看出，在流场参数的插值精度方面主要依赖于已经求得的无黏外流解的精度及插值变量间距(即马赫数与迎角的增量值)。因此在实际计算中根据弹道参数特征合理建立插值状态数据库及采用高精度无黏外流求解器将有利于提高流场参数的插值精度，能够最大程度地满足工程设计中的误差要求。该方法在插值精度与误差分析方面的详细分析参见文献[23]。

4 高温化学非平衡效应估算

计算模型中，将无黏外流解得到的物面流动参数作为边界层外缘参数，即该无黏流的物面流动参数不是实际黏性流场中飞行器的物面参数。由前述可知，实际物面的热流密度是由Fay-Riddell热流密度计算公式(1)得到的[29]。需要说明的是，该公式是在平衡边界层的前提条件下推导得出的，因而未考虑高温条件下的空气化学非平衡效应。在实际高超声速飞行条件下，空气的高温非平衡效应对气动热影响显著，因此对于实际高温化学非平衡流动，需要对Fay-Riddell驻点热流密度计算公式进行修正。式(10)给出了高温空气非平衡边界层驻点传热与平衡边界层驻点传热表面热通量的关系式[30]：

(10)

式中，下标O、N分别表示氧原子和氮原子，φ为壁面催化因子，h0为生成焓，CO,s和CN,s分别为氧原子和氮原子的质量浓度，Le为路易斯数。

通过式(10)可看出，求解高温化学非平衡流状态下的驻点热流，需要确定驻点位置氧原子和氮原子的质量浓度。本文主要发展的是一种面向工程应用的快速气动加热计算技术，因此在组元参数特性的获取方面，使用目前国内外广泛认可的组元参数简化计算模型，而不使用耗费资源巨大的精确多组元N-S方程数值解。依据Dunn和Kang[31]发展的各组元浓度计算模型，在流场温度9000 K以下时，可以只考虑O2、O、O+、N2、N、N+、NO、NO+、e-等九种组元及以下六个化学反应式：

其中Ki为摩尔密度平衡常数，对于不同的化学反应式，取值可由文献[31]附表查得。

由高温空气化学反应特性可知：当高温空气在9000 K以下时化学反应以离解反应为主，电离反应可忽略，在混合气体组成中，NO、NO+、e-、N+、O+的浓度比O2、O、N2、N的浓度小得多，即：

(11)

(12)

因此，在化学反应效应对气动加热的影响特性分析方面，暂不考虑NO、NO+、N+、O+、e-等5种组元。同时，氧原子和氮原子的壁面催化因子有如下关系成立：

(13)

(14)

其中[O2]0、[N2]0分别表示为空气中氧分子和氮分子的摩尔密度。联立式(13)、式(14)与反应方程①②求解，可以分别得到氧原子和氮原子的摩尔密度为：

(15)

(16)

由于φO、φN远小于1，因此在一级近似时可以取φO=φN=0。然后联立反应方程①～③，可得到[O2]、[N2]和[NO]的摩尔密度计算式为：

(17)

最后联立反应方程④～⑥和电荷守恒方程[19]可以得到电子的摩尔密度为：

(18)

在气动加热计算中，各组元质量分数需根据外流流动特征采用时间历程相关的迭代方法计算获得。由此，考虑高温化学非平衡效应下的气动加热计算方法可修正为：首先由Fay-Riddell驻点热流密度计算式(1)，得到平衡条件下热流密度qeq，然后根据化学反应计算模型得到各组元摩尔密度和质量分数，最后由驻点化学非平衡边界层与平衡边界层热流密度关系式(10)，求得高温非平衡条件下的热流密度qne。

研究中将如上所述的高温空气化学非平衡多组元效应嵌入到所发展的全机外形飞行器弹道飞行状态下的气动加热与结构传热数值计算方法中，使得计算模型更接近于高超声速飞行器真实飞行条件下的热环境。

5 算例与分析

根据以上理论分析与数值建模过程，发展了一种复杂飞行器气动加热与结构传热耦合计算技术。该技术主要由无黏外流解、气动热工程计算和结构传热三部分构成，如图2所示。

图2 高超声速气动加热方法示意图

在实际计算中，无黏外流解模块需根据计算任务要求完成所需无黏流动状态的计算，并形成后续计算所需的外流解数据库。边界层加热和结构传热的计算可根据飞行器实际飞行时间进行耦合求解。图3给出了计算技术总体流程示意图。

图3 求解流程示意图

下面采用三类典型的高超声速气动加热与结构传热算例来对所发展的耦合计算方法的效率、功能及精度进行具体分析。

5.1 典型外形气动加热验证

5.1.1 钝头体瞬态气动加热

本算例计算模型(图4)选取NASA报告[32]中的三维钝锥模型，模型长为0.352 m，头部半径0.027 94 m，锥角30°。无黏外流求解所用的计算网格为块结构网格，网格总单元数约32万，物面网格单元约1.1万。来流条件设定为：马赫数Ma∞=10.6，迎角α=20°，压强P∞=132.1 Pa，温度T∞=47.3 K；设定钝头体壁温Tw=294.44 K。

该算例在普通微机(3.3 GHz/4 GB，下同)上约在8 s CPU机时(不包含无黏外流解计算时间，下同)内即可完成瞬态气动加热的计算并输出飞行器表面热流分布。

图4 计算模型

图5为该算例物面热流密度计算结果。图5(a)为模型表面总热流密度分布云图。可见在有迎角情况下，计算所得的沿子午面热流密度分布具有很好的对称性，而且热流密度等值线图分布均匀，没有出现明显的锯齿、跳跃等现象。图5(b)为采用驻点热流值进行归一化处理后的子午面(截面方位角0°与180°)上热流密度分布与试验值[32]的比较，可以看到计算结果与试验值相符。由于飞行器大面积气动热数据的试验值难以获得，因此在热流具体数值的对比方面，取驻点作为对比验证点。表2给出了当前流动状态下驻点热流密度的计算结果与参考试验值对比，相对误差小于10%，基本满足气动热工程设计的误差要求。

(a)表面热流密度分布

(b)子午面热流对比

表2 驻点热流值对比

5.1.2 RAM-CII巡航状态长时加热

计算外形采用典型高速飞行器RAM-CII球锥体试验模型，几何参数为：头部曲率半径0.152 m，半顶角9°，模型长度1.3 m。无黏外流场解的计算网格为块结构网格，总单元数约40万，物面单元约1.5万。设定的巡航状态为：飞行高度H=60 km，迎角α=0°，Ma∞=6.0，初始壁温Tw=247 K，飞行时间t=1000 s。长时加热中考虑结构传热，飞行器热防护结构材料为金属Ti，厚度2 mm，表面发射率为0.8。结构传热计算中，时间步长取为0.05 s(即总时间推进步数为2万步)，采用热薄壁传热模型。该算例在普通微机上计算耗时约1020 s(17 min)CPU机时。

(a)外表面温度分布

(c)子午线温度分布对比

通过图6(c)中的对比可以看出，在子午线温度分布上，计算结果与辐射平衡温度吻合较好，当时间足够长时，表面温度计算结果愈趋近于辐射平衡温度，与理论分析一致。

5.2 复杂飞行器巡航状态长时传热

计算模型为如图7所示的类X-37B外形，图中同时给出了计算所用的表面网格。计算状态为巡航状态长时传热：Ma∞=5，迎角α=3°，飞行高度H=40 km，飞行时间t=1000 s。用于无黏外流求解的块结构网格总单元数约512万，物面网格单元约25.1万。在热防护结构方面，该算例进行了更接近于工程设计的复杂方案，即飞行器头部区域设定为热防护区1，以 TPS1表示，全机其他部位设为热防护区2，以TPS2表示(见图7)，不同的热防护区域使用表3给出的不同的热防护方案，其中热防护结构材料的最外层(飞行器表面)和最内层温度初值设为飞行高度上的大气环境温度245 K。根据该算例的飞行条件，不考虑高温非平衡效应影响特性。该算例耗时约1680 s(28 min)CPU机时。

图7 计算模型与TPS方案

表3 热防护系统参数设置

图8给出了计算结果，其中图8(a)和图8(b)分别为第1000 s时飞行器防热层外表面(Tw,surf)和内表面(Tw,in)温度分布，飞行器在计算初始时刻(第0 s)和计算结束时(第1000 s)的防护层外表面热流密度分布如图8(c)、图8(d)所示。由计算结果可见，巡航1000 s时，类X-37B热防护层内表面温度在不同的热防护区域出现明显差异(图8(b))，而热防护层整个外表面温度分布均匀。这是由于TPS1区使用的是热薄壁与隔热层的组合防热结构，并且隔热层采用热沉性好的二氧化硅(SiO2)材料，故在该热防护区域内表面温度出现显著降低，大部分区域温度在410 K左右，说明防热层隔热效果良好。而TPS2区仅采用了热薄壁结构，防热材料为厚度0.002 m的金属Ti，其热传导性很好，故防护层内外表面温度很快达到平衡，与外表面温度差别很小，且防护区域热流密度很小，该区域温度大部分在600 K以上。图8(e)、8(f)分别为驻点内外层温度与热流密度随时间变化曲线，可以发现TPS1区域内外表面温度变化差别较大，计算结果显示驻点区域防热层外表面最高温度约为979 K，内表面最高温度约为521 K，造成此温度值差异的原因仍然是由于在这个区域使用了SiO2隔热层。算例计算与分析结果表明，发展的复杂飞行器气动加热与结构传热耦合计算方法可用于复杂外形及多种热防护方案问题的求解。

(a)第1000 s 时外表面温度

(b)第1000 s 时内表面温度

(c)初始时刻外表面热流

(d)第1000 s外表面热流

(e)驻点内/外层温度

(f)驻点热流

5.3 复杂飞行器弹道状态长时加热

文中的弹道状态长时加热指的是飞行器按照给定的弹道状态(变高度、迎角、马赫数)飞行时的气动加热与结构传热特性，弹道计算状态更接近飞行器实际飞行包线状态。计算模型采用算例5.2的模型及计算网格。由于考虑的是弹道状态，因此使用了前文中的弹道状态动态插值方法。具体计算方法流程如图9所示。

由于缺乏相关弹道参数数据，本文根据相关文献设定的弹道参数如表4所示。弹道参数考虑了马赫数、飞行高速及迎角的变化，弹道飞行时间1000 s，在此期间飞行高度由60 km 降至30 km，飞行马赫数由Ma=6降至Ma=4，迎角变化范围为±5°。该算例全机表面采用同一种组合热防护方案模型，具体参数与算例5.2中 TPS1区相同。该算例耗时约1860 s(31 min)CPU机时。

表4 弹道参数

图10给出了该算例计算结果。从图10中可以清晰地看到飞行器表面沿弹道飞行时逐渐加热到平衡状态的过程，并且随迎角从-5°至5°变化过程中,飞行器头部与机翼前缘处的高温区往下表面移动，第1000 s时，最高温出现在飞行器头部下表面，约为1050 K，这与高超声速飞行器防热结构材料特性理论分析相符。由于没有在公开发表文献中找到与此类似的算例进行对比，因此这里仅给出本文算例计算结果的分析。该算例表明所发展的计算方法可对复杂高超声速飞行器在弹道飞行状态下的气动加热与结构传热耦合进行快速高效计算与分析，可为高速飞行器的热防护设计与初期选型提供技术参考。