CN106682392A - 复杂高超声速飞行器烧蚀效应快速计算技术 - Google Patents

Abstract

本发明公开了一种复杂高超声速飞行器烧蚀效应快速计算技术，属飞行器气动设计领域。该技术采用了基于碳基材料的烧蚀模型和基于温度分区的烧蚀效应计算方法，结合普朗特理论的无粘外流解与边界层理论，发展了一种碳基材料烧蚀效应快速数值计算技术，给出了复杂外形高超声速飞行器热防护系统表面烧蚀的热流密度分布，烧蚀率以及热防护结构温度分布的时变特性。本发明同时增加了高超声速飞行器在高速飞行弹道过程中，热防护系统表面发生烧蚀情况的模拟能力。该方法计算效率高，鲁棒性好，结果准确，填补了国内在高超声速飞行器发生烧蚀情况下，高效、快速、准确的给出数值模拟结果方面的空白，为高超声速飞行器初期设计阶段提供有效的技术支持。

Description

复杂高超声速飞行器烧蚀效应快速计算技术

技术领域：

本发明涉及一种复杂高超声速飞行器烧蚀效应快速计算技术，属飞行器气动设计领域。

背景技术：

随着高超声速飞行器的迅速发展，飞行速度越来越快，飞行条件越来越严苛，碳-碳复合材料逐渐取代了早期防护材料。热防护系统方案设计的关键问题之一是防热材料发生烧蚀的情况下，对气动加热的快速准确预测。

目前，对于高超声速飞行器的烧蚀防热问题，工程方法、数值方法和工程与数值相结合的方法的发展已相对成熟。实验是研究空气动力学问题最早的手段，各种研究方法所得的结果，一般以实验结果为准，其地位无可取代。目前，国内外仍是以实验方法为主要研究手段，因此，发展数值模拟计算预测变得十分重要。

本发明采用了基于碳基材料的烧蚀模型和基于温度分区的烧蚀效应计算方法，结合普朗特理论的无粘外流解与边界层理论，发展了一种碳基材料烧蚀效应快速数值计算技术，给出了复杂外形高超声速飞行器热防护系统表面烧蚀的热流密度分布，烧蚀率以及热防护结构温度分布的时变特性。

该方法增加了高超声速飞行器在高速飞行过程中，热防护系统表面发生烧蚀情况的模拟能力。该方法计算效率高，鲁棒性好，结果准确，填补了国内在高超声速飞行器发生烧蚀情况下，高效、快速、准确的给出数值模拟结果方面的空白，在高超声速飞行器初期设计阶段，提供有效的技术支持。

发明内容：

本发明主要采用数值计算手段，围绕高超声速气动加热和烧蚀效应，发展了一种复杂高超声速飞行器烧蚀效应快速计算技术。

本发明采用如下技术方案：一种复杂高超声速飞行器烧蚀效应快速计算技术，包括如下步骤：

步骤1、高超声速飞行器气动加热；

步骤2、材料的烧蚀模型；

步骤3、烧蚀效应快速计算；

步骤4、利用动态插值技术，对弹道状态进行准动态模拟。

进一步地，步骤1具体包括：

步骤1-1、无粘外流，采用欧拉方程作为控制方程求解无粘流动，取物面参数近似作为边界层的外缘参数，其中，所需要的边界层外缘参数有：网格信息、速度分量、马赫数、静压、密度、静温以及总能；

步骤1-2、高速边界层传热，计算表面热流密度需分为两个部分：驻点区和非驻点区，在驻点区，在任何情况下，驻点的流动都存在精确相似解，采用Fay-Riddell方法，驻点热流密度计算公式表示如下：

式中，Pr＝0.71，Le＝1.0；式中ρ_w为物面密度，μ_w为物面粘性系数,h_w为物面焓值，ρ_s为驻点处密度，μ_s为驻点处粘性系数，h_D为平均空气电离焓。

步骤1-3、气动加热流固耦合计算，对于结构传热，根据热防护结构的毕奥数Bi大小，防热层可以分为热薄壁和热厚壁两种类型；

对于热薄壁，采用一维热传导方程，计算公式如下：

初始条件：T_w|_t＝0＝T₀

式中，ρ_δ为材料密度；c_δ为材料比热容；δ为材料厚度；α为传热系数；ε为材料辐射系数；

