功能梯度材料结构热力耦合概率分析方法

刘 健， 张 琨， 李 滕， 辛健强， 姚建尧*,3

(1.重庆大学 航空航天学院，重庆 400044；2.中国航天科技集团公司第一研究院 系统工程业务部，北京 100076；3.重庆大学 非均质材料力学重庆市重点实验室，重庆 400044)

1 引 言

功能梯度材料FGM(functional graded material)是一种由两种或多种成分组合而成的新型复合材料。通过改变材料的相对体积分数、物理状态和几何构型，功能梯度材料可以表现出材料性质的连续空间变化[1]。

在功能梯度材料的设计过程中，最为关键的问题之一是如何描述混合后结构的宏观属性。Wakashima等[2]提出了一个自洽的平均场微观力学W-T(Wakashima-Tsukamoto)模型，Cho等[3]证明了该模型能得到更准确的结果。Kapuria等[4]结合该模型并通过实验验证了该模型的有效性。

已经证明FGM具有优化温度场、降低热应力和提高热阻的能力[5]，但在高温高压的工作环境中，功能梯度材料会出现热应力失配、界面破坏和热腐蚀退化等失效形式。因此，功能梯度材料的热传导和热力耦合问题是学者的研究重点，包括数值算法和解析算法。Chen等[6]用有限元分析了FGM稳态和瞬态热传导问题的灵敏度，给出了响应相对于设计变量的变化率，可以用于FGM的优化问题。Sankar等[7]获得了材料特性随厚度呈指数变化的FGM梁中热应力分布的精确解。

然而，功能梯度材料在设计和加工过程中存在许多不确定性和随机性[8，9]。概率分析方法是一种考虑热不确定性的有效方法。屈强等[10]结合有限元和蒙特卡洛模拟建立了刚性陶瓷瓦热防护系统概率设计方法。Zhang等[11]利用蒙特卡洛模拟研究了典型的多层热防护结构的概率温度响应。辛健强等[12]针对辐射式热防护系统结合蒙特卡洛直接抽样建立了一种概率瞬态热分析方法。实践证明，在热分析中，蒙特卡洛模拟是功能梯度材料考虑参数不确定性的一种准确有效的方法。近年来，代理模型已取代传统的数值模拟广泛应用于概率热传导和热应力分析。Xin等[13]提出一种结合拉丁超立方抽样和响应面法的多层热防护结构概率分析方法。

为确保功能梯度材料在耐热结构中的可靠性和安全性，对功能梯度材料进行热力耦合概率分析是很有必要的。Sugano等[14]用解析法研究了初始温度随机变化的FGM板的温度和热应力统计。Chiba等[15]分析了导热系数和热膨胀系数不确定性对温度和热应力场的影响。Yang等[16]用几个独立的随机变量对FGM板进行了弹性屈曲分析。Ferrante等[17]研究了组成成分的体积分数和孔隙率的随机变化对功能梯度金属/陶瓷板在恒定热边界条件下响应的影响。Chiba等[18]分析了用随机过程表示的热载荷作用下功能梯度板的热弹性问题。然而，上述不确定性分析大多都是针对模型某一特定的随机变量，很少有考虑两个以上的随机变量对FGM热响应的影响。而在实际工程问题中，多种不确定性是同时存在的，同时热响应具有时间历程的性质，其输出是一个多输出问题。因此，建立一种能够考虑多种不确定性及多输出的热力耦合概率设计分析方法对功能梯度材料的设计尤为重要。

当只对单一输出或点进行预测时，代理模型能够有比较好的预测效果。而对于一些时间历程多输出问题，则需要对每个时刻上的输出(设时间维度为N)进行预测。为提高代理模型的效率和精度，本文采用本征正交分解方法[19]对输出变量进行减缩。本征正交分解是一种将高维数据投影到低维空间以实现数据降阶的方法。在构建代理模型时，只需对包含绝大部分能量的前几阶正交基(模态)的投影系数进行预测，从而减少计算量[20]，并能够提高模型的预测精度和稳定性。

本文基于W-T模型建立了材料属性与温度相关的参数化有限元模型，给出了功能梯度材料空间分布的方式和整体性能估计的方法。将随机参数作为输入对模型进行热力耦合分析，以温度和热应力作输出，建立了功能梯度材料的概率分析方法。验证了Kriging代理模型的可行性和准确性，以及POD降阶方法对解决时间历程多输出问题的有效性，并分析了不确定性参数的影响及其灵敏度。最后讨论了该方法的有效性。

2 功能梯度材料确定性建模方法

