烧蚀热响应计算中的不确定性传播分析方法研究

2020-12-07陈江涛胡星志国义军崔鹏程

宇航学报 2020年11期

章超,刘骁,陈江涛,胡星志,国义军,崔鹏程

(1.空气动力学国家重点实验室,绵阳 621000；2.中国空气动力研究与发展中心,绵阳 621000)

0 引言

物理问题的数值模拟存在着大量的不确定性，主要来源于建模和模拟两个过程。建模的不确定性主要是由于人们对现实世界认知的局限性，使得模型不能准确反映真实的物理过程。模拟的不确定性主要体现在数值计算过程中没有准确的描述模型内容，存在输入参数、网格、计算格式、收敛等因素引入的误差和不确定性[1-2]。以再入式飞行器的烧蚀热响应计算为例。该类飞行器以高超声速再入地球时，由于受到严酷的气动加热，往往需要引入烧蚀热防护系统。而提升热防护系统设计的可靠性和精细化程度依赖于对再入过程中系统结构内部温度分布的准确预测[3]。过去五十年来，与烧蚀材料内部热响应计算相关的数学建模和计算方法得到了持续发展，但由于热响应计算所使用的材料物性参数存在大量不确定性，且相互关系不明，综合影响难以确定。为准确获得不确定性影响下的目标变量(例如关键部位温度)的概率分布，需要一套不确定性传播分析方法。

在常用的不确定度量化方法中，蒙特卡洛(Monte Carlo，MC)方法[4]构造简单、易于实现，被广泛使用。但是MC方法收敛速度慢，对计算资源需求巨大。多项式混沌(Polynomial chaos，PC)方法[5]是一种更加高效的不确定度量化方法，通过构建目标变量和不确定性输入参数多项式形式的响应关系来分析目标的不确定性信息，已经有了很多成功的应用[6-8]。但PC展开的项数随着输入变量的维数和展开阶数的增加成几何级数式急剧增加，易导致“维数灾难”现象产生。

在有限的计算资源下，如何通过规模大小可以接受的样本数据构建相对准确的响应模型，分析不确定性参数的传播规律，是各国学者们一直努力的方向。在机器学习和优化设计等领域广泛使用的代理模型为不确定度量化分析提供了新的研究思路。代理模型利用有限的样本，构建系统输出与输入变量的响应关系，最终目的之一是代替耗时的数值模拟或者试验过程，预测新输入变量下的输出。常用的代理模型包括：多项式回归模型，人工神经网络，支持向量机和高斯过程回归(Gauss process regression，GPR)等。本文使用在工程领域应用广泛的GPR模型。GPR模型，也被称为Kriging模型,是由南非采矿工程师Krige[9]于1951年提出，最初应用在地质储量估计，Sacks等[10]将其进一步推广到确定性计算机试验的设计和分析领域。GPR模型的构造与优化是基于贝叶斯理论开展的，对非线性函数具有良好近似能力和独特的误差估计功能。目前GPR模型已经被广泛应用到不确定性建模和参数优化的研究中[11-13]。本文发展了基于GPR模型的不确定度量化方法，并以烧蚀热响应模拟中的不确定性输入参数量化分析为例，验证了该方法的有效性和高效率。

1 烧蚀热响应数值模拟方法

炭化烧蚀材料广泛应用于再入式高超声速飞行器的热防护系统。该类材料在受到外部加热时，其中的树脂会发生热解反应并产生热解气体，热解通常会在一定温度范围内进行，形成热解区。材料热解后留下多孔的炭化层并在表面发生烧蚀，热解气体流经炭化层引射到表面上，与表面烧蚀产物一起对气动加热起阻塞作用，材料的这种热解反应不仅影响表面烧蚀速度，而且还影响材料内部的温度分布。

假设热解型防热材料由三个组分构成,密度分别为ρA,ρB,ρC，体积分数为XA,XB,XC，则固体材料的密度为：

ρs=XAρA+XBρB+XCρC

(1)

每一个组分的热解过程各由一个化学动力学方程式描述：

(2)

假设完全未热解的原始材料的密度为ρv，完全炭化后的材料密度为ρc，引入炭化率：

(3)

防热材料为多孔介质，其孔隙率为：

φ=φv(1-α)+φcα

(4)

热解产生的热解气体在多孔介质中的流动符合达西定理：

(5)

其中，Γ为材料的渗透率，υg为热解气体的黏度，P为气体压力，Vg为热解气体流动的速度。热解气体的流动满足动量守恒方程：

(6)

材料内部的能量输运过程满足能量守恒方程：

(7)

本文中烧蚀热响应计算所采用的求解器是基于开源有限元计算框架MOOSE开发的，以有限元方法作为控制方程的数值求解方法。接下来以能量方程为例对控制方程的离散做简要介绍。有限元方法以变分或加权余量法为基础，通过分块逼近思想，对偏微分方程进行求解。对于能量守恒方程，先将待求变量T写为插值函数与待求系数乘积的形式

(8)

