CN115544675B - 高超声速飞行器防热材料表面催化特性多尺度预测方法 - Google Patents

Abstract

本发明公开了一种高超声速飞行器防热材料表面催化特性多尺度预测方法，包括步骤：将高超声速来流远场条件输入CFD求解器中进行高超声速飞行器外流场的数值模拟；CFD计算收敛之后，通过后处理软件CFD‑Post提取壁面组分质量分数和温度分布；通过RMD求解器计算得到用于统计复合成分子的原子数量的产物文件；通过复合成分子的原子个数与撞击表面原子个数的比值计算表面催化复合速率，将其输入CFD求解器中进行迭代计算，得到壁面热流密度和流场温度分布。该方法针对高超声速飞行器防热材料表面催化复合特性进行多尺度预测，提高了对气固界面材料演化机理的表征能力和非烧蚀型热防护系统表面热预测的精度。

Description

高超声速飞行器防热材料表面催化特性多尺度预测方法

技术领域

本发明涉及航空宇航科学领域中对被动式热防护系统（非烧蚀型防热材料）进行流固耦合跨尺度数值模拟的研究，更具体地说，涉及一种高超声速飞行器防热材料表面催化特性多尺度预测方法。

背景技术

高超声速飞行器在临近空间以5马赫及以上速度飞行时，由于激波压缩和粘性阻滞而产生严重的气动加热现象。为了保证高超声速飞行器的飞行安全，必须克服飞行过程中产生的极端热载荷问题。热防护系统及其技术领域的发展已成为制约高超声速技术发展瓶颈的世界难题。热防护系统的可靠性在很大程度上依赖于对高超声速飞行器表面复杂气动热载荷的准确预测。在高温、高压、高焓的环境条件下，热防护系统表面会发生包括多相催化复合反应、烧蚀反应、热氧化反应和热裂解反应等复杂的物理和化学反应过程，这些表面复杂的物理化学反应，又进一步使得气动热载荷难以准确预测。

目前，高超声速飞行器的风洞实验研究能够实现简单的气动热载荷预测，但只能满足马赫数的要求，焓值和雷诺数等参数都无法在地面实验中复现临近空间高超声速环境，尤其是热防护系统表面发生的复杂多相催化复合反应，而数值模拟的方法则能够一定程度上弥补实验的缺陷。针对高超声速气动热数值模拟方法，中国发明专利CN108255781A公开了一种高超声速目标表面动态温度建模方法，将气动加热和材料结构瞬态热传导耦合起来，提高温度场仿真计算精度和效率；但该专利仅考虑了较为极端的完全催化壁面边界条件，然而实际流场环境下，壁面催化反应的催化重组系数会受材料、温度和来流撞击壁面组分影响，并不是该专利中所采用的完全催化壁。针对界面计算方法，中国发明专利CN115270534A公开了一种基于离散粒子-有限元网格耦合模型的界面计算方法及装置，对界面附近的吸附、扩散和传热等过程进行模拟，可用于研究材料催化特性；然而该专利仅局限于壁面研究，并未将其应用到高超声速气固界面热预测。Cinquegrana等人在论文(Cinquegrana D, Votta R, Purpura C, et al. Continuum breakdown and surfacecatalysis effectsin NASA arc jet testing at SCIROCCO[J]. Aerospace Scienceand Technology, 2019,88: 258-272.)对NASA的实验模型进行了深入分析，发现实验测量值和数值模拟预测值之间的重大差异是由于局部稀薄效应致使连续性假设失效，从而导致物理精度降低，并针对材料表面催化特性通过微观DSMC方法和宏观CFD方法进行多尺度模拟，以提高数值模拟精度，减小与实验测量值的差距；但该方法未能表征气固界面相互作用的规律，且无法将控制因素转化为气动热环境计算的边界模型，以提高宏观气动热环境的预报精度。因此，在高超声速非烧蚀型热防护系统领域，尚未形成一套完整的多尺度耦合数值模拟计算方法，既能够对气固界面非均相反应机理进行表征，又能够对宏观热防护系统进行高精度的热预测。

发明内容

