三维复杂外形高速飞行器流‑固‑热快速计算方法

摘要

本发明公开了一种三维复杂外形高速飞行器流‑固‑热快速计算方法，属飞行器气动计算领域。本发明依据边界层理论，采用CFD数值计算与工程算法相结合的方法，静气动热的计算简化为绕飞行器的无粘外流数值解和边界层内热流求解两份部分，同时耦合防热系统结构传热计算模型。本发明结合了全长数值模拟与工程算法各自的优势，可以快速高效对三维复杂外形高速飞行器的多种飞行状态下气动加热与结构耦合传热特性的计算方法，给出飞行器全机表面热流密度与防热结构温度场时变特性，弥补了直接数值模拟方法计算代价高、效率低、周期长等缺陷，同时拓展了工程算法的应用范围。

说明书

技术领域

本发明属飞行器气动计算技术领域，具体指代一种三维复杂外形高速飞行器流-固- 热快速计算方法。

背景技术

高超声速飞行器在大气层中持续飞行时，飞行器表面将承受巨大的气动热(王江峰, 伍贻兆.高超声速复杂气动问题数值方法研究进展[J].航空学报.2015,33(1):159-175)。气 动热会导致飞行器结构温度升高，改变结构内温度场分布，从而改变结构材料的物理属 性，而结构材料特性的改变会引起结构刚度和结构模态的变化，对飞行器的飞行安全造 成极大隐患。因此热防护问题对于高超声速飞行器设计显得尤为重要(程克明,吕英伟. 飞行器持续气动加热的耦合性分析[J].南京航空航天大学学报,2000,32(2):150-155)，越 来越受到世界各国的高度重视。

气动加热问题(卞荫贵,徐立功.气动热力学[M].中国科技大学出版社.1997:15-20)的 主要研究方法包括：数值模拟方法(SinhaK,K.Reddy D S.Effect of ChemicalReaction Rates on Aeroheating Predictions of Reentry Flows[J].Journal ofThermophysics and Heat Transfer,2011,25(1):21-33)、工程近似计算方法(Hamilton HH,WeilmuensterKJ, DeJarnette F R.Ap-proximate Method for Computing ConvectiveHeating on Hypersonic Vehicles Using Unstructured Grids[J].Journal ofspacecraft and rockets.2014,1288:1305)、 地面风洞实验(田旭昂,王成鹏,程克明.Ma5斜激波串动态特性实验研究[J].推进技 术.2014,35(8):1030-1039)及自由飞实验等，其中后两种方法因代价高昂等因素不适于工 程设计的初期与选型、改型阶段。

数值模拟方法主要是对Euler方程、N-S(Navier-Stokes)方程及相关简化形式的控制 方程进行求解，具有计算精度高、可处理复杂流动和全机外形等优点，但在气动热与结 构传热耦合问题的求解方面对计算资源需求巨大且非常耗时。国内外在求解高超声速气 动热相关问题的数值计算方法方面做了大量研究，主要工作集中在计算效率与精度等方 面。

气动加热工程计算方法(Zhao J S,Gu L X,Ma H Z.A rapid approach toconvective aeroheating prediction of hypersonic vehicles[J].Science ChinaTechnological Sciences,2013: 1-15)对简单外形的气动热求解具有高效、准确的特点，因此在实际工程应用中率先得 到发展。但在复杂气动外形的气动热分析方面适应性较差，并且需要基于大量实验数据 对计算方法与结果进行人工修正。在这方面比较著名的是NASA兰利研究中心研发的一 套气动加热预测程序(MINIVER)(李会萍.高超声速飞行器气动加热特性及其计算方法 研究[D].上海：上海交通大学,2010)，该程序中驻点区域使用经典的Fay-Riddle公式，采 用参考焓方法来计算高速流动压缩性效应，另外可对过渡流区气动热进行计算。 Hamilton(Hamilton H H,Weilmuenster KJ,DeJarnette FR.Approximate method for computing laminar and turbulent convective heatingon hypersonic vehicles using unstructured grids[C].The 41st Annual AIAAThermophysics Conference,AIAA Paper.2009, 4310:2009.)针对三维钝头体发展了一种适用于空气平衡气体的高超声速流动热流密度 计算方法，该方法可计算不同高速流动状态(层流、转捩、湍流)下的热流密度。Zoby(Zoby E V,Simmonds A L.Engineeringflowfield method with angle-of-attack applications[J]. Journal of spacecraftand rockets,1985,22(4):398-404)等人开发了LATCH方法，该方法 基于参考焓和修正雷诺比拟计算热流密度，可用于有化学反应参与的钝头体气动热计 算。李建林(李建林,唐乾刚,霍霖,程兴华.复杂外形高超声速飞行器气动热快速工程估算 [J].国防科技大学学报.2012,34(6):89-93)等人对气动热工程计算方法进行拓展应用，对乘 波体构型高速飞行器进行气动热快速估算，得到较好结果。总体而言，工程计算方法存 在需要大量人工干涉及以庞大实验数据为基础等缺陷，在处理复杂外形飞行器方面的通 用性仍需完善。