首先给出材料成分在空间的分布和材料属性随温度变化的函数，并构建了基于W-T模型的金属基陶瓷功能梯度材料模型，用于考虑材料成分的变化和结构整体性能的描述。

2.1 组成成分的空间分布

假设陶瓷体积分数和金属体积分数由幂关系给出,

Vc=(z/t)n,Vm=1-Vc

(1,2)

式中下标m和c分别表示金属和陶瓷成分，z为沿厚度方向的坐标，t为平板厚度，n为幂指数，n≥0。这个幂指数关系描述了平板是如何通过厚度分级的。图1给出了相对于不同幂指数值的陶瓷体积分数Vc。

图1 陶瓷体积分数在板厚方向上的变化

2.2 温度对材料属性的影响

ZrO2:

(3)

Ti-6Al-4V:

(4)

材料参数对温度的敏感程度和变化规律不同，为保证结果的可靠性，准确定义材料与温度相关的函数在高温热弹性研究中显得尤为重要。

2.3 基于W-T模型的性能描述

在Mori-Tanaka模型的基础上，Wakashima等[2]提出的W-T模型广泛应用在功能梯度材料宏观弹性参数预估中。它采用了一种平均场的方法，基于基体中的平均应力的概念，近似估计组成成分的相互作用。导出了有效体积模量K和剪切模量G与热导率k、比热C和热膨胀系数α之间的关系,

(5)

(6)

式中

(7)

(8)

热导率k、比热C和热膨胀系数α的表达式为

(9)

C=C1V1+C2V2

(10)

(11)

对于梯度层的区域来说，某种材料可以作为基体或夹杂物，这就会产生两个不同的结果，这两个结果与HS模型[21]的上下界相似。采取如式(12)所示的基于体积的加权平均来估计，这种平均方法能够在两个界限之间平滑过渡。

K=KW T 1V1+KW T 2V2

(12)

式中KW T 1是将材料1作为基体代入式(5)得到的结果，KW T 2是将材料2作为基体代入式(5)得到的结果。剪切模量G、热导率k、比热C和热膨胀系数α可通过同样的方法计算。

3 功能梯度材料热力耦合概率分析

本文考虑热导率、比热及空间分布等13个不确定性参数对热应力的影响，热力耦合概率分析流程如图2所示，步骤如下。

(1) 识别输入变量及其不确定性，包括温度相关的材料属性和空间分布。

(2) 采用拉丁超立方抽样方法生成独立的输入变量。利用参数化有限元模型和热载荷条件，计算功能梯度平板的瞬态热特性，得到其温度分布。

(3) 将每个时刻温度场以载荷的形式施加到结构上进行瞬态结构分析，输出每个时刻的最大热应力。

(4) 通过输入变量和输出的热应力构建Kriging代理模型。

(5) 利用Kriging代理模型进行蒙特卡洛模拟，得到功能梯度平板的热应力概率特性，并获得每个输入变量对输出结果的灵敏度。

图2 FGM热力耦合概率分析方法流程

analysis method flow

3.1 不确定性有限元模型

本算例中，平板的尺寸大小为30 mm×16 mm，保持底面恒定300 K，上表面施加热流载荷q=300 kW/m2，左右边界保持绝热。首先进行一次确定性热分析，结果表明，功能梯度平板在300 s之前就已达到稳定，上表面最高温度为1281.84 K，如图3所示。在瞬态热分析的基础上进行瞬态结构分析，设置下表面简支的边界条件，输出每个时刻最大von Mises stress。

根据式(3,4)，将材料属性随温度变化的函数中的常数项作为随机参数输入变量，随机参数的离散程度通过变异系数(COV)来控制，其定义为变量标准差与均值之比。根据文献[13]的研究，材料属性的变异系数设置为0.02，另一随机输入变量为材料的空间分布幂指数n，变异系数为0.01，不确定性输入参数列入表1。

3.2 POD降阶的Kriging模型及其精度分析

由于热力耦合分析的输出是时间历程上的热应力，即多输出问题。因此每一个样本的输出为一向量，总输出为一矩阵。样本矩阵可表示为

(13)

式中Pi j为来自第i个训练样本的材料属性和空间分布的第j个随机输入参数，子矩阵Q为训练样本的输出，Qi j为第i个训练样本的第j时刻的最大热应力。m为随机参数的个数，b为训练样本的个数，nt为时间采集步数。因此，矩阵S的每一行代表一个训练样本。