本发明的目的在于提供一种高超声速飞行器防热材料表面催化特性多尺度预测方法，该方法耦合宏观尺度和微观尺度的数值模拟方法对气固界面间演化机理进行表征和表面热预测，提高对气固界面材料演化机理的表征能力和非烧蚀型热防护系统表面热预测的精度。

为实现上述目的，本发明提供如下技术方案：

一种高超声速飞行器防热材料表面催化特性多尺度预测方法，包括以下步骤：

S1，无壁面催化反应情况下，将高超声速来流远场条件输入CFD求解器中进行高超声速飞行器外流场的数值模拟；

S2，CFD计算收敛之后，通过后处理软件CFD-Post提取壁面组分质量分数和温度分布，并输入RMD求解器；

S3，对目标非烧蚀防热材料气固界面进行建模，得到目标非烧蚀材料的微观模型，通过RMD求解器对微观模型进行数值模拟，以CFD计算得到的壁面温度为加热温度，以CFD计算得到的壁面组分质量分数作为分子和原子的入射比例对微观模型中的固相原子模型进行撞击，得到用于统计复合成分子的原子数量的产物文件；

S4，待微观模型达到稳定时，通过复合成分子的原子个数与撞击表面原子个数的比值计算表面催化复合速率；

S5，将催化复合反应速率输入CFD求解器中进行迭代计算，得到考虑气固界面催化复合效应情况下高超声速飞行器壁面热流密度和流场温度分布。

进一步，所述步骤S1具体为：

将远场温度、压强、速度和远场组分质量分数的边界条件输入CFD求解器中，通过求解非平衡NSF方程和双温模型，得到结果文件，包括记录边界条件、湍流模型、组分输运模型、差分格式和网格的.cas文件和记录每个网格物理量数值的.dat文件。

进一步，所述CFD求解器采用ANSYS Fluent。

进一步，所述RMD求解器采用ReaxFF反应力场支持下的LAMMPS。

进一步，所述步骤S2具体为：

将.cas和.dat文件输入后处理软件CFD-Post中，通过计算功能，对壁面所有网格上存储的参数进行求和取平均，计算获得壁面温度T _w 、壁面N原子质量分数ω_N、壁面N₂分子质量分数ω_N2、壁面NO分子质量分数ω_NO、壁面O原子质量分数ω_O和壁面O₂分子质量分数ω_O2。

进一步，所述步骤S3中，通过Materials Studio和Packmol对目标非烧蚀防热材料气固界面进行微观建模，得到目标非烧蚀材料的微观模型，包括真空层、入射层、气固界面、加热层和固定层，加热层和固定层构成固相原子模型。

进一步，所述步骤S4包括以下子步骤：

S41，读取产物文件，所述产物文件中奇数行为输出的分子产物名称，偶数行为分子产物对应的分子数量，以字符串数组的形式导入MATLAB中；

S42，催化复合系数的计算公式为

遍历定位所有奇数行中N2字符串所在数组的位置，输出该奇数行对应的偶数行的时间步和N₂分子的数量信息，将N₂分子的数量乘以分子中单原子个数2再除以入射N原子的个数，得到催化复合系数γ随时间变化的趋势；

S43，选取催化复合系数稳定后的时间段，对该时间段内的催化复合系数取平均，得到宏观催化复合系数γ _h ，所对应的催化复合反应速率为

其中，k _w 为催化复合反应速率，γ _h 为宏观催化复合系数，M为气体分子的质量，R为气体常数，T _w 为壁面温度。

进一步，所述步骤S5具体为：

在CFD求解器中重新读取步骤S1中的.cas和.dat文件，在原有计算结果的基础上，加入壁面催化复合反应，反应速率与步骤S43计算得到的催化复合反应速率k _w 一致；CFD计算完成后，通过后处理软件CFD-Post获取壁面热流密度和外流场温度云图。

本发明与现有技术相比所具有的有益效果：

（1）本发明采用宏观和微观相耦合的多尺度预测方法，相比传统的单一宏观壁面热流分布预测手段精度更高，对界面处催化复合现象的表征更为准确；相比单一的微观分子动力学模拟，本发明弥补了其尺度过小、难以在宏观尺度应用的缺点。

