一种基于概率技术的热防护系统设计方法与流程

本发明涉及一种基于概率技术的热防护设计方法，属于高超声速热防护设计相关领域。

背景技术：

热防护系统(Thermal Protection System,TPS)是高超声速飞行器的最重要的系统之一，在飞行器气动加热的高温恶劣环境下保障系统的结构完整性和飞行器内部的温度条件。热防护系统的可靠性和结构完整性对飞行器的安全至关重要，由其故障而造成的事故在航天/飞行器任务失败中占相当比例，因此热防护技术是高超飞行器设计中最为关键的环节。

在热防护系统的设计过程中，最为关键的问题是合理确定热防护系统的材料种类、结构形式及其尺寸，在保证热防护系统防隔热效果的前提下降低系统的重量。然而热防护系统的设计输入条件的存在很多不确定性，如飞行轨道偏差、分析模型误差、数据测量误差、材料属性分散性、系统加工装配偏差、外部载荷环境等，这些不确定参数及其分布会导致热防护设计结果存在不确定性。在传统的确定性热防护系统设计方法中，通过设置安全系数来考虑设计过程中不确定性的不利影响，这种方法得到的设计结果往往过于保守，导致系统重量过大，而在某些情况下又由于未考虑系统性能对设计参数或不确定参数的灵敏度，导致系统可靠性降低。

技术实现要素：

本发明的技术解决问题是：克服现有技术的不足，提供一种基于概率技术的科学、合理、可靠的飞行器热防护设计方法。

本发明的技术解决方案是：一种基于概率技术的飞行器热防护系统设计方法，所述热防护系统由多个热防护结构件拼接而成，处于飞行器最外侧，对其内部结构进行防隔热保护，其外形尺寸与飞行器的外形尺寸一致，该方法包括如下步骤：

(1)、根据飞行器的外形尺寸和热环境条件，初步确定每个热防护结构的设计参数，所述设计参数包括材料选型、结构形式、几何尺寸，所述热环境条件是指飞行器飞行过程中的气动热流密度；

(2)、根据热防护系统中每个热防护结构件的设计参数，确定热防护系统的热随机变量及每个热随机变量的概率分布，所述热随机变量包括几何尺寸，材料的密度、比热容、热导率、不同结构件之间界面热阻、热防护系统表面辐射率；

(3)、按照每个热随机变量的概率分布规律，将每个热随机变量的取值范围分为m个等概率区间，对于每个热随机变量，在各等概率区间内随机选取一个值，组成一组样本序列，共得到N个热随机变量的样本序列；

(4)、从每个热随机变量样本序列中随机抽取一个样本值，形成所有随机变量的一个样本组合，再从每个热随机变量样本序列剩余的数据中随机抽取一个样本值，形成又一个样本组合，以此类推，进行m次随机抽取，共得到所有热随机变量的m个组合样本，任意一个热随机变量的值在m个组合样本中出现且仅出现一次；

(5)、针对步骤(4)所得到的每个热随机变量样本组合，先根据样本组合中几何尺寸热随机变量，对热防护系统的瞬态热分析模型进行网格化处理，得到有限元网格，再根据样本组合中其他热随机变量采用有限元分析工具对热防护系统的瞬态热分析模型进行空间离散，共得到m个以时间t为独立变量的线性常微分方程组；

(6)、采用时间积分方法方法求解步骤(5)得到的m个线性常微分方程组，得到m个温度场向量φ(x,y,z,t)，所述温度场向量φ(x,y,z,t)为整个热防护系统空间内所有有限元节点在时间历程内每个计算时刻的温度计算结果；

(7)、从每个温度场向量φ(x,y,z,t)中提取指定区域整个时间历程内的温度特征量，得到指定区域的m个热随机变量组合样本所对应的随机输出温度特征量；

(8)、针对步骤(7)中所得到指定区域的m个热随机变量组合样本及其所对应的随机输出温度值，利用响应面方法建立指定区域温度特征量的高精度多项式响应面模型；

(9)、对步骤(8)中所得到的指定区域温度特征量的高精度多项式响应面模型进行n次蒙特卡洛抽样仿真，从仿真结果中提取n个指定区域随机输出温度特征量，根据n个指定区域随机输出温度特征量计算该指定区域的热防护可靠度，作为热防护系统的可靠度；