使用传统的Kriging模型预测时间历程上的多输出问题时，无法保证热应力在时域上的连续性。针对这一问题，本文首先通过POD方法提取时间历程上热应力的特征信息，从而降低待预测数据的维度，进而构建Kriging代理模型。

3.2.1 本征正交分解

通过POD将大量相互依赖的变量减少到少量的不相关变量，同时尽可能多地保留原始变量的变化。POD对输出矩阵Q降阶过程如下。

令矩阵V=QT,V的第i列向量vi为第i组参数下热应力在时间历程上的变化。热应力随参数和时间变化的关系可以由式(14)逼近，

(14)

(15)

式中(·,·)为内积运算。Φ为向量集合V张成空间的一组规范正交基，那么这组基中任何一个基都可以用原始向量之间的线性叠加来表示，即

(16)

特征值的大小表征了该特征值对应的POD基所包含矩阵V特征的多少。一般情况下，往往前几个POD基就能包含90%甚至99%以上的广义能。因此，对于本文多输出问题，只需要从中抽取一定数量快照进行本征正交分解，分析少数广义能比重较大的POD基就能得到系统的主导特征。

取Φ的前J基向量按列归一化得到Ψ，降阶系数可通过式(17)求得

A=ΨTV

(17)

因此通过POD降阶后，式(13)的形式可转换为

(18)

式中J≪nt，即输出矩阵经过POD降阶后维度大大降低。不仅减少了后续Kriging模型的构建过程，降阶后的输出矩阵还包含了各时间历程输出之间的关系，这对与时间相关的结果预测很重要。

3.2.1 代理模型构建及精度分析

通过计算得到500个样本构建降阶的Kriging模型，流程如图4所示。

对比POD降阶的Kriging模型和传统的Kriging模型的结果，如图5所示。可以看出，POD降阶的Kriging模型预测时间历程上的多输出问题效果相较传统的要好。传统的Kriging导致了预测结果存在明显的非正常波动。而通过POD降阶后的模型很好地解决了这一问题，降阶后的矩阵提取了输出矩阵的绝大部分信息，考虑了相邻向量之间的关系，从而能够与原始数据吻合得比较好。

图4 降阶的Kriging模型构建流程

图5 降阶的Kriging模型与Kriging及真实值的结果对比

with Kriging and original values

为了评价代理模型的精度，定义归一化均方根误差NRMSE为

(19)

式中s为样本数，PV为预测值，RV为真实值。

如图6所示，当训练样本达到360时，模型开始收敛。通过构建的Kriging模型对整个加载时间历程中的热应力进行预测，结果如图7所示。

3.3 热力耦合概率分析

由于每个样本在时间历程上均有一个多输出结果，在应力-时间图中表现为每一条曲线代表一个样本结果。所有样本各个时刻最大热应力值如图8所示，图中粗曲线为该样本的最大和最小热应力值，星条线为确定性结构分析的结果，其他样本结果均分布在这两条极值曲线中间。可以看出，局部最大热应力的最大值和最小值相差约60 MPa，而稳定时达到了近100 MPa的偏差。由此可知，输入参数的微小随机性能会引起较大的热应力偏差，在工程应用中尤其需要重点关注不确定性对热应力的影响。

图6 NRMSE随着训练样本的变化

图7 预测热应力和原始热应力的对比

图8 500个样本各时刻最大热应力

同样地，为了研究不确定性输入参数对热应力的影响，分别计算每个随机参数与输出的相关系数，如图9所示。可以看出，在整个时间历程中，对最大热应力响应影响最大的是K1和α1，即ZrO2的热导率和热膨胀系数，其中ZrO2的热膨胀系数对热应力的影响是正相关的，而热导率的影响则是负相关的。其次则为E1,α2和n，其他变量都集中在0附近，对结果影响较小。

图9 随机输入变量的灵敏度

4 结 论

根据本文研究，可以得到以下结论。

(1) 本文提出的概率分析方法能够同时考虑多种随机因素对功能梯度材料热响应的影响，并验证了该方法的准确性和有效性，给功能梯度材料的结构设计和灵敏度分析提供一定指导。

(2) 通过灵敏度分析可知，对于各时刻的最大热应力，影响最大的是陶瓷材料的热导率和热膨胀系数，其次为金属材料的热膨胀系数和陶瓷材料的弹性模量以及空间分布幂系数。因此，在功能梯度材料的工程应用中需要重点关注材料参数以及材料空间分布的不确定性对结果的影响。

(3) 基于POD的Kriging模型对多输出的结构响应的预测具有较好的结果。POD降阶法对时间历程上的多输出问题能够大大减少计算量，优化预测结果。