发明内容

针对于上述现有技术的不足，本发明的目的在于提供一种三维复杂外形高速飞行器 流-固-热快速计算方法，以解决直接数值模拟方法计算代价高、效率低、周期长等缺陷，同时拓展了工程算法的应用范围。

为达到上述目的，本发明采用的技术方案如下：

一种三维复杂外形高速飞行器流-固-热快速计算方法，其特征在于，包括步骤如下：

步骤1、无粘外流场数值解；

步骤2、工程表面热流率计算；

步骤3、工程结构传热耦合计算；

步骤4、弹道状态动态插值方法；

步骤5、高温化学非平衡效应。

步骤1具体包括：

控制方程采用三维欧拉方程对无粘外流场求解，采用块结构网格机翼基于分布式并 行计算技术，取无粘流场解结果中的物面参数作为边界层外缘参数提供给工程方法中的 边界层简化算法。输出的物面参数包括：坐标(m)、物面速度分量(m/s)、马赫数、压 强(Pa)、密度(kg/m3)、静温(K)、总能(J/m3)。而边界层的工程算法可以得到用于 防热系统中结构传热计算所需的飞行器表面热流及传热系数等参数。

步骤2具体包括：

本发明针对的是全机外形，根据热流密度工程计算方法将飞行器表面热流的计算划 分为驻点区和非驻点区两个区域。驻点区热流计算采用目前广泛使用的Fay-Riddell公式 (中国人民解放军总装备部军事训练教材编辑工作委员会.高超声速气动热和热防护[M].北京:国防工业出版社,2003:46-136.)：

式中：ρ_w、μ_w、h_w分别表示物面密度、物面粘性系数和物面焓，ρ_s、μ_s分别为驻点密>D为平均空气电离焓，计算中取普朗特数Pr＝0.71、路易斯数>

非驻点区的热流密度计算公式(卞荫贵,钟家康.高温边界层传热[M].科学出版社.1986.)，本发明采用的是工程中常用的平板传热模型。

步骤3具体包括：

热防护结构内表面采用绝热壁边界条件，依据防热结构材料的毕奥数B_i的大小，采>

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

当B_i＜0.1时，采用如下式所示的热薄壁传热模型：

初始条件为：

T_w|_t＝0＝T₀>

式中ρ_δ、c_δ、ε分别为表示分别表示防热材料的密度、比热容和表面辐射系数。采用差>

当B_i＞0.1时，采用热厚壁传热模型。将热厚壁由内而外分为j层进行逐层推进求解：

表层：

最内层：

中间层：

初始条件为：

T_j|_t＝0＝T₁|_t＝0＝T_n|_t＝0＝T₀>

式中下标m为为材料类型，n表示层数，λ_m为材料m的热传导系数。利用差分离散可计>

气动加热与结构传热过程的耦合特性(吴洁,阎超.气动热与热响应的耦合研究[J]. 导弹与航天运载技术,2009(004):35-39.)以物面温度T_w作为物面边界输入条件，然后通>n；而在计算结构传热时，又以热流密度>n作为热边界条件，计算得到物面温度。在时间步t_n的求解中，以上一时间步t_n-1的物>作为本时间步的物面温度边界条件，采用工程算法进行边界层的气动加热计 算，得出t_n时间步的热流密度然后以作为热边界输入条件提供给结构传热计算 模块，计算得到t_n时间步的物面温度如此循环以实现结构传热和气动热环境的耦合 计算。

步骤4具体包括：

计算非定常弹道状态的气动加热与结构传热耦合计算需要足够多的无粘外流的数值 解作为其输入参数。