(10)、当热防护系统的可靠度大于等于预设门限时，结束；否则，调整热防护系统中热防护结构件的设计参数，重复步骤(2)～步骤(10)，直到热防护系统的可靠度大于等于预设门限。

所述m为正整数，m的取值满足如下条件：

m≥2×(N+1)×(N+2)。

所述蒙特卡洛抽样仿真次数n的取值不低于10⁶。

所述步骤(9)中热防护系统的可靠度的计算公式为：

式中，n为蒙特卡洛抽样次数，Tlim为指定区域随机输出温度特征量许用值，y为指定区域随机输出温度特征量计算结果，N(Tlim>y)表示指定区域温度特征量计算结果低于许用值的次数。

所述步骤(5)中根据几何尺寸进行网格化处理的具体过程为：

初始时，根据几何尺寸随机变量均值，将热防护系统几何模型划分为包含n1个节点、n2个单元的有限元网格，并采用节点号数组NP和单元信息数组E表示，其中，节点号数组NP为n1行1列一维数组；单元信息数组E为n2行9列二维数组，单元信息数组每行包括单元号和对应单元内8个节点编号信息；根据几何尺寸随机变量均值，获取网格所有节点的基准坐标写入坐标数组P，坐标数组P 为n1行3列的二维数组，每行数组元素包括每个节点3个自由度的坐标向量x,y,z；

后续，每次从热随机变量组合样本中得到新的几何尺寸随机变量之后，根据几何尺寸随机变量，计算所有节点坐标相对于基准坐标在空间3个自由度方向的增量{hx,hy,hz}，更新节点坐标数组为P'，P'每行数组元素分别为x+hx,y+hy,z+hz；根据节点号数组NP和节点坐标数组P'生成参数化节点，再根据单元信息数组E生成单元，得到新的有限元网格。

根据步骤(9)所得的指定区域的n次蒙特卡洛仿真结果，能够得到指定区域热防护系统随机输出温度特征值的热随机变量灵敏度系数，具体步骤为：

(9.1)根据n次蒙特卡洛仿真结果计算输出温度特征值的均值y和标准差σ的估计值；

(9.2)、将n次蒙特卡洛仿真的各热随机变量输入与温度特征值输出进行线性回归分析，计算每个热随机变量的相关系数；

(9.3)、根据相关系数计算得到每个随机变量对于输出的灵敏度系数。

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

(1)、本发明提出一种全新的基于概率技术的热防护设计和可靠性评估方法，有效减少热防护系统设计余量，降低热防护系统重量，降低飞行器的功耗；

(2)、本发明先在小样本基础上进行有限元分析、拟合响应面模型，再在响应面模型上进行大样本的蒙特卡罗抽样，大幅度减少了运算量，提高了效率。

(3)、本发明采用时间积分方法方法求解热防护分析模型的线性常微分方程组，节约了解算时间。

(3)、本发明建立了包含几何偏差的热随机变量。突破现有的防护系统参数化建模的技术难题。

(4)、本发明考虑了包括几何尺寸，材料的密度、比热容、热导率、不同结构件之间界面热阻、热防护系统表面辐射率的多个热随机变量对热防护的影响，建立了热分析模型。

附图说明

图1本发明基于概率的热防护设计方法流程图；

图2本发明实施例典型热防护结构模型示意图；

图3本发明实施例典型热防护结构有限元模型；

图4本发明实施例背温时间历程曲线；

图5a本发明实施例输出变量的概率分布曲线；

图5b本发明实施例输出变量的累积概率分布曲线；

图6a本发明实施例背温与厚度H1的散点图(相关系数为-0.034)；

图6b本发明实施例背温与厚度H2的散点图(相关系数为-0.27)；

图6c本发明实施例背温与厚度H3的散点图(相关系数为-0.65)；

图6d本发明实施例背温与与厚度H4的散点图(相关系数为-0.21)；

图7本发明实施例输入参数灵敏度分析结果。

具体实施方式

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

飞行热防护系统处于飞行器最外侧，对其内部结构进行防隔热保护。一般由多个热防护结构件拼接而成，其总的外形尺寸与飞行器的外形尺寸一致。

如图1所示，一种基于概率技术的飞行器热防护系统设计方法，包括如下步骤：