将式(8)代入到能量守恒方程中，并将式中所有项移动到等号同一端，再乘以检验函数后在控制单元内积分，令检验函数等于插值函数，即Wi=Ni，可得到：

(9)

对式(9)中的扩散项运用高斯积分公式：

(10)

并将温度和温度梯度的表达式写成矩阵乘积的形式：

(11)

可得到能量守恒方程的Galerkin弱格式：

(12)

接下来再基于等参变换，构造三维插值函数。通过坐标变换，将任意全局坐标系下的非标准单元变换为局部坐标下的标准单元，为此需要建立坐标转换关系式：

(13)

根据等参变换，待求变量与坐标变换使用同样的插值函数，以变量T为例，在局部坐标系下：

(14)

本文采用拉格朗日插值多项式来构造插值函数。插值函数的形式确定之后，可将全局坐标下的单元有限元方程变换到局部坐标系下。先写出插值函数的变换式：

(15)

同时，全局坐标下的积分项dxdydz变换为：

dxdydz=|J|dξdηdζ

(16)

其中，|J|为雅可比矩阵的行列式。最后可得到单元有限元方程在局部坐标下的表达式：

(17)

对于动量方程和质量守恒方程，可根据同样的方式进行离散。获得单元有限元方程后，将其中的矩阵按照一定的规则放入到总体矩阵，可得到总体方程。对上述控制方程组，最终获得的总体方程为瞬态非线性方程组。本文对其在时间方向采用二阶Euler后差(BDF2)进行时间积分，并结合Newton-GMRES迭代对所得到的大型稀疏矩阵进行求解。上述计算方法及求解器的验证与考核详见文献[3]。

在烧蚀热响应数值模拟中包含了模型构造、不确定性输入参数、数值求解过程等诸多不确定性因素。本文主要分析与材料物性相关的输入参数的不确定性传播。

2 不确定性传播分析方法

针对不确定输入参数的不确定性传播问题，本文发展了基于代理模型的分析流程，包括以下步骤：

(1)不确定性因素的量化；

(2)通过确定性计算生成样本；

(3)构建代理模型；

随着数码技术和互联网应用的进步，摄影已深深的融入社会生活，数码音乐播放器，移动电话等数码产品开始配备摄影功能，拍摄的图片可以无线传播，摄影开始多元化发展。

(4)代理模型校验及可能的模型参数优化调整；

(5)敏感性分析及输出响应统计特性。

流程如图1所示，对模型的精度有要求时可以迭代步骤(2)～(4)。

第2.1节和2.2节将简单介绍使用的Kriging代理模型以及模型校验与敏感性分析方法，结合概率统计方法，得到输出响应的统计特性。

2.1 Kriging代理模型及校验

Kriging代理模型是一种插值模型，插值结果是已知样本的线性加权，

(18)

其中，y(i)为样本点，数量为n，给出加权系数w=[w(1),w(2),…,w(n)]T的表达式，即可得到其他位置的预估值。引入统计学假设：将未知函数看成是某个高斯静态随机过程的具体实现。该静态随机过程定义为

Y(x)=β0+Z(x)

(19)

其中，β0为未知常数，也称为全局趋势模型，Y(x)表示数学期望值；Z(x)表示的是均值为零、方差为σ2(σ2(x)≡σ2,∀x)的静态随机过程。在参数空间的不同位置，这些随机变量存在一定的相关性，通过协方差表示为：

Cov[Z(x),Z′(x)]=σ2R(x,x′)

(20)

式中：R(x,x′)为只与空间距离相关的“相关函数”，距离为零时取值为1，距离无穷大时取值为0，相关性的大小随距离的增大而减小。基于以上假设，Kriging模型寻找最优加权系数w，使得均方差

(21)

最小，并满足无偏差条件

(22)

(23)

式(23)中，R是由相关函数组成的“相关矩阵”，r为未知点与已知样本点之间相关函数组成的“相关矩阵”，F是全1的列矩阵，全局趋势模型的表达式为β0=(FTR-1F)-1FTR-1ys。相关函数的选择及模型参数的训练可以参考文献[14]。

对代理模型的校验可以采用留出法或交叉验证法，模型精度可以采用R2指标来评判，

(24)

如果存在额外的数据量足够的校验数据，可以采用分布比较的方法。假设(x,yv)是校验数据，相应的代理模型响应是(x,ySM)，定义均值指标dE测量它们之间的差异，

dE(yv,ySM)=E(|yv-ySM|)

(25)

其中，E为期望算子。还可以通过将校验数据和代理模型响应值表示成累计概率分布的形式，直观地显示出两者差异。

2.2 输出响应的统计特性和敏感性分析

(26)

式中:z为标准正态分布分位数，s为样本标准偏差，置信水平表示为100(1-α)。α通常取值0.1或者0.05。

敏感性分析是不确定度量化传播的一个重要步骤，其中全局敏感性分析能够衡量模型输入变量的不确定性在其整个分布域内对输出响应的不确定性的影响[15]。本文基于代理模型对输入不确定性参数进行敏感性分析。基本思路如下：