本发明利用无粘外流解动态插值方法提高耦合算法计算效率：通过已经预先完成的 有限个无粘外流解的流场结果，插值得到当前弹道时间点上飞行条件下的流场解，以供气动加热与结构传热计算所需，具体步骤如下：

1)根据飞行器物面上的当地流场参数与来流参数的比值在不同飞行高度下几乎保持不 变(Bova S W,Howard M A.Coupling strategies for high-speed aeroheatingproblems[R]. Sandia National Laboratories,2011.)这一规律，对同一飞行器，已经某高度H₀下的>0以及该高度下的流动参数P_inf0，同时已知任意高度下的大气参数P_infx，可以通过P_x＝P_infx×P₀/P_inf0的简单换算获得该飞行器在任意高度上的边界层>

2)由1)，对于整个弹道状态，无粘外流解只需计算设定参考飞行高度下的两个马赫数 和两个迎角，即4个飞行状态，其他飞行状态则可根据插值方法得到。这四个计算 状态分别为设定飞行高度下的两个马赫数和两个迎角的组合，即：>1,α₁),(Ma₁,α₂),(Ma₂,α₁),(Ma₂,α₂)。与这四个状态所对应的飞行器表面任一点>11、P₁₂、P₂₁、P₂₂，则某个弹道时间点上飞行条件>x,α_x)状态下的飞行器表面同一位置处的流动参数P_xx可通过如下插值算法得到：

步骤5具体包括：

在实际高超声速飞行条件下，空气的高温非平衡效应对气动热影响显著(欧阳水吾, 谢中强.高温非平衡空气扰流[M].国防工业出版社,2003)。高温空气非平衡边界层驻点 传热与平衡边界层驻点传热表面热通量的关系式表示为：

式中，下标O、N分别表示氧原子和氮原子，φ为壁面催化因子；h⁰为生成焓；C_O.s和C_N.s分别表示氧原子和氮原子的质量浓度；Le为路易斯数，取值范围1～2。

对于组元参数特性，使用组元参数简化计算模型，在流场温度9000K以下时，只考虑O₂、O、O⁺、N₂、N、N⁺、NO、NO⁺、e^-等9种组元及以下6个化学反应方程式：

其中K_i为摩尔密度平衡常数。本发明暂不考虑NO、NO⁺、e^-、N⁺、O⁺。

氧原子和氮原子的壁面催化因子有如下关系成立：

其中，[O₂]₀，[N₂]₀分别表示为空气中氧分子和氮分子的摩尔密度。联立式(11)、(12)>

取一级近似φ_O＝φ_N＝0。联立反应方程①②③，可得到[O₂]、[N₂]和[NO]的摩尔密度计算>

联立反应方程④⑤⑥和电荷守恒方程(唐国庆.高超声速多体分离动态气动加热计算技 术[D].南京：南京航空航天大学,2013.)，可以得到电子的摩尔密度：

考虑高温化学非平衡效应下的气动加热计算方法可修正为：首先由Fay-Riddell驻点 热流密度计算公式(1)得到平衡条件下热流密度q_eq，然后根据化学反应计算模型得到各组元摩尔密度和质量分数，最后由驻点化学非平衡边界层与平衡边界层热流密度关系式(10)求得高温非平衡条件下的热流密度q_ne

本发明的有益效果：

本发明依据边界层理论，采用CFD(Computational Fluid Dynamics)数值计算与工 程算法相结合的方法，将气动热的计算简化为绕飞行器的无粘外流数值解和边界层内热 流求解两个部分，同时耦合防热系统结构传热计算模型，发展了可用于快速计算三维复杂外形高速飞行器气动加热与结构传热特性的计算方法，实现了复杂飞行条件下飞行器全机表面热流密度与防热结构温度场时变特性的计算。

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

附图说明

图1a为沿子午线的组元摩尔浓度分布(t＝1000s)。

图1b沿子午线表面温度分布(t＝1000s)。

图1c为表面温度：不考虑化学非平衡效应(t＝1000s)。

图1d为表面温度：考虑化学非平衡效应(t＝1000s)。

图2为计算模型及物面网格和TPS方案。

图3a为第1000s时得飞行器防热层外表面温度分布。

图3b为第1000s时得飞行器防热层内表面温度分布。

图3c为飞行器在计算初始时刻(t＝0s)的防护层外表面热流密度分布。

图3d为飞行器在计算结束时刻(t＝1000s)的防护层外表面热流密度分布。