(1)、根据飞行器的外形尺寸和热环境条件，初步确定每个热防护结构的设计参数，所述设计参数包括材料选型、结构形式、几何尺寸，所述热环境条件是指飞行器飞行过程中的气动热流密度；材料选择时要考虑到防热、隔热要求，每一种材料不能超过使用温度极限。结构形式和结构尺寸设计过程中需要考虑到防隔热要求。

(2)、根据热防护系统中每个热防护结构件的设计参数，确定热防护系统的热随机变量及每个热随机变量的概率分布，所述热随机变量包括几何尺寸，材料的密度、比热容、热导率、不同结构件之间界面热阻、热防护系统表面辐射率；在实际应用中，通常假设不确定参数符合给定均值和标准差的正态分布。

(3)、按照每个热随机变量的概率分布规律，将每个热随机变量的取值范围分为m个等概率区间，对于每个热随机变量，在各等概率区间内随机选取一个值，组成一组样本序列，共得到N个热随机变量的样本序列：

{ei1,ei2,ei3,...,eim},i＝1,2,3,....,N。

(4)、从每个热随机变量样本序列中随机抽取一个样本值，形成所有随机变量的一个样本组合，再从每个热随机变量样本序列剩余的数据中随机抽取一个样本值，形成又一个样本组合，以此类推，进行m次随机抽取，共得到所有热随机变量的m个组合样本，任意一个热随机变量的值在m个组合样本中出现且仅出现一次。具体为：

(4.1)、按照第1个热随机变量x1的概率分布，将x1的取值范围分为m个等概率区间，x1在m个区间取值的概率均为1/m，在x1的每个等概率区间内随机选取一个代表值组成一个热随机变量样本序列，为了求解方程,m为正整数，m的取值满足如下条件：

m≥2×(N+1)×(N+2)

式中，N为热防护系统中热随机变量总数。

如果区间数m>1000，可取区间的中点值代替随机从区间中取值；

(4.2)、按照(4.1)的方式，依次对x2-xN随机取值，共得到N个随机变量样本序列；

(4.3)、从N个随机变量样本序列中每个随机选取一个代表值进行第一次组合；从N个随机变量各剩余(m-1)个代表值中随机选取一个代表值进行第二次组合，接着从N个随机变量剩余m-2个代表值中随机选取一个代表值进行第三次组合，以此类推，完成N次组合。得到m个热随机变量组合样本pi(i＝1,2,3,.....,m)。

该抽样方式的优势有两个方面：一方面有效降低了抽样次数，例如N个随机变量每个有m个表达值，这样就有m^N种随机变量的组合次数，利用该抽样方式仅进行了m次抽样；另一方面保证了样本覆盖率，在整个样本分布空间内均进行了有效抽样。

(5)、针对步骤(4)所得到的每个热随机变量样本组合，先根据样本组合中几何尺寸热随机变量，对热防护系统的瞬态热分析模型进行网格化处理，得到有限元网格，再根据样本组合中其他热随机变量采用现有的有限元分析工具(如ANSYS或者ABAQUS)对热防护系统的瞬态热分析模型进行空间离散，共得到m个以时间t为独立变量的线性常微分方程组。

通用的热结构三维瞬态热分析模型可以描述为

式中，ρ是材料密度(kg/m³)；c是材料比热容(J/kg·K)；t是时间(s)，kx、ky、kz是材料沿物体3个方向(x，y，z方向)的热传导系数(W/m·K)。热防护设计的最主要内容即是求解热结构瞬态温度场，从而确定其温度。求解瞬态温度场问题，是在求解初始条件下满足瞬态热传导方程及边界条件的温度场函数φ。

边界条件是：

式中，nx、ny、nz是边界外法线的方向余弦；是边界上的温度；q是边界上的热流密度(W/m²)；h是对流换热系数(W/m²·K)；φa，在自然对流条件下，φa是外界环境温度；在强迫对流条件下，是边界层的绝热壁温度。

初始条件为

φ＝φ0

式中，φ0为温度的空间分布向量。

利用有限元方法可以得到用来确定n个有限元节点温度的求解方程：

这是一组以时间t为独立变量的线性常微分方程组。式中C是热容矩阵，K是热传导矩阵，C和K都是对称正定阵；P是温度载荷向量，φ是有限元节点温度向量，是有限元节点温度对时间的导数向量。

矩阵K、C和P的元素由有限元单元的相应矩阵元素集成，即