对上式构建差分方程进行求解，并使用代替如下：

对于热厚壁，在计算中，由内向外热厚壁分为j层，简化后的热传导方程为：

(1)表层：

(2)最内层：

(3)中间层：

初始条件：T_j|_t＝0＝T₁|_t＝0＝T_n|_t＝0＝T₀

其中下标m为材料类型，n为第n层，λ_m为热传导系数；

同样构建差分方程，并做近似处理为的近似代替，为的近似代替，如下：

(1)表层：

(2)最内层：

(3)中间层：

在t_n时间步的求解中，先以上一时间步t_n-1的壁面温度计算结构作为本时间步的壁面热边界条件，进行边界层的气动加热计算，得到t_n时间步的热流密度接着，以该热流密度作为热边界条件，进行结构传热的计算，得到t_n时间步的壁面温度然后又作为下一时间步t_n+1流场计算的物面热边界条件。

进一步地，步骤2具体包括：

碳材料在较低温度下首先发生氧化，随着温度升高，氧化急剧增加，进入氧化扩散控制区，氧化率受边界层内氧气的输运快慢程度的控制，在更高温度下，碳-氮反应和碳的升华反应出现，引入质量烧蚀率和B无因此质量烧蚀率系数，下列计算公式中，C代表组元浓度，其下标表示组元的成分，同时下标中e与w分别表示边界层外缘条件和壁面条件，MO表示氧原子的分子量，Tw为表面温度，对于层流，对于湍流，q_or为无烧蚀情况下的热流，h_r为总焓，氧化受速率控制的过程(T_w＜1700K)，忽略该温度范围内的碳氮反应以及碳的升华反应，按照下式计算：

其中，为氧化速率系数，

为当量扩散系数，α_c(0)＝ψq_or/h_r，

氧化受氧气运输快慢控制(1700K≤T_w＜3300K)，忽略该温度范围被的碳氮反应以及碳的升华反应，按照下式计算：

在升华过程中(3300K≤T_w)，可按照下式计算：

根据能量守恒原则可以得到：

式中ΔH_CO表示一氧化碳的生成热，h表示空气焓，下标r表示总焓，

在升华烧蚀阶段，空气温度较高，碳的氧化反应、氮化反应以及升华反应都已经出现，同样根据能量守恒原则得到：

其中，

ΔH表示反应热，表示平均升华热，下标V表示碳蒸汽。

进一步地，步骤3具体包括：在步骤2中引入q_ab表示考虑烧蚀效应的热流密度，在获得了q_ab之后，将其作为第二类边界条件带入结构传热计算当中

对于热薄壁，可以改写为：

对于热厚壁，可以改写为：

表层：

式中，对于氧化烧蚀：

对于升华烧蚀：

进一步地，步骤4具体包括：在不同时刻，沿整个弹道飞行仅需要四次无粘外流解的计算状态，即4个插值文件，假设在该时刻各个状态对应的马赫数和攻角分别为：(Ma₁,α₁),(Ma₁,α₂),(Ma₂,α₁),(Ma₂,α₂)，对应的飞行器表面某点处某参数的值分别为A11,A12,A21,A22,则在所需要求解的马赫数攻角状态(Ma_x,α_x)下在参考高度上的该点处的该参数可以用以下一组公式进行计算：

本发明具有如下有益效果：

1.发明了碳基防热材料烧蚀模型及数值计算方法

基于温度分区的方法，建立了快速高效的碳基材料烧蚀工程计算方法。针对不同的材料，只需要给出相关物理参数，采用本发明的计算方法，就可对高速流动下其烧蚀特性进行研究。

2.发明了弹道状态动态插值技术

在不同时刻，沿整个弹道飞行仅需要四次无粘外流解的计算状态，即4个插值文件，采用该弹道状态动态插值技术，可以求解所需要求解的任意马赫数攻角状态下相同位置所对应的参数值。该技术增加了高超声速飞行器在高速飞行弹道过程中，热防护系统表面发生烧蚀情况的模拟能力。

3.发明了复杂高超声速飞行器烧蚀效应快速计算技术