图3e为驻点内/外层温度随时间的变化曲线。

图3f为驻点热流密度随时间的变化曲线。

图4为弹道状态动态插值方法计算流程图。

图5a为弹道状态t＝0s时的表面温度分布。

图5b为弹道状态t＝100s时的表面温度分布。

图5c为弹道状态t＝350s时的表面温度分布。

图5d为弹道状态t＝600s时的表面温度分布。

图5e为弹道状态t＝1000s时的表面温度分布。

图6为本发明计算技术总体流场示意图

具体实施方式

为了便于本领域技术人员的理解，下面结合实施例与附图对本发明作进一步的说明，实施方式提及的内容并非对本发明的限定。

1、RAMC-II巡航状态长时加热(高温非平衡效应)

1.1技术参数

计算外形采用典型高速飞行器RAMC-II球锥体试验模型，其几何参数为：头部曲率半径0.152m，半顶角9°，模型长度1.3m。无粘外流场解的计算网格为总单元数约 40万、物面单元约1.5万的块结构网格。设定的巡航状态为：Ma_∞＝12，迎角α＝0°，飞>

长时加热中考虑结构传热，飞行器热防护结构材料为金属Ti，厚度2mm，表面发 射率为0.8。结构传热计算中，时间步长取为0.05s(即总时间推进步数为2万步)，采用 热薄壁传热模型.该算例在普通微机上计算耗时约17minCPU机时。

2.2数值模拟结果

图1给出了该算例数值计算结果。图1a为九组元模型O、N、O₂、N₂、NO、e^-、>+、O⁺、NO⁺在长时加热第1000秒时刻沿子午线分布的组元摩尔浓度分布。由于在本>2为主，O2基本分解完毕，N₂部分分解，计算中O₂与N₂的分解特性与>

2、X-37B巡航状态及弹道状态长时传热

2.1、技术参数

本发明的计算模型选取类X-37B外形，如图2所示，图2同时给出了计算所用的表面网格。计算状态如表1所示。

表1

用于无粘外流求解的块结构网格总单元数约512万，物面网格单元约25.1万。在热防 护结构方面，该算例进行了更接近于工程设计的复杂方案，即飞行器头部区域设定为热防护区1，以TPS1表示，全机其他部位设为热防护区2，以TPS2表示(见图1)，不同的热 防护区域使用表2给出的不同的热防护方案，其中热防护结构材料的最外层(飞行器表面) 和最内层温度初值设为飞行高度上的大气环境温度245K。根据该算例的飞行条件，不考 虑高温非平衡效应影响特性。该算例耗时约28minCPU机时。

表2

2.2、数值模拟结果

图3a和3b分别为第1000秒时飞行器防热层外表面(TWsurf)和内表面(TWin)温 度分布，飞行器在计算初始时刻(第0s)和计算结束时(第1000s)的防护层外表面热流 密度分布如图3c、3d所示。由计算结果可见，巡航1000秒时，类X-37B热防护层内 表面温度在不同的热防护区域(TPS)出现明显差异(图3b)，而热防护层整个外表面 温度分布均匀。这是由于TPS1区使用的是热薄壁与隔热层的组合防热结构，并且隔热 层采用热沉性好的二氧化硅材料(SiO2)，故在该热防护区域内表面温度出现显著降低， 大部分区域温度在410K左右，说明防热层隔热效果良好，与工程实际情况相符。而TPS2 区仅采用了热薄壁结构，防热材料为0.002m的金属Ti，其热传导性很好，故防护层内 外表面温度很快达到平衡，与外表面温度差别很小，且防护区域热流密度很小，该区域 温度大部分在600K以上。图3e、3f分别为驻点内外层温度与热流密度随时间变化曲线， 可以发现TPS1内外表面温度变化差别较大，计算结果显示驻点区域防热层外表面温度 为最高温约为979K，内表面最高温度约为521K，造成此温度值差异的原因仍然是由于 在这个区域使用了SiO₂隔热层。

2.3、弹道状态长时加热

在以上算例的基础上，考虑弹道状态，即变高度、迎角、马赫数，弹道状态更加接近飞行器实际飞行包线状态。计算流程如图4所示。

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

表3

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

本发明具体应用途径很多，以上所述仅是本发明的优选实施方式，应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以作出若干改进，这些改进也应视为本发明的保护范围。