一种应用于烧蚀型热防护系统的设计方法

技术领域

本发明属于高超声速热防护设计相关领域，具体涉及一种应用于烧蚀型热防护系统的设计方法。

背景技术

高速飞行器返回再入及其在大气层内飞行的过程中面临多样化的气动加热环境，热防护系统作为高速飞行器最为重要的子系统之一，在保障人员装备安全和飞行器结构完整性方面具有不可替代的作用。为应对不同的热环境特征，目前发展了烧蚀型和非烧蚀型两大类热防护系统，其中烧蚀型热防护系统所用材料是一种低密度，在高温下会发生热解反应吸收热量并产生热解气体的热解炭型材料，其防热机制主要包括热辐射、热容吸热、热解反应吸热以及热阻塞效应，可通过材料热解反应吸热、热解气体释放以及表面烧蚀后退达到减少热量传导的作用，从而保护飞行器内部结构。烧蚀型热防护系统结构与成型工艺较为简单，具有高效、轻质、工艺简单、成本低等优点，从近程弹道导弹到远程弹道导弹、飞船、返回式卫星以及低升力再入体均采用烧蚀防热，从长时间到几十秒时间、高气动加热率都能适用。经过几十年的发展，烧蚀材料目前已成为应用最成熟、最广泛的热防护材料体系之一。

由于弹道偏差、来流条件的不确定性、材料属性(特别是烧蚀型热防护材料的热物理和热化学材料属性)的分散性、飞行器外形参数偏差等各类因素，热防护系统的设计和使用中通常存在不同程度的不确定性，给系统的可靠性设计与评估带来极大挑战。由于在分析方法、试验测量和设计手段等方面存在不足，传统设计中通常采用安全系数或设计裕度来考虑不确定性的影响，但这种方法不能准确考虑不确定性分布对热防护系统的影响，通常会导致偏保守的设计；在某些特殊情况下，又不能充分考虑敏感参数不确定性的影响，使得系统可靠性降低甚至不能满足工程需求。

发明内容

本发明意在针对现有设计方法的不足，提供一种基于概率技术的精细化、符合实际情况的飞行器烧蚀型热防护系统设计方法。

本发明中的应用于烧蚀型热防护系统的设计方法，包括如下步骤：

S1、根据所给飞行器的外形尺寸和飞行历程设计参数以及热防护材料选型，初步确定烧蚀型热防护系统的各设计参数，其中包括烧蚀型热防护系统的热响应相关参数；

所述热响应相关参数包括几何尺寸参数、材料属性参数、飞行弹道参数和来流条件参数。

S2、根据所述热响应相关参数确定其中具有不确定性的参数，并对这些参数进行不确定性量化表征；

S3、对N个具有不确定性的参数，按照每个参数的不确定性量化表征所对应的概率分布规律进行M次随机抽样，得到N组，每组有M个值的样本序列；

S4、从这N组样本序列中各随机选取一个样本值，得到一个随机输入样本，再从各样本序列剩下的值中再各随机选取一个样本值得到新的一个随机输入样本，以此类推可得到M个随机输入样本且任意不确定性参数的值在M个随机输入样本中仅出现一次；

S5、建立烧蚀型热防护系统气动热以及热响应的参数化分析模型；

对所述分析模型进行网格化处理，根据设计尺寸分别做出用于气动热计算的标称面网格以及用于热响应计算的标称体网格；

根据网格变形算法从标称面网格/体网格得到每组随机输入样本对应的面网格/体网格；

S6、依据从步骤S4中所得到的随机输入样本中的来流条件参数和其他具有确定性的来流条件参数，采用有限元分析工具对所述热防护系统气动热分析模型进行空间离散和分析计算，并代入当前时刻下的来流条件参数得到M个的线性微分方程组，求解得到M个当前时刻下的冷壁热流

对所给飞行器整个飞行历程下的t个时刻求解得到M个飞行历程下的冷壁热流q

S7、依据从步骤S4中所得到的随机输入样本中的的材料属性参数和几何尺寸参数，采用有限元分析工具对热防护系统的热响应分析模型进行空间离散和分析计算，其中的热载荷依据步骤S6计算同一样本下的冷壁热流转换为的热壁热流，共得到M个包含时间t为独立变量的非线性微分方程组；