发明了一种无粘欧拉解与工程方法相结合的数值计算方法用于快速计算与分析，复杂高超声速飞行器复杂高速流动(弹道状态)状态下，烧蚀效应对气动加热与防热结构材料传热特性的影响，给出了复杂外形高超声速飞行器热防护系统表面烧蚀的热流密度分布，烧蚀率以及热防护结构温度分布的时变特性。所发明的计算方法具有快速高效及工程设计可接受的计算精度，鲁棒性好等特点，填补了国内在高超声速飞行器发生烧蚀情况下，高效、快速和准确的给出数值模拟结果方面的空白，可为高超声速飞行器的热防护系统设计及热环境特性分析等提供技术参考与支持。

附图说明：

图1(a)为第0.05s时石墨材料表面热流密度分布。

图1(b)为第0.05s时509材料表面热流密度分布。

图2(a)为第100s时石墨材料表面热流密度分布。

图2(b)为第100s时509材料表面热流密度分布。

图3(a)为第200s时石墨材料表面热流密度分布。

图3(b)为第200s时509材料表面热流密度分布。

图4(a)为第300s时石墨材料表面热流密度分布。

图4(b)为第300s时509材料表面热流密度分布。

图5(a)为第400s时石墨材料表面热流密度分布。

图5(b)为第400s时509材料表面热流密度分布。

图6(a)为第500s时石墨材料表面热流密度分布。

图6(b)为第500s时509材料表面热流密度分布。

图7(a)为第0.05s时石墨材料表面温度分布。

图7(b)为第0.05s时509材料表面温度分布。

图8(a)为第100s时石墨材料表面温度分布。

图8(b)为第100s时509材料表面温度分布。

图9(a)为第200s时石墨材料表面温度分布。

图9(b)为第200s时509材料表面温度分布。

图10(a)为第300s时石墨材料表面温度分布。

图10(b)为第300s时509材料表面温度分布。

图11(a)为第400s时石墨材料表面温度分布。

图11(b)为第400s时509材料表面温度分布。

图12(a)为第500s时石墨材料表面温度分布。

图12(b)为第500s时509材料表面温度分布。

图13为普朗特边界层理论示意图。

图14为第500s时对称面内外层温度分布。

图15为RAMC-II对称面热流密度分布。

图16为RAMC-II对称面温度分布。

图17为RAMC-II驻点内外层温度随时间的变化。

图18为烧蚀气动加热和结构传热耦合计算步骤。

具体实施方式：

本发明复杂高超声速飞行器烧蚀效应快速计算技术，其具体包括如下步骤：

步骤1、高超声速飞行器气动加热问题。本发明采用工程方法与数值模拟计算相结合的算法：无粘外流解、高速边界层传热以及结构传热三者耦合的方法。其大致分为两个部分：其一是分析飞行器外部高超声速流场的传热，得到飞行器表面热流密度的分布；其二是表面热流进入飞行器被捕之后，分析在热防护结构内的传热问题。

步骤2、材料的烧蚀模型。本发明针对热防护材料的烧蚀问题，采用了基于碳基材料的烧蚀模型和基于温度分区的烧蚀效应计算方法，该方法能正确地表征在烧蚀过程的传质和传热机理，给出烧蚀效应较好的模拟结果。

步骤3、烧蚀效应快速计算方法。将烧蚀效应的计算方法耦合到气动加热计算方法之中，提出一种能够考虑碳基材料烧蚀效应的气动加热计算技术，从而可以得到高超声速飞行器热防护系统表面烧蚀的热流密度分布，烧蚀率以及防热结构的温度分布。

步骤4、利用动态插值技术，对弹道状态进行准动态模拟。

进一步地，步骤1具体包括：

步骤1-1、无粘外流。由于高超声速边界层厚度很薄，在求解边界层外的无粘流动时忽略边界层的厚度，采用Euler方程作为控制方程求解无粘流动，取物面参数近似作为边界层的外缘参数。其中，所需要的边界层外缘参数有：网格信息(物面网格点坐标和网格单元信息)、速度分量(物面网格点上三个方向的速度分量)、马赫数、静压、密度、静温以及总能。

步骤1-2、高速边界层传热。计算表面热流密度需分为两个部分：驻点区和非驻点区。

在驻点区，在任何情况下，驻点的流动都存在精确相似解，本发明采用Fay-Riddell方法，驻点热流密度计算公式表示如下：