（2）传统的高焓热环境与防热材料间界面气固非均相反应模型一般通过求解壁面附近质量守恒方程的方法构建，但现有的有限化学反应速率常数通常来源于采用Arrhenius公式对试验数据的线性拟合，过度简化的经验模型在精度和可靠性方面都有较大缺陷。而微观的RMD方法是通过牛顿运动方程计算模拟系统内每个粒子受到的力，再通过积分方法计算每个粒子的运动轨迹，最后通过对模拟系统内大量的粒子进行统计学平均的方法得到宏观流体性质，能够更准确地表征壁面催化反应，通过真实模拟得到壁面催化复合系数。

附图说明

图1是高超声速飞行器防热材料表面催化特性多尺度预测方法的流程图；

图2是气固界面原子分子入射原理示意图；

图3是非烧蚀材料的微观模型数据文件示意图；

图4是微观系统中分子统计示意图；

图5是催化复合系数随时间变化曲线示意图；

图6是壁面热流密度示意图；

图7是外流场温度云图示意图；

图8是建立壁面催化复合反应示意图。

具体实施方式

下面结合附图和实施例对本发明进行进一步的详细介绍。

一种高超声速飞行器防热材料表面催化特性多尺度预测方法，如图1所示，包括以下步骤：

S1：无壁面催化反应情况下，将高超声速来流远场条件输入宏观计算流体力学(CFD) 求解器中进行高超声速飞行器外流场的数值模拟。

具体地，将远场温度、压强、速度和远场组分质量分数的边界条件输入宏观CFD求解器中，通过求解非平衡NSF（Navier-Stokes-Fourier）方程和双温模型，得到结果文件，包括记录边界条件、湍流模型、组分输运模型、差分格式和网格的.cas文件和记录每个网格物理量数值如温度、压强、速度、组分质量分数等的.dat文件。

优选地，所述CFD求解器采用ANSYS Fluent。Fluent提供的网格特性十分灵活，应用范围非常广，可用于二维平面和三维流动分析，可完成定常与非定常流动分析、可压流与不可压流计算、层流和湍流模拟、传热和热混合分析、流固耦合传热分析和多相流分析等。除此之外，Fluent还能支持MPI多核并行运算，具有良好的前处理与后处理接口，前处理接口支持ICEMCFD、Pointwise等网格划分软件，后处理接口支持Tecplot、CFD-Post等后处理软件，以输出残差曲线和多种云图等。

其中，NSF方程通常用于计算连续介质假设成立的瞬态可压缩反应流，是宏观流场计算的基础。设高超声速流场包含N _s 种组分，组分以原子和分子的形式存在，N _s 种组分中包含N _m 种分子。在空间零通量散度形式下，NSF方程可以写成

其中，

和U为矩阵，能够表示为连续性方程、N _s 个组分输运和反应方程、动量方程、N _m 个振动能量方程和总能量方程。

为源项，t为时间。守恒量U可以写成

式中u，v，w是速度矢量在坐标系中三个方向的分量，ρ为总密度，ρ _s 为组分s的密度，E _ve,m 为分子m的振动电子能，E为总能量。通过该式可求得流场各位置处速度、密度和温度等参数。

源项

包含了质量源项和能量源项，用于满足质量和能量守恒的要求，可表示为

其中

由下式得到的

式中，Q _m,V-T 是通过与分子m反应的振动-平动能，Q _m,V-V 是通过与分子m反应的振动-振动能，Q _m,C-T 是通过与分子m反应而增加或减少的振动-电子能，Q _m,e-v 是通过与分子m反应的振动-电子能。该项为能量源项，目的是在于满足能量守恒，并用于求解得到温度等热力学参量。

是组分s的净质量产量，质量作用定律规定组分s的净来源是由化学反应引起的，计算表达式是由该组分参与反应的反应源项之和得到的，如下式所示

式中

和