S8、求解步骤S7中得到的M个非线性微分方程组，得到M个温度场向量T(x,y,z,t)，即整个热防护系统空间内所有有限节点在时间历程内每个计算时刻的温度结果；

S9、从这M个温度场向量T(x,y,z,t)中提取出热防护系统的背温区域在整个时间历程内的最高温度作为这M个不确定性输入样本所对应的输出温度特征量；

S10、由步骤S9所得到的M个不确定性输入参数组合样本及其对应的输出温度特征量，构建和训练用于拟合随机输入样本变量与输出温度特征量间函数关系的代理模型，并对代理模型的精度进行评价，若不满足精度要求，则返回并重复步骤S3-S10，直至所得模型的精度达到精度要求；

S11、再次对各个具有不确定性的参数进行n次抽样，并利用步骤S10中得到代理模型得到n个目标随机输出温度特征量，进而计算该设计参数情况下的热防护系统的可靠度；

S12、当热防护系统的可靠度大于或等于所要求的数值时则结束；

否则调整热防护系统中的设计参数，重复步骤S2～S11，直到热防护系统的可靠度大于或等与所要求的数值。

进一步的，所述材料属性参数中具有不确定性的参数包括原始材料的密度，原始材料与炭化材料的比热容、热导率，孔隙率，气体渗透率以及材料热解热；

所述几何尺寸参数的不确定性量化表征根据光学扫描仪对热防护涂层的测量获得；

所述来流条件参数中具有不确定性的参数包括飞行器的整个飞行历程中的高度、速度、攻角。

优选的，所述代理模型为Kriging模型。

进一步的，步骤S3和步骤S11中所述抽样的方法为基于拉丁超立方的蒙特卡洛抽样方法，具体包括，对某个具有不确定性的参数，在其取值范围内分为M个概率相同的区间，从这M个概率相同的区间内随机选取一个值，获得的M个值组成该参数的一组样本序列，共得到N组不确定性参数的样本序列。

进一步的，步骤S8中，采用牛顿迭代法与隐式时间积分方法求解非线性微分方程组。

进一步的，步骤S7中，对所述热响应模型的分析计算中边界条件处的气动热流加载方式为；

式中C

进一步的，步骤S7中，所述热壁热流密度计算公式为；

式中q

进一步的，步骤S11中热防护系统的可靠度的计算公式为：

式中T

进一步的，步骤S11中，还包括对具有不确定性的参数进行灵敏度分析的步骤。

进一步的，所述灵敏度分析采用Sobol指标计算，包括，基于抽样生成n×2N的样本矩阵，前N列为矩阵A，后N列为矩阵B；

构建矩阵

Var(Y)＝Var(f(A)+f(B))

/>

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

该设计方法

(1)、本发明提出一种可针对烧蚀型热防护系统的概率设计和可靠性评估方法，可有效减少烧蚀型热防护系统的设计余量，从而降低热防护系统的重量。

(2)、本发明考虑了多种不确定性参数，包括几何偏差、弹道偏差、来流偏差以及热防护材料属性的不确定性，综合全面的设计分析该烧蚀型热防护系统。

(3)、本发明使用工程算法计算瞬态气动热载荷，可快速近似求解气动热，从而快速获得热防护系统的热流载荷。

因此，本发明对参数的不确定性进行了准确量化和表征，基于概率设计方法合理确定烧蚀材料的厚度，在保证可靠性的前提下尽最大可能减轻热防护系统重量，提高了热防护系统效率，满足了飞行器轻量化和高效防热的要求。

附图说明

图1为本发明烧蚀型热防护系统概率设计方法分析流程图；

图2为本发明实施例中的气动热计算有限元网格示意图；

图3为本发明实施例中的热响应二维轴对称有限元网格示意图；

图4为本发明实施例中的材料相关属性曲线；

图5为本发明实施例中的材料相关属性曲线；

图6为本发明实施例中的材料相关属性曲线；

图7为本发明实施例中的材料相关属性曲线；

图8为本发明实施例中的材料相关属性曲线；

图9为本发明实施例中代理模型的计算结果与确定性计算结果的对比图；

图10为本发明实施例中的概率特性分布与累计概率分布曲线图。

具体实施方式

下面结合附图和具体实施方式对本发明进行详细说明。