1)重采样生成大量的样本;

2)基于代理模型快速计算输出响应值；

3)基于敏感性分析方法分析每个输入参数对输出响应的贡献。

本文使用了Sobol指标[16]来量化每个输入变量对系统输出方差的贡献程度。Sobol指标通过Sobol分解得到，是一种基于方差分解的全局敏感性分析方法，其优点是考虑了随机输入参数在整个取值范围内对输出响应的贡献以及随机输入参数的交互作用。

许多时候，工程上对于输出响应符合某种特定的分布存在期待，可以引入假设检验中的拟合优度检验来解决这个问题。常用的拟合优度检验方法有χ2检验、柯式检验法等。

至此给出了基于代理模型的不确定度量化的理论基础和主要过程。接下来，将以一维烧蚀热响应问题为例，分析其不确定性的传播过程。

3 烧蚀热响应数值模拟中的不确定性传播分析

3.1 问题描述及样本数据处理

本节选取了一维烧蚀热响应算例，算例来自4th Abaltion Workshop[17]官方网站提供的标准算例。本算例计算所使用的材料，TACOT，是Abaltion Workshop主办方为方便数据对比而提供的一种假想材料，其物性参数的不确定度无法通过试验手段获得，计算条件如图2所示。

图2 计算条件

选择25 mm厚度位置的温度值作为目标输出，引入9个材料物性参数作为不确定性输入，如表1所示。

表1 输入参数及其不确定度范围

图3 确定性数值仿真结果直方图

从图3可以看出，30个样本点的数据呈现较大的散布，分布无明显规律，500个样本点的数据呈现出较强的正态分布特征，对500个样本数据进行拟合优度检验。

由表2得到检验统计量

表2 拟合优度检验数据表

(27)

从结果对比可以发现，30个样本数据和500个样本数据的统计特性相差较大。需要指出的是这里并不说明500个样本数据的结果一定是更准确的，但充分说明MC方法统计结果非常依赖于样本数据的规模。在工程问题中，计算资源有限的情况下需要更加高效的不确定度量化方法。

3.2 代理模型构建及校验

本文基于GPR模型利用30个样本数据构建代理模型，模型基函数为常数，核函数为“Matern 5/2”，在使用代理模型进行不确定性分析之前，首先验证代理模型的准确度，由第2.1节所述，采用两种方式。

1)30样本点自身的交叉验证；采用留一交叉验证，R2=0.99。显示泛化误差很小，模型的准确度很高。

代理模型响应与确定性计算结果的差值分布如图5所示，90%的校验数据与代理模型输出响应的差值在1.5 K以内。图6给出了代理模型响应输出与确定性计算结果的累计概率分布比较，可以看出差别较大的输出结果主要集中在输出温度小于350 K的区域。综合图4～6的校验结果，30个样本生成的GPR代理模型具有很高的精度，能够很好地代替数值模拟过程。

图4 代理模型响应与确定性计算结果比较

图5 代理模型响应与确定性计算结果差值

图6 代理模型响应与确定性计算累计分布比较

3.3 敏感性分析及输入参数优化

采用拉丁超立方[18](Latin hypercube sampling，LHS)采样10000次，采用代理模型进行快速计算，再利用Sobol指标进行敏感性分析。从表3可以看出,对输出响应影响最大的是输入参数1和5，参数2,3影响明显，参数4和7有一定的影响，参数6,8,9对输出响应完全没有影响。

表3 不同物性参数Sobol指标

在本文的一维烧蚀热响应问题中，若将9个参数优化为6个参数(参数1，2，3，4，5，7)，对计算结果没有影响。若进一步优化为4个(参数1，2，3，5)，采用同样的30个样本生成新的代理模型，R2=0.95，在500个确定性计算中的校验中，均值误差为0.03%，期望指标dE=6.00 K,为均值的1.6%。减少不确定性参量的个数，可以减少构建代理模型所需要的样本数量，从而减少确定性计算的耗费。通过对一维烧蚀热响应输入参数的敏感性分析可以得到对此类问题影响最大的输入参数为原始材料导热率、炭化材料导热率、原始材料比热和原始材料密度。需要说明的是，实际上孔隙率对于材料热导率和密度具有关联关系，在本例中不显著可能和响应温度位置、参数不确定度的范围等因素有关。这也说明，代理模型的快速计算分析并不考虑物理本质，是一种高效的替代方法。

3.4 输出响应统计特性

图7给出了输出响应的直方图，与500个样本数据结果也很接近。因此可以认为由30个样本数据构建的代理模型给出的目标变量不确定性统计信息与500个样本数据的结果基本一致，但计算量大大减少。由式(26)可以给出响应均值置信水平为99%的置信区间为(355.64 K，356.08 K)。图8给出了输出响应的概率图，从图中可以看出置信水平为90%时的置信区间为(341.78 K，370.62 K)，区间宽度占比为8%。