是反应r中组分s的正向和反向化学计量系数，N _r 是反应r中组分的数量，M _s 为组分s的分子质量，M _k 和ρ _k 分别为第k个参与反应分子的分子质量及密度，k _f,r 为正向速率常数，k _b,r 为反向速率常数，

和

为第k个参与反应分子的正向和反向化学计量系数。该项为质量源项，目的是在于满足质量守恒，并用于求解得到密度和各组分质量分数等参数。

假定反应速率k _f 遵循阿伦尼乌斯定律：

其中A是指前因子，β是温度指数，T _a 是由活化能得出的活化温度，T _c,f 是正向反应的控制温度，通过该式可以计算得到各化学反应的反应速率，用于计算不同组分的质量分数。

与传统的能量方程不同的是，双温模型中的总能量可以被视为7种不同能量的总和，每种能量都由一个特定的温度表征，并相互交换能量，公式如下：

其中的第一项和最后一项能量表达式分别表示动能和化学能。其余各项E _t ,E _r ,E _v ,E _e 和E _el 分别代表了平动、转动、振动、电子和电能。在动能表达式中，ρ是总密度，u _i 表示三个速度分量，指数i从1到3的分别代表x，y，z三个方向。ρ _s 和

在化学能项中，指的是组分s的密度和标准生成焓。通过能量方程可以计算得到速度分量、组分密度和温度等参量。下式以转动能为例，说明了某一模态的总能量与组分s的单位质量能量之间的关系。

其中，e _r,s 为组分s单位质量转动能。

CFD求解器中采用的双温度Park模型假定平动温度等于转动温度，这是由于平动能和转动能之间能量交换过于迅速。该模型还假设分子的振动温度、电子的平动温度以及分子和原子的电子激发温度这三个温度是相同的。这一假设主要是基于分子的振动运动与自由电子的平动运动之间的高能量转移率，以及原子和分子的低能电子态与基态之间的高平衡速度。下面列出了每种模式下每单位质量物质的不同能量的表达式。

其中，e _t ,e _r ,e _v,s ,e _e ,e _el,s 分别为单位质量平动能、转动能、振动能、电子能和电能，T _tr 为平动-转动温度，T _ve 为振动-电子温度，R _s 为组分s的气体常数，θ _v,s 是组分s的特征振动温度，g _i,s 和θ _el,i,s 是组分s特定电子态i的简并度和特征电子温度。通过两个温度T _tr 和T _ve 可以对每种能量形式进行数值描述，代入总的能量方程后，即可通过迭代求解出两个温度T _tr 和T _ve 。简化过后的双温模型公式，如下所示：

。

S2：宏观CFD计算收敛之后，通过后处理软件CFD-Post提取壁面组分质量分数和温度分布，并输入微观Reactive Molecular Dynamics (RMD)求解器。

具体地，将.cas和.dat结果文件输入后处理软件CFD-Post中，通过计算功能，对壁面所有网格上存储的参数进行求和取平均，计算获得壁面温度T _w 、壁面N原子质量分数ω_N、壁面N₂分子质量分数ω_N2、壁面NO分子质量分数ω_NO、壁面O原子质量分数ω_O和壁面O₂分子质量分数ω_O2。

S3：通过Materials Studio和Packmol对目标非烧蚀防热材料气固界面进行微观建模，如图2所示，该模型为由C、H和O原子组成的酚醛树脂模型，沿纵向划分为真空层、入射层、气固界面、加热层和固定层，其中加热层和固定层共同构成固相原子模型，并输出模型数据文件，如图3所示。通过微观RMD求解器对微观模型进行数值模拟，以宏观CFD计算得到的壁面温度为加热温度对加热层进行等温加热，以宏观CFD计算得到的壁面组分质量分数作为分子和原子的入射比例在入射层射出，对固相原子模型进行撞击，得到用于观察壁面原子催化复合成分子过程的原子轨迹文件和用于统计复合成分子的原子数量的产物文件。