高超声速飞行器的热防护系统处于最外侧，烧蚀型热防护系统材料主要分为硅基、碳基以及热解碳化型。在近年来，航天返回舱以及弹头再入大气层过程中广泛的使用硅橡胶这类热解炭化材料，其主要热耗散机制为热辐射、热容吸热、热解反应吸热以及热阻塞效应，对其内部结构进行防隔热保护。

如图1所示，本实施例中的应用于烧蚀型热防护系统的设计方法，包括如下步骤：

S1、根据飞行器的外形尺寸和飞行历程设计参数，以及热防护材料选型，初步确定烧蚀型热防护系统的各设计参数，其中包括烧蚀型热防护系统的热响应相关参数；包括几何尺寸产出、飞行弹道参数、来流条件参数、材料属性参数等。

S2、根据上述的热防护系统相关参数，确定其中具有不确定性的参数并对相应参数进行不确定性量化表征，其中材料属性参数中具有不确定性的参数主要为原始材料的密度，原始材料与炭化材料的比热容(焓值)、热导率，孔隙率，气体渗透率，材料热解热等，可采用试验或相关文献对其进行不确定性量化表征；

几何尺寸参数的不确定性表征根据光学扫描仪对热防护涂层的测量获得，并进行不确定性量化表征。来流条件参数为飞行器的整个飞行历程中的高度、速度、攻角等，本实施例中结合相关文献对这些不参数进行不确定性量化表征。

S3、按照上述所表征的每个具有不确定性的参数的概率分布规律进行随机抽样，采用拉丁超立方抽样，对N个参数，将每个参数的取值范围分为M个概率相同的区间，从这M个概率相同的区间内随机选取一个值，获得的M个值组成对应不确定性参数的一组样本序列，共得到N组不确定性参数的样本序列。

S4、从这N组样本序列中各随机选取一个样本值，得到一个随机输入样本，再从各样本序列剩下的值中再各随机选取一个样本值得到新的一个随机输入样本，以此类推可得到M个随机输入样本且任意参数的值在M个随机输入样本中仅出现一次。

S5、由步骤S4所得随机输入样本以及其他具有确定性的参数建立烧蚀型热防护系统气动热以及热响应的参数化模型，首先对热防护系统的气动热以及热响应分析模型进行网格化处理，根据设计尺寸分别做出分别如图2和图3所示的用于气动热计算的面网格以及热响应计算的体网格，再根据各随机输入样本中的几何尺寸参数以及网格变形算法得到每组随机输入样本各自的有限元网格。

S6、根据飞行弹道、来流条件以及飞行器表面网格，采用有限元分析工具进行热防护系统的气动热分析。

首先，利用牛顿修正理论计算高超声速飞行器表面压力分布，采用非结构网格对飞行器模型表面进行划分，在每个三角形面元内利用修正牛顿理论，假设来流在面元法向上的动量分量完全损失掉，根据牛顿第二运动定律，求解面元内的气动压力，最后叠加获得飞行器模型表面上的气动压力分布。基于牛顿最速下降方法，计算高超声速飞行器表面流线轨迹，利用之前所计算的压力分布及参考焓法和考虑高温气体效应情况下气体参数的拟合公式，计算流线轨迹上各数据点表面参考焓条件下的气体参数，沿流线积分获得目标点处的雷诺数，再把相关参数带入计算公式，获得目标点处的热流密度。

本实施例中，其具体过程为：

S6.1用非结构网格对所研究的模型进行网格划分，输出所有节点坐标并按一定顺序输出每个三角形面元三个节点的节点编号。

S6.2给定来流条件等输入参数，使其能够满足气动力计算程序的输入需求，然后读入计算气动力子程序，采用修正牛顿理论进行计算。气动力计算程序的计算是以面元为对象的，计算输出结果为每个节点上的气动压力。

S6.3以S6.1输出的节点及面元信息和S6.2计算得到的气动压力为输入项，带入热流密度计算程序。

热流密度计算程序的计算流程如下：

S6.3.1计算飞行器表面经过目标点的流线轨迹；

S6.3.2利用正激波关系式计算驻点处激波波后的焓，根据上一步计算流线轨迹确定驻点位置及驻点处压强，确定驻点处的熵值；