式中，Pr＝0.71，Le＝1.0。

在非驻点区，热力学参数以及输运参数沿物面方向不变化的假设不成立，从而得不到精确的相似解。但是在高冷壁条件下，仍可假设，沿法向各热力学参数的变化远大于沿物面方程的变化，作局部相似。本发明主要是应用平板传热模型计算传热系数和热流密度，对于较复杂外形的计算进行适当的简化处理。

步骤1-3、气动加热流固耦合计算。飞行器表面以及热防护层的温度处在不断变化之中，而热流密度的计算需要物面温度等参数作为条件，两个重要物理量的求解各自需要彼此作为求解条件，这样就使得表面温度与热流密度直接存在耦合关系。

对于结构传热，根据热防护结构的毕奥数Bi大小，防热层可以分为热薄壁和热厚壁两种类型。

对于热薄壁，采用一维热传导方程，计算公式如下：

初始条件：T_w|_t＝0＝T₀

式中，ρ_δ为材料密度；c_δ为材料比热容；δ为材料厚度；α为传热系数；ε为材料辐射系数。对上式构建差分方程进行求解，并使用代替如下：

对于热厚壁，在计算中，由内向外热厚壁分为j层，简化后的热传导方程为：

(1)表层：

(2)最内层：

(3)中间层：

初始条件：T_j|_t＝0＝T₁|_t＝0＝T_n|_t＝0＝T₀

其中下标m为材料类型，n为第n层，λ_m为热传导系数。

同样构建差分方程，并做近似处理(为的近似代替，为的近似代替)，如下：

(1)表层：

(2)最内层：

(3)中间层：

而气动加热和结构传热要进行耦合计算——气动加热与结构传热的交叉迭代，即在t_n时间步的求解中，先以上一时间步t_n-1的壁面温度计算结构作为本时间步的壁面热边界条件，进行边界层的气动加热计算，得到t_n时间步的热流密度接着，以该热流密度作为热边界条件，对结构传热进行计算，得到t_n时间步的壁面温度然后又作为下一时间步t_n+1流场计算的物面热边界条件。如此循环往复，实现结构传热和气动热环境计算的交替进行和相互影响，达到综合考虑气动加热与结构传热的目的。

进一步地，步骤2具体包括：

碳材料在较低温度下首先发生的化学反应是氧化，起初氧化是速率控制的过程，表面反应动力学条件决定其氧化率。随着温度逐渐升高，氧化反应急剧增加，氧气逐渐供应不足，进入氧化扩散控制区，这时，氧化率受边界层内氧气的输运快慢程度的控制。随着温度进一步升高，碳-氮反应和碳的升华反应逐渐出现。引入质量烧蚀率和B无因此质量烧蚀率系数。下列计算公式中，C代表组元浓度，其下标表示组元的成分，同时下标中e与w分别表示边界层外缘条件和壁面条件，MO表示氧原子的分子量，Tw为表面温度。对于层流，对于湍流，q_or为无烧蚀情况下的热流，h_r为总焓。

氧化受速率控制的过程(T_w＜1700K)，忽略该温度范围内的碳氮反应以及碳的升华反应。可按照下式计算：

其中，为氧化速率系数。

为当量扩散系数。α_c(0)＝ψq_or/h_r。

氧化受氧气运输快慢控制(1700K≤T_w＜3300K)，忽略该温度范围被的碳氮反应以及碳的升华反应。可按照下式计算：

在升华过程中(3300K≤T_w)，可按照下式计算：

作为简化假设，认为在碳基材料的氧化烧蚀阶段，氧化产物只有，而此时空气温度较低，没有出现氧气和氮气等的解离反应。因此，根据能量守恒原则可以得到：

式中ΔH_CO表示一氧化碳的生成热，h表示空气焓，下标r表示总焓。

在升华烧蚀阶段，空气温度较高，碳的氧化反应、氮化反应以及升华反应都已经出现，同样根据能量守恒原则得到：

其中，

ΔH表示反应热，表示平均升华热，下标V表示碳蒸汽。

进一步地，步骤3具体包括：

在步骤2中提出了基于温度分区基于温度分区的烧蚀效应计算方法，给出了基于能量守恒原则表面热流密度计算的表达式。引入q_ab表示考虑烧蚀效应的热流密度，在获得了q_ab之后，可以将其作为第二类边界条件带入结构传热计算当中。