具体地，通过微观计算软件LAMMPS(Large-scale Atomic/Molecular MassivelyParallel Simulator)对微观模型进行区域划分，划分为真空层、入射层、气固界面、加热层和固定层；设定固相原子模型加热层加热温度为宏观CFD计算得到的壁面温度；设定入射层组分入射比例为宏观CFD计算得到的壁面组分质量分数比例；设定原子轨迹文件的输出时间间隔，并通过分子间距离判断原子是否形成分子产物，统计所有时刻的分子产物数量，得到不同时刻下的分子产物文件。

优选地，所述微观RMD求解器采用在ReaxFF (Reactive Force Field)反应力场支持下的LAMMPS。LAMMPS是一基于分子动力学的开源数值模拟软件，支持多种势函数，具有良好的扩展性，可用于模拟气态、液态或者固态相形态下各种系综高达百万级的原子分子体系。ReaxFF反应力场是基于键级的化学反应力场，为研究原子分子体系中可能发生的化学反应提供了可能，并且能够在不给定化学反应路径的情况下模拟可能发生的化学反应，大大提高了分子模拟的合理性与准确性。ReaxFF反应力场覆盖范围很广，能够很好的模拟大分子聚合物热裂解的过程。

S4：待微观模型达到稳定时，通过复合成分子的原子个数与撞击表面原子个数的比值计算表面催化复合速率。

具体地，包括以下子步骤：

S41，读取LAMMPS分子动力学模拟输出的“species.out”产物文件，如图4所示，该文件中奇数行为输出的分子产物名称，偶数行为分子产物对应的分子数量，以字符串数组的形式导入MATLAB中；

S42，催化复合系数的计算公式为

遍历定位所有奇数行中“N2”(即N₂)字符串所在数组的位置，输出该奇数行对应的偶数行(即下一行)的时间步和N₂分子的数量信息，将N₂分子的数量乘以分子中单原子个数2再除以入射N原子的个数，即可得到催化复合系数γ随时间变化的趋势，如图5所示；

S43，选取催化复合系数稳定后的时间段，即该时间段内催化复合系数变化幅度小于0.02，如图5中被灰色方框框选出来的200-300 ps区间，对该时间段内的催化复合系数取平均，得到步骤S5中所需的宏观催化复合系数γ _h ，所对应的催化复合反应速率为

其中，k _w 为催化复合反应速率，γ _h 为宏观催化复合系数，M为气体分子的质量，R为气体常数8.314 J/(mol·K)，T _w 为壁面温度。

S5：将微观方法计算得到的催化复合反应速率输入宏观CFD求解器中进行迭代计算，得到考虑气固界面催化复合效应情况下高超声速飞行器壁面热流密度和流场温度分布，分别如图6和图7所示。

具体地，在宏观CFD求解器中重新读取步骤S1中的.cas和.dat文件，在原有计算结果的基础上，加入壁面催化复合反应，如N + N → N₂等，反应速率与步骤S43计算得到的催化复合反应速率k _w 一致，如图8所示，在组分输运模型中建立N + N → N₂反应，并将反应类型定义为壁面反应 (Wall Surface),阿伦尼乌斯速率 (Arrhenius Rate)中的指前因子(Pre-Exponential Factor)即为催化复合反应速率k _w ，阿伦尼乌斯定律如下式所示

其中k _f 是反应速率，A是指前因子，β是温度指数，T _a 是由活化能得出的活化温度，T _c,f 是正向反应的控制温度，通过该式可以计算得到各化学反应的反应速率，用于计算不同组分的质量分数和反应热，从而计算得到温度和热流分布。由于步骤S43中通过催化重组系数计算得到的反应速率公式与ANSYS Fluent软件中通过阿伦尼乌斯定律求得的反应速率公式的表现形式不一致，故直接将步骤S43中的反应速率公式计算得到的催化复合反应速率k _w 作为指前因子A，温度指数β和活化温度T _a 均设为0，即此时满足

。

宏观CFD计算完成后，通过后处理软件CFD-Post获取壁面热流密度和外流场温度云图等所需参数。

以上所述仅为本发明的具体实施方式，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (8)

1.一种高超声速飞行器防热材料表面催化特性多尺度预测方法，其特征在于，包括以下步骤：