S6.3.3利用等熵假设，利用已求得的飞行器表面压力分布，结合高温气体效应，计算流线轨迹上各离散点处的焓值分布；

S6.3.4根据已知条件和以求结果，计算流线轨迹上各离散点处的气流速度、壁面焓、绝热壁面焓和参考焓；

S6.3.5利用参考焓和压力分布，计算流线轨迹上各离散点处的气流密度、温度、声速及粘性系数；

S6.3.6根据之前的计算结果由下式计算目标点处的雷诺数：

其中，x

S6.3.7将雷诺数带入计算式，获得目标点出的热流密度。

需要说明的是，这样计算出的热流是冷壁热流，在后续的有限元计算中的气动热载荷部分，需要将其转换为热壁热流，同时考虑辐射的影响来进行热响应分析。

气动热载荷的求解采用每个时刻的稳态来近似计算瞬态气动热，具体分析流程如下：

①输入弹道数据；

②读取弹道当前时间步的高度、速度、攻角到气动热计算程序，计算稳态情况下的模型整体气动热，输出当前时步计算结果；

③推进到下一时间步，重复流程②；

④整个弹道时间步计算完成，对热流数据在时间维度上进行插值，得到整体冷壁热流密度矩阵q

S7、从上述步骤中所得到随机输入样本，依据其中的材料属性参数、几何尺寸参数采用有限元分析软件(SAMCEF-Thermal/Amaryllis)对热防护系统进行热响应计算分析，其中所需的热载荷依据步骤S6计算得到的同一样本下的冷壁热流转换为热壁热流来近似，共得到M个包含时间t为独立变量的非线性偏微分方程组。

以硅橡胶这类热解炭化型热防护材料的热结构瞬态热分析模型为例，主要涉及经典热传导、多孔介质流体、材料热解反应和表面烧蚀后退，相应的机理方程如下所示。

1)、经典热传导：热平衡方程，傅里叶定律

当满足热平衡方程时，流入和流出控制体积的传导热的差异必须等于源项加上由材料加热或冷却引起的容性热。瞬时温度分布T(x，t)在固体占据体积V并满足热平衡方程时由下式给出：

在这个方程中，ρh是材料的焓值也可以通过热容量来表示，Q是体积源项，q是傅里叶定律与局部温度梯度相关的热流通量矢量。

其中q和

对于热解状态下的材料参数，如热导率可视为从原始状态线性到炭化状态如下式所示：

λ＝λ

式中λ

2)、多孔介质流体：达西定律

热解气体的质量流量和压力之间的关系通过下式给出：

气体流动项对热方程的影响：

其中K

M

3)、材料热解反应：阿伦尼乌斯定理

热解炭化材料在高温下会进行热化学分解吸收热量并产生热解气体，热解动力学方程一般采用阿伦尼乌斯方程，单一阿伦尼乌斯定律的热解反应速率如下所示：

其中密度ρ当前状态下的密度，下标c和v分别代表"碳化"和"原始"状态。A是反应的频率(指前因子)，E是活化能，N是反应的阶数，R是气体常数。材料热解吸收对热平衡方程的影响如下：

式中

热解气体的释放会产生热阻塞效应，进一步减少热流的加载，减少的热流如下式所示：

式中η

4)、表面烧蚀后退

烧蚀是由于材料所承受的热和/或机械载荷而从结构表面去除的过程，因为烧蚀现象会导致材料的去除，这将需要结构网格重新划分。主要分为三种不同类型的表面烧蚀，即机械烧蚀、化学烧蚀、相变烧蚀。

①机械烧蚀

在这种情况下，由于与周围"空气"的机械摩擦，材料被移除，不会吸收带走热量。表面烧蚀率一般取决于温度、压力和剪切应力如下式所示：

其中a，b和T

②化学烧蚀

指材料表面在高温高压的条件下发生氧化、升华以及碳氮之间的反应等现象造成表面材料发生质量损失，从而带走大量的热量，同时化学烧蚀会进一步的产生阻塞效应减小热量的加载

③相变烧蚀

当表面材料达到相变温度T

热响应分析模型的气动热载荷采用的冷热壁面转换公式如下：

式中q

SAMCEF-Thermal/Amaryllis软件中边界条件处的气动热流加载方式通过下式转换得到；

式中C