式中，是有限元单元对热传导矩阵的贡献，是有限元单元热交换边界对热传导矩阵的修正，是有限元单元对热容矩阵的贡献，是有限元单元给定热流边界的温度载荷，是有限元单元的对流换热边界的温度载荷。

至此，已将时间域和空间域的偏微分方程问题在空间域内离散为个n节点温度φi(t)的常微分方程的初值问题。

在有限元工具的具体实现过程中，需要将几何尺寸和其他热随机变量进行参数化描述，这是实现概率热分析的基础。对于密度、热导率、比热容等其他热随机变量，通过参数化描述较为简单，不再赘述。对于几何尺寸在有限元分析中的参数化描述，存在一定困难：在有限元分析工具中，几何模型处理、有限元网格划分过程通常需要人机交互式处理，很难将上述建模过程通过程序控制。因此，本发明从有限元网格本身出发，通过几何随机变量直接得到有限元节点坐标，进而生成有限元节点和单元，从而实现了热防护系统几何随机参数的有限元模拟。

根据样本组合中几何尺寸热随机变量，对热防护系统的瞬态热分析模型进行网格化处理的具体方法为：

初始时，根据几何尺寸随机变量均值，将热防护系统几何模型划分为包含n1个节点、n2个单元的有限元网格，并采用节点号数组NP和单元信息数组E表示，其中，节点号数组NP为n1行1列一维数组；单元信息数组E为n2行9列二维数组，单元信息数组每行包括单元号和对应单元内8个节点编号信息；根据几何尺寸随机变量均值，获取网格所有节点的基准坐标写入坐标数组P，坐标数组P为n1行3列的二维数组，每行数组元素包括每个节点3个自由度的坐标向量x,y,z；

后续每次处理时，从热随机变量组合样本中得到新的几何尺寸随机变量之后，根据几何尺寸随机变量，计算所有节点坐标相对于基准坐标在空间3个自由度方向的增量{hx,hy,hz}，更新节点坐标数组为P'，P'每行数组元素分别为x+hx,y+hy,z+hz；根据节点号数组NP和节点坐标数组P'生成参数化节点，再根据单元信息数组E生成单元，得到新的有限元网格。例如，厚度的变化可以根据热防护系统的结构，推算出在其引起的每个节点坐标空间3个自由度方向的增量。

(6)、采用时间积分方法方法求解步骤(5)得到的m个线性常微分方程组，得到m个温度场向量φ(x,y,z,t)，所述温度场向量φ(x,y,z,t)为整个热防护系统空间内所有有限元节点在时间历程内每个计算时刻的温度计算结果。

时间积分方法求解常微分方程的计算过程如下：

将常微分方程变换为如下形式：

其中

H＝-C^-1·K

r＝C^-1·P

由常微分方程的理论，方程积分形式的解答为

将求解域划分为若干时间步长τ＝tk+1-tk，则在初始条件φ0已知的情况下各时刻的解可以由如下的递推关系得到

其中，T＝e^τH。假设非齐次项在时间步长τ内是线性的，即

r(t)＝r0+r1(t-tk)

那么递推关系可以表示为

φk+1＝T[φk+H^-1(r0+H^-1r1)]-H^-1[r0+H^-1r1+r1τ]

从上面的推导过程来看，采用时间积分方法求解方程在每一时间点上得到的是解析解，该方法的最大的特点是算法在步长很大的情况下仍能保持极高的精度，在递推过程中不需求解线性方程组，而只涉及矩阵向量的乘积运算，并且递推矩阵T只需计算一次。

(7)、从每个温度场向量φ(x,y,z,t)中提取指定区域整个时间历程内的温度特征量，得到指定区域的m个热随机变量组合样本所对应的随机输出温度特征量；

温度特征量可以选择计算时间内某一指定位置的温度最大值，如热防护系统的表面温度最大值和背部温度最大值，通常这两个温度值可以表征系统的防热、隔热性能；此外，温度特征量也可以也可以选择某一指定时刻、指定位置的温度值，表征某一种材料或者结构形式的热性能。

(8)、针对步骤(7)中所得到指定区域的m个热随机变量组合样本及其所对应的随机输出温度值，利用响应面方法建立指定区域温度特征量的高精度多项式响应面模型。具体为：

响应面法假设随机输入变量对于随机输出变量的影响可以用数学函数来表达。通常，该函数是一个二次多项式，拟合函数为：