对于热薄壁，可以改写为：

对于热厚壁，可以改写为：

表层：

式中，对于氧化烧蚀：

对于升华烧蚀：

在给出初始时刻温度和材料烧蚀物理属性后，按照时间步长推进，就可以得到考虑烧蚀效应的气动加热下的结构温度分布。

进一步地，步骤4具体包括：

弹道状态是指飞行器在沿着弹道飞行的过程中，飞行器随着时间的变化改变其高度、迎角、马赫数三个参数(可将侧滑角忽略)。对整个弹道状态，选取不同的时间，通过对烧蚀效应的计算，得到各物理量在不同材料的飞行器表面的变化。下面，简要描述对无粘外流解弹道状态插值技术：

在不同时刻，均对高度效应作简化处理，沿整个弹道飞行仅需要四次无粘外流解的计算状态(根据弹道的时间、马赫数和迎角参数等的变化)，即4个插值文件。假设在该时刻各个状态对应的马赫数和攻角分别为：(Ma₁,α₁),(Ma₁,α₂),(Ma₂,α₁),(Ma₂,α₂)，对应的飞行器表面某点处某参数的值分别为A11,A12,A21,A22,则在所需要求解的马赫数攻角状态(Ma_x,α_x)下在参考高度上的该点处的该参数可以用以下一组公式进行计算：

下面结合算例对本发明进行说明。

1.技术参数

采用的算例的工程背景是巡航状态下的高超声速钝头体气动加热问题。计算模型是RAMC-II，是一个总长1.295m，头部半径为0.1525m，半锥角9°的顿锥。设定巡航状态：马赫数15、巡航高度70km、巡航时间500；热防护材料：热薄壁(石墨、厚度2mm)+隔热层(SiO₂、厚度20mm)，表面发射率0.8。该算例所得带烧蚀效应的气动热结果，将对比不同烧蚀材料对于烧蚀效应计算的影响。热防护设置如下表所示：

2.数值模拟结果

整个算例的数值模拟在普通个人PC上运行，算例模拟计算时间为CPU耗时29分钟，国内暂无出现同等高效的高超声速飞行器烧蚀效应计算技术，属国内首次，下面为计算结果分析：

从图1(a)～图6(b)中可以发现，不论石墨材料还是509材料，由于碳基材料以及隔热层的导热性较差，表面热流密度在很短的时间内就达到了平衡；对比相同时间的表面热流密度，石墨材料的表面热流密度要高于509材料的。从图7(a)～图12(b)可以发现，在烧蚀过程中，石墨材料的烧蚀过程是一个放热过程，固态燃烧产生一氧化碳，释放热量，导致表面热流密度高于没有烧蚀的表面热流密度。相反，509材料的烧蚀过程是一个吸热过程，材料热解吸收的热量不仅超过了燃烧产生的热量，还包括了没有烧蚀的气动加热量，导致了表面热流密度的减小。

从图14中可以看出，在对称面上，表面温度明显高于内层温度，头部区域的最高温差达到1700K，在头部与身部的过渡区，温度梯度最大。

从图15和图16中可以看出，热流密度分布的不同导致了表面温度分布的不同。在初始时刻，两种烧蚀材料的表面温度都很低，随着时间的增加，加热效果明显，高温区逐步向下游地区延伸；在任意相同的时刻、相同位置，由于石墨烧蚀是放热过程，而509烧蚀是吸热过程，石墨材料的表面温度要高于509材料。

从图17可以看出，驻点表层温度在巡航的初始阶段就达到了平衡，石墨材料表面温度升高的趋势较509材料更剧烈一些，石墨材料的平衡温度在2500K左右，而509材料的平衡温度较低，在2000K左右。两种材料的内层温度升高趋势较表层温度要缓和得多，且都还没有到达平衡值。因为509材料的表层温度要低于石墨材料，所以509材料的内层温度也低于石墨材料。

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

Claims (5)

1.一种复杂高超声速飞行器烧蚀效应快速计算技术，其特征在于：包括如下步骤

步骤1、高超声速飞行器气动加热；

步骤2、材料的烧蚀模型；

步骤3、烧蚀效应快速计算；