S8、采用牛顿迭代法与时间隐式积分方法求解步骤S7中得到的M个非线性偏微分方程组，得到M个温度场向量T(x,y,z,t)即整个热防护系统空间内所有有限节点在时间历程内每个计算时刻的温度结果，具体求解方法如下所示：

求解的离散化方程为：

其中C(q)是离散化系统的热容矩阵，K(q)是热传导矩阵，q是未知的温度/压力矢量；g

以一阶隐式时间积分法求解传热方程，在此方法中，节点变量的变化率在一个时间步长中被认为是恒定的。

求解方案按以下方式进行；

时域被划分为n和n+1之间的连续时间间隔，其中单个时间：

t

解向量可以以相同的方式表示，得到以下表达式：

q

节点变量的变化率为：

将上述表达式代入传热方程，可以得到以下表达式：

使用牛顿迭代法求解，采用如下方式计算：

i+1

校正的计算方法如下：

残余值定义为：

迭代计算到达收敛性，使用以下表达式来计算时间n+1处的解。

S9、从这M个温度场向量T(x,y,z,t)中提取出热防护系统的背温区域在整个时间历程内的最高温度作为这M个不确定性的随机输入样本所对应的输出温度特征量。

S10、由步骤S9所得到的M个随机输入样本及其对应的输出温度特征量，构建并训练用于拟合随机输入样本与输出温度特征量间的函数关系的代理模型，本例中采用但不限于Kriging模型。

本实施例中的Kriging模型由全局模型与局部偏差迭加而成，可表示为

Y(x)＝f(x)+Z(x) (25)

式中Y(x)为未知的近似模型，f(x)已知的多项式函数，Z(x)是均值为零、方差为σ

Cov[Z(x

式中R为相关矩阵，R(x

/>

式中n

此处y是长度为n

式中

全局模型的方差估计值为

通过极大似然估计确定相关参数θ

当θ

代理模型的精度评价公式如下所示：

式中

可将M组随机输入样本分为训练样本和测试样本，如果测得代理模型拟合的精度不满足精度要求(低于设定的门限值)，则还需返回重新进行步骤S3-10，直至精度满足精度要求；

S11、利用步骤S10中构建的Kriging代理模型进行n次基于拉丁超立方抽样的蒙特卡洛模拟，得到n个目标随机输出温度特征量即最高背温，用于计算该设计参数情况下的热防护系统的可靠度；热防护系统的可靠度的计算公式为：

式中T

通过Kriging代理模型进行n次基于拉丁超立方抽样的蒙特卡洛模拟还可以计算得到热防护系统最高背温的概率分布特性，其均值和标准差可估计为：

通过代理模型进行不确定性参数灵敏度分析，采用Sobol指标计算。基于拉丁超立方的蒙特卡洛抽样方法生成n×2N的样本矩阵，前N列为矩阵A，后N列为矩阵B。进一步构建矩阵

Var(Y)＝Var(f(A)+f(B)) (37)

S12、判断热防护系统的可靠度是否大于或等于所要求的数值，当热防护系统的可靠度大于或等于所要求的数值时则结束；否则，调整热防护系统中设计参数，重复步骤(2)～(11)，直到热防护系统的可靠度大于或等与设置的最低要求。

本实施例中，以一球形构件的烧蚀型热防护系统设计为例，利用本发明所提供的方法进行概率分析和参数设计。其中热防护结构模型采用二维轴对称网格如图3所示，考虑来流条件，材料属性以及几何尺寸的不确定性参数。

其中来流条件参数中的具有不确定性的参数为来流静温、密度和速度，通过查阅相关文献资料可得来流条件参数的分散性可达到20％，于是设置来流条件参数的标准差为0.05。通过试验以及相关文献对材料属性参数进行不确定性量化表征如表1所示。球头模型的几何偏差主要考虑基体外形半径以及涂层材料的厚度如表2所示。

表1材料属性参数及其不确定性量化表征

表2几何尺寸参数及其不确定性定量表征

利用实施例中的方法进行随机抽样仿真，进行100次随机抽样，得到100组随机不确定性参数与最高背温样本数据。以50组为训练样本构建代理模型，剩余50组为测试样本，对代理模型进行精度评价R

表3球头模型中具有不确定性的参数的Sobol指标

/>