式中，β0是常数项，βi(i＝1,2,...,m)是线性项系数，βij(i＝1,2,...,m,j＝1,2,...,m)是二次项系数，针对热随机变量样本组合pi(i＝1,2,3,.....,m)和对应的m个响应值，即每个热随机变量组合样本所对应随机输出温度值，利用最小二乘法开展回归分析，可以得到这些系数。

(9)、对步骤(8)中所得到的指定区域温度特征量的高精度多项式响应面模型进行n次蒙特卡洛抽样仿真，从仿真结果中提取n个指定区域随机输出温度特征量，根据n个指定区域随机输出温度特征量计算该指定区域的热防护可靠度，作为热防护系统的可靠度；

热防护系统热可靠度的定义为：

式中，PF为失效概率，Tlim为系统所能承受的极限温度，通常作为设计输入给出，y为系统的最高温度，即随机输出变量，f(s)为随机输出变量的概率密度函数。

对于所采用的蒙特卡洛法，可直接从统计输出结果得到：n次蒙特卡洛模拟：

式中，n为蒙特卡洛抽样次数，Tlim为指定区域随机输出温度特征量许用值，通常作为设计输入给出，y为指定区域随机输出温度特征量计算结果，N(Tlim>y)表示指定区域温度特征量计算结果低于许用值的次数。

通过对高精度多项式响应面模型进行n次蒙特卡洛抽样仿真结果还可以计算得到热防护系统温度值的概率分布和热随机变量灵敏度，其中随机输出结果Y，随机输出参数的均值y和标准差σ的估计为：

为了确定各随机输入对随机输出的不确定度贡献，需要将各个输入与输出进行线性回归分析，计算样本相关系数，其中第i个(i＝1,2,3,.....,N)热随机变量的相关系数表达式为：

式中，Xij表示第j次蒙特卡洛抽样仿真中i第个热随机变量的输入值，表示n次蒙特卡洛抽样仿真中第i个热随机变量的平均值，

Yj表示第j次蒙特卡洛抽样仿真输出值，表示n次蒙特卡洛抽样仿真中输出的平均值。

同时，完成灵敏度系数的求解，衡量各设计变量变异性对系统可靠度的影响程度，作为热防护系统可靠性优化设计确定最优解的搜索方向。其中第i个随机变量对于输出的灵敏度系数定义为：

(10)、当热防护系统的可靠度大于等于预设门限时，结束热防护系统设计过程；否则，调整热防护系统的设计参数，重复步骤(2)～步骤(10)，直到热防护系统的可靠度大于等于预设门限，所述预设门限大于99.9937％。

实施例：

以某飞行器辐射式多层热防护系统为例，利用本发明所提供的方法进行概率分析。其中热防护结构模型如图2所示，各层材料的热随机变量如表1所示，有限元模型如图3所示。

表1随机输入变量统计参数

在典型再入气动加热环境下，完成热防护系统确定性瞬态防隔热仿真，设计输入参数均取均值。其中表面温度最高为1183℃。背面温度最高为114℃，背面温度的时间历程，如图4所示。

利用本发明提出的方法进行随机抽样仿真，首先，进行700次随机抽样，利用抽样随机输入数据与背温随机输出结果拟合响应面模型；然后，在响应面模型上进行10⁶次随机抽样，针对结果进行概率特性分析。图5a、图5b为背温的概率分布特性，其均值为114.7℃，标准差为9.9℃。当Tlim为145℃时，系统的可靠度指标计算结果为99.9945％。图6a、图6b、图6c、图6d为厚度随机输入与背温的散点图，由图中可以看出，随机输出变量与H3的相关性最强，且呈现明显的负相关性，也就是说，增大H3可以最高效的降低背温。图7为背温的输入参数灵敏度图，由图中同样可以看出，H3的灵敏度值最大。

为了利用传统设计方法与本发明提出的方法进行比较，针对上述算例，利用传统方法进行分析。Tlim同样为145℃，按照热防护设计经验，给定20℃设计裕度，此时热防护背温限制为125℃，计算得到热防护系统厚度。将两种方法计算得到的厚度和重量列在表2中，在本算例中，利用本方法比传统方法减重8.8％。

表2两种设计方法比较

本发明说明书未公开的技术属本领域公知技术。