S1，无壁面催化反应情况下，将高超声速来流远场条件输入CFD求解器中进行高超声速飞行器外流场的数值模拟；

S2，CFD计算收敛之后，通过后处理软件CFD-Post提取壁面组分质量分数和温度分布，并输入RMD求解器；

S3，对目标非烧蚀防热材料气固界面进行建模，得到目标非烧蚀材料的微观模型，通过RMD求解器对微观模型进行数值模拟，以CFD计算得到的壁面温度为加热温度，以CFD计算得到的壁面组分质量分数作为分子和原子的入射比例对微观模型中的固相原子模型进行撞击，得到用于统计复合成分子的原子数量的产物文件；

S4，待微观模型达到稳定时，通过复合成分子的原子个数与撞击表面原子个数的比值计算表面催化复合速率；

S5，将催化复合反应速率输入CFD求解器中进行迭代计算，得到考虑气固界面催化复合效应情况下高超声速飞行器壁面热流密度和流场温度分布。

2.根据权利要求1所述的高超声速飞行器防热材料表面催化特性多尺度预测方法，其特征在于，所述CFD求解器采用ANSYS Fluent。

3.根据权利要求1所述的高超声速飞行器防热材料表面催化特性多尺度预测方法，其特征在于，所述RMD求解器采用ReaxFF反应力场支持下的LAMMPS。

4.根据权利要求1-3任一项所述的高超声速飞行器防热材料表面催化特性多尺度预测方法，其特征在于，所述步骤S1具体为：

将远场温度、压强、速度和远场组分质量分数的边界条件输入CFD求解器中，通过求解非平衡NSF方程和双温模型，得到结果文件，包括记录边界条件、湍流模型、组分输运模型、差分格式和网格的.cas文件和记录每个网格物理量数值的.dat文件。

5.根据权利要求4所述的高超声速飞行器防热材料表面催化特性多尺度预测方法，其特征在于，所述步骤S2具体为：

将.cas和.dat文件输入后处理软件CFD-Post中，通过计算功能，对壁面所有网格上存储的参数进行求和取平均，计算获得壁面温度T _w 、壁面N原子质量分数ω_N、壁面N₂分子质量分数ω_N2、壁面NO分子质量分数ω_NO、壁面O原子质量分数ω_O和壁面O₂分子质量分数ω_O2。

6.根据权利要求5所述的高超声速飞行器防热材料表面催化特性多尺度预测方法，其特征在于，所述步骤S3中，通过Materials Studio和Packmol对目标非烧蚀防热材料气固界面进行微观建模，得到目标非烧蚀材料的微观模型，包括真空层、入射层、气固界面、加热层和固定层，加热层和固定层构成固相原子模型。

7.根据权利要求6所述的高超声速飞行器防热材料表面催化特性多尺度预测方法，其特征在于，所述步骤S4包括以下子步骤：

S41，读取产物文件，所述产物文件中奇数行为输出的分子产物名称，偶数行为分子产物对应的分子数量，以字符串数组的形式导入MATLAB中；

S42，催化复合系数的计算公式为

遍历定位所有奇数行中N2字符串所在数组的位置，输出该奇数行对应的偶数行的时间步和N₂分子的数量信息，将N₂分子的数量乘以分子中单原子个数2再除以入射N原子的个数，得到催化复合系数γ随时间变化的趋势；

S43，选取催化复合系数稳定后的时间段，对该时间段内的催化复合系数取平均，得到宏观催化复合系数γ _h ，所对应的催化复合反应速率为

其中，k _w 为催化复合反应速率，γ _h 为宏观催化复合系数，M为气体分子的质量，R为气体常数，T _w 为壁面温度。

8.根据权利要求7所述的高超声速飞行器防热材料表面催化特性多尺度预测方法，其特征在于，所述步骤S5具体为：

在CFD求解器中重新读取步骤S1中的.cas和.dat文件，在原有计算结果的基础上，加入壁面催化复合反应，反应速率与步骤S43计算得到的催化复合反应速率k _w 一致；CFD计算完成后，通过后处理软件CFD-Post获取壁面热流密度和外流场温度云图。

Images