步骤4、利用动态插值技术，对弹道状态进行准动态模拟。

2.如权利要求1所述的复杂高超声速飞行器烧蚀效应快速计算技术，其特征在于：步骤1具体包括

步骤1-1、无粘外流，采用欧拉方程作为控制方程求解无粘流动，取物面参数近似作为边界层的外缘参数，其中，所需要的边界层外缘参数有：网格信息、速度分量、马赫数、静压、密度、静温以及总能；

步骤1-2、高速边界层传热，计算表面热流密度需分为两个部分：驻点区和非驻点区，在驻点区，在任何情况下，驻点的流动都存在精确相似解，采用Fay-Riddell方法，驻点热流密度计算公式表示如下：

$q_{ws}=0.763{\mathrm{Pr}}^{-0.6}{\left(\frac{{\mathrm{\ρ}}_w{\mathrm{\μ}}_w}{{\mathrm{\ρ}}_s{\mathrm{\μ}}_s}\right)}^{0.1}\sqrt{{\mathrm{\ρ}}_s{\mathrm{\μ}}_s{\left(\frac{{\mathrm{du}}_e}{dx}\right)}_s}\·\[1+({\mathrm{Le}}^{0.52}-1)\frac{h_D}{h_s}\](h_s-h_w)$

式中，Pr＝0.71，Le＝1.0；式中ρ_w为物面密度，μ_w为物面粘性系数,h_w为物面焓值，ρ_s为驻点处密度，μ_s为驻点处粘性系数，h_D为平均空气电离焓；

步骤1-3、气动加热流固耦合计算，对于结构传热，根据热防护结构的毕奥数Bi大小，防热层可以分为热薄壁和热厚壁两种类型；

对于热薄壁，采用一维热传导方程，计算公式如下：

${\mathrm{\ρ}}_{\mathrm{\δ}}{c}_{\mathrm{\δ}}\mathrm{\δ}\frac{dT}{dt}=\mathrm{\α}(T_{\mathrm{\α}w}-T_w)-{\mathrm{\ϵ\σT}}_w^4$

初始条件：T_w|_t＝0＝T₀

式中，ρ_δ为材料密度；c_δ为材料比热容；δ为材料厚度；α为传热系数；ε为材料辐射系数；

对上式构建差分方程进行求解，并使用代替如下：

$T_w=\frac{\[{\mathrm{\αT}}_{\mathrm{\α}w}-\mathrm{\ϵ}\mathrm{\σ}{\left(T_w^{k-1}\right)}^4\]\frac{\mathrm{\Δ}t}{{\mathrm{\ρ}}_{\mathrm{\δ}}{c}_{\mathrm{\δ}}\mathrm{\δ}}+T_w^{k-1}}{1+\mathrm{\α}\frac{\mathrm{\Δ}t}{{\mathrm{\ρ}}_{\mathrm{\δ}}{c}_{\mathrm{\δ}}\mathrm{\δ}}}$

对于热厚壁，在计算中，由内向外热厚壁分为j层，简化后的热传导方程为：

(1)表层：

(2)最内层：

(3)中间层：

初始条件：T_j|_t＝0＝T₁|_t＝0＝T_n|_t＝0＝T₀

其中下标m为材料类型，n为第n层，λ_m为热传导系数；

同样构建差分方程，并做近似处理为的近似代替，为的近似代替，如下：

(1)表层：

(2)最内层：

(3)中间层：

在t_n时间步的求解中，先以上一时间步t_n-1的壁面温度计算结构作为本时间步的壁面热边界条件，进行边界层的气动加热计算，得到t_n时间步的热流密度接着，以该热流密度作为热边界条件，进行结构传热的计算，得到t_n时间步的壁面温度然后又作为下一时间步t_n+1流场计算的物面热边界条件。

3.如权利要求2所述的复杂高超声速飞行器烧蚀效应快速计算技术，其特征在于：步骤2具体包括：

碳材料在较低温度下首先发生氧化，随着温度升高，氧化急剧增加，进入氧化扩散控制区，氧化率受边界层内氧气的输运快慢程度的控制，在更高温度下，碳-氮反应和碳的升华反应出现，引入质量烧蚀率和B无因此质量烧蚀率系数，下列计算公式中，C代表组元浓度，其下标表示组元的成分，同时下标中e与w分别表示边界层外缘条件和壁面条件，MO表示氧原子的分子量，Tw为表面温度，对于层流，对于湍流，q_or为无烧蚀情况下的热流，h_r为总焓，氧化受速率控制的过程(T_w＜1700K)，忽略该温度范围内的碳氮反应以及碳的升华反应，按照下式计算：

$\left\{\begin{array}{c}{\stackrel{\·}{m}}_C=\frac{{\tilde{C}}_{O_{2}e}}{1+B}/\left(\frac{1}{1/{\stackrel{\‾}{\mathrm{\α}}}_{ch}+1/{\mathrm{\α}}_{d1}}\right)\\ B=\frac{{\stackrel{\·}{m}}_C}{{\mathrm{\ψq}}_{or}/h_r}\\ {\stackrel{\‾}{\mathrm{\α}}}_{ch}=\frac{{\mathrm{\α}}_{ch}(1+B)}{\left(C_{\stackrel{\‾}{O}e}-\frac{M_O}{M_C}B\right)}\\ {\mathrm{\α}}_{d1}=\frac{M_C}{M_O}(1+B){\mathrm{\α}}_c\left(0\right)\end{array}\right.$

其中，为氧化速率系数，

为当量扩散系数，α_c(0)＝ψq_or/h_r，

氧化受氧气运输快慢控制(1700K≤T_w＜3300K)，忽略该温度范围被的碳氮反应以及碳的升华反应，按照下式计算：

$\left\{\begin{array}{c}B=\frac{M_C}{M_O}{\tilde{C}}_{O_{2}e}=0.1725\\ {\stackrel{\·}{m}}_C=B\·\frac{{\mathrm{\ψq}}_{or}}{h_r}\end{array}\right.$

在升华过程中(3300K≤T_w)，可按照下式计算：

$B=\frac{{\[\underset{i=1}{\overset{3}{\Σ}}{C}_i+\frac{M_C}{M_{CO}}{C}_{CO}+\frac{M_C}{M_{CN}}{C}_{CN}+\frac{M_C}{M_{C_{2}N}}{C}_{C_{2}N}\]}_w}{1-{\[\underset{i=1}{\overset{3}{\Σ}}{C}_i+\frac{M_C}{M_{CO}}{C}_{CO}+\frac{M_C}{M_{CN}}{C}_{CN}+\frac{M_C}{M_{C_{2}N}}{C}_{C_{2}N}\]}_w}$

${\stackrel{\·}{m}}_C={\mathrm{B\ψq}}_{or}/h_r$

根据能量守恒原则可以得到：

$q_b={\mathrm{\ψq}}_{or}1-h_w/h_r-{\mathrm{\ϵ\σT}}_w^4+{\stackrel{\·}{m}}_C{\mathrm{\ΔH}}_{CO}$

式中ΔH_CO表示一氧化碳的生成热，h表示空气焓，下标r表示总焓，

在升华烧蚀阶段，空气温度较高，碳的氧化反应、氮化反应以及升华反应都已经出现，同样根据能量守恒原则得到：

$\begin{array}{c}{q}_b=\frac{{\mathrm{\ψq}}_{or}}{h_r}\[h_r-h_w+(1+B)\frac{M_O}{M_{CO}}{C}_{{\mathrm{CO}}_w}({\mathrm{\ΔH}}_{CO}\\ -\frac{M_C}{M_O}\mathrm{\Δ}{\stackrel{\‾}{H}}_V-{\mathrm{\ΔH}}_{O_2})+(1+B)\frac{M_N}{M_{CN}}{C}_{{\mathrm{CN}}_w}({\mathrm{\ΔH}}_{CN}\\ -\frac{M_C}{M_N}\mathrm{\Δ}{\stackrel{\‾}{H}}_V-{\mathrm{\ΔH}}_{N_2})-B\mathrm{\Δ}{\stackrel{\‾}{H}}_V+(1+B)\frac{M_N}{M_{C_{2}N}}\\ C_{C_2{N}_w}({\mathrm{\ΔH}}_{C_{2}N}-\frac{M_C}{M_N}\mathrm{\Δ}{\stackrel{\‾}{H}}_V-{\mathrm{\ΔH}}_{N_2})\]-{\mathrm{\ϵ\σT}}_w^4\end{array}$

其中，

$\mathrm{\Δ}{\stackrel{\‾}{H}}_V=\underset{i=1}{\overset{3}{\Σ}}{\mathrm{\λ}}_i{\mathrm{\ΔH}}_{V_i}$

${\mathrm{\λ}}_i=\frac{C_{C_i}}{\underset{i=1}{\overset{3}{\Σ}}{C}_{C_i}}$

ΔH表示反应热，表示平均升华热，下标V表示碳蒸汽。

4.如权利要求3所述的复杂高超声速飞行器烧蚀效应快速计算技术，其特征在于：步骤3具体包括：在步骤2中引入q_ab表示考虑烧蚀效应的热流密度，在获得了q_ab之后，将其作为第二类边界条件带入结构传热计算当中

对于热薄壁，可以改写为：

$T_w=\frac{\[q_{ab}-\mathrm{\ϵ}\mathrm{\σ}{\left(T_w^{k-1}\right)}^4\]\mathrm{\Δ}t}{{\mathrm{\ρ}}_{\mathrm{\δ}}{c}_{\mathrm{\δ}}\mathrm{\δ}}+T_w^{k-1}$

对于热厚壁，可以改写为：

表层：

式中，对于氧化烧蚀：

$q_{ab}={\mathrm{\ψq}}_{or}1-h_w/h_r+{\stackrel{\·}{m}}_C{\mathrm{\ΔH}}_{CO}$

对于升华烧蚀：

$\begin{array}{c}{q}_{ab}=\frac{{\mathrm{\ψq}}_{or}}{h_r}\[h_r-h_w+(1+B)\frac{M_O}{M_{CO}}{C}_{{\mathrm{CO}}_w}({\mathrm{\ΔH}}_{CO}-\frac{M_C}{M_O}\mathrm{\Δ}{\stackrel{\‾}{H}}_V-{\mathrm{\ΔH}}_{O_2})\\ +(1+B)\frac{M_N}{M_{CN}}{C}_{{\mathrm{CN}}_w}({\mathrm{\ΔH}}_{CN}-\frac{M_C}{M_N}\mathrm{\Δ}{\stackrel{\‾}{H}}_V-{\mathrm{\ΔH}}_{N_2})-B\mathrm{\Δ}{\stackrel{\‾}{H}}_V+\\ (1+B)\frac{M_N}{M_{C_{2}N}}{C}_{C_2{N}_w}({\mathrm{\ΔH}}_{C_{2}N}-\frac{M_C}{M_N}\mathrm{\Δ}{\stackrel{\‾}{H}}_V-{\mathrm{\ΔH}}_{N_2})\]\end{array}.$

5.如权利要求4所述的复杂高超声速飞行器烧蚀效应快速计算技术，其特征在于：在不同时刻，沿整个弹道飞行仅需要四次无粘外流解的计算状态，即4个插值文件，假设在该时刻各个状态对应的马赫数和攻角分别为：(Ma₁,α₁),(Ma₁,α₂),(Ma₂,α₁),(Ma₂,α₂)，对应的飞行器表面某点处某参数的值分别为A11,A12,A21,A22,则在所需要求解的马赫数攻角状态(Ma_x,α_x)下在参考高度上的该点处的该参数可以用以下一组公式进行计算：

$A_{x1}=A_{11}*\frac{M_2-M_x}{{\mathrm{Ma}}_2-{\mathrm{Ma}}_1}+A_{21}*\frac{M_x-M_1}{{\mathrm{Ma}}_2-{\mathrm{Ma}}_1}$

$A_{x2}=A_{12}*\frac{M_2-M_x}{{\mathrm{Ma}}_2-{\mathrm{Ma}}_1}+A_{22}*\frac{M_x-M_1}{{\mathrm{Ma}}_2-{\mathrm{Ma}}_1}$

$A_{xx}=A_{x1}*\frac{{\mathrm{\α}}_2-{\mathrm{\α}}_x}{{\mathrm{\α}}_2-{\mathrm{\α}}_1}+A_{x2}*\frac{{\mathrm{\α}}_x-{\mathrm{\α}}_1}{{\mathrm{\α}}_2-{\mathrm{\α}}_1}.$