一种基于蒙特卡罗的高速旋转热障涂层可靠性评估方法

技术领域

本发明属于航空发动机涡轮叶片热障涂层可靠性分析技术领域，具体涉及一种基于蒙特卡罗的高速旋转热障涂层可靠性评估方法。

背景技术

热障涂层(Thermal barrier coatings，TBCs)是由隔热的陶瓷层、抗氧化且增强结合力的粘结层以及镍基高温合金基体组成。此外，在制备和高温服役过程中，由于氧气的扩散并与金属元素的反应，在粘结层和陶瓷层之间形成第四层，称为热生长氧化物层(TGO)。各层成分和界面微观结构的复杂性、热力学性能的差异性，再加上涡轮叶片等热端部件形状的复杂、服役环境的恶劣，造成涂层在没有先兆的情况下发生开裂、剥落，最终导致失效。为了保证TBCs在航空发动机上的安全应用，必须对涂层发生开裂、剥落失效的时间即服役寿命进行准确的预测。

为此，研究者对TBCs的寿命预测进行了大量的探索。He等基于断裂理论提出当能量释放率达到断裂韧性时，就会在涂层的拉伸区域产生裂纹。Busso等基于连续损伤力学建立了TBCs的寿命预测模型。其次，还有学者关注某一特定参量的变化开展TBCs的寿命预测，如裂纹扩展长度、损伤面积等。也有研究人员认为界面氧化是导致TBCs失效的重要因素并以此为基础评估TBCs的服役寿命。这些模型所考虑的损伤参量是单一的，不能全面反映出TBCs多样、复杂的失效模式，也没有考虑TBCs实际服役环境的复杂性与不确定性。因此，各种寿命预测模型的适用范围是非常有限的。

现阶段TBCs的寿命预测并未考虑高速旋转状态，而高速旋转是涡轮叶片TBCs的真实服役工况，也是造成热障涂层剥落的重要因素。同时，材料参数和微观结构的不均匀、服役环境的复杂多变，TBCs的使用寿命充满了不确定性，试图用单一损伤参数的演化来预测出一个确定的服役寿命必然会与实际情况存在很大的误差。因此，实际寿命与目前的寿命预测模型之间仍然存在较大差距。

发明内容

为了克服上述现有技术存在的不足，本发明的目的在于提供一种基于蒙特卡罗的高速旋转热障涂层可靠性评估方法，能够准确预测出燃气热冲击下高速旋转热障涂层的失效概率及失效位置，并获得影响涂层失效的敏感性因子，在服役环境参数不可控的情况下，通过调整涂层制备中可控的敏感性参数来提高热障涂层的可靠性。

为了实现上述目的，本发明采用的技术方案是：

一种基于蒙特卡罗的高速旋转热障涂层可靠性评估方法，包括以下步骤；

步骤一(S1)：分析燃气热冲击下高速旋转热障涂层的失效行为，获得热力化耦合主要破坏模式；

步骤二(S2)：根据所述主要破坏模式建立失效准则，得到功能函数；

步骤三(S3)：确定热力化耦合服役环境、微结构、性能的各参量的随机统计特性；

所述功能函数中存在材料参数的分散性、几何参数的随机性和载荷参数的不确定性，利用随机变量统计学获得每个参数的均值、标准差及分布类型；

步骤四(S4)：利用蒙特卡罗法程序进行热障涂层可靠性计算，将所述(S2)中的功能函数和(S3)中每个随机参数的均值、标准差及分布类型输入到蒙特卡罗模型中，即可得到涡轮叶片上单个节点的失效概率；

步骤五(S5)：通过有限元软件获得热障涂层表面的温度等状态参数分布，将其代入(S4)中的蒙特卡罗循环程序中，即可得到涡轮叶片所有节点的失效概率，最后采用Tecplot可视化软件将该计算结果投影成分布云图；

步骤六(S6)：拟合所述失效概率和各参量的二次函数，通过所述(S4)计算得到随机参数取不同值时的失效概率，然后以随机参数为横坐标、失效概率为纵坐标在origin中作出散点图，再用加权最小二乘法拟合出失效概率和各随机参数均值的二次函数g(x)，最后采用下式计算出随机参数的平均值敏感性因子w

步骤七(S7)：根据敏感性分析，获得高速旋转热障涂层的关键影响因素。

将所述S6计算出来的敏感性因子采用origin绘图软件作成柱状图，根据柱状图的高低得到影响高速旋转热障涂层的关键因素，柱状图越高，影响越大。

所述步骤一具体为：

为获取热力化耦合主要破坏模式，首先在模拟试验平台上对涡轮叶片热障涂层做寿命实验(也称循环实验)，并用高清相机拍摄不同循环次数下的表面形貌，分析实验过程中服役参数(如温度、转速、曲率半径等)对寿命的影响以及获得热障涂层的宏观失效位置；

其次，对不同参数下热障涂层的截面和表面进行扫描电镜(SEM)分析，从微观上具体分析热障涂层的失效行为；所述失效行为包括TGO异常生长、陶瓷层与过渡层之间的热失配应力、柱状晶自身断裂造成的涂层剥落；

最后，根据涡轮叶片热障涂层的服役环境，建立相对应的模型，采用有限元软件如COMSOL、ANSYS进行模拟，从理论和实验两方面验证得出热障涂层的失效模式；所述失效模式包括热循环、界面氧化、高温冲蚀、腐蚀以及陶瓷层破坏失效。

所述步骤二具体为：

所述失效准则的函数公式为功能函数：

Z＝g(X)＝R(X)-S(X) (1)

其中X＝X

根据步骤一对燃气热冲击下高速旋转热障涂层失效模式的分析，当能量释放率G等于陶瓷层的断裂韧性Γ

Z＝Γ

在简单的应变条件，刚性基底上涂层从自由边界开裂时的能量释放率G满足：

式中，h为热障涂层的厚度，E

因此，由(2)(3)式得：

燃气热冲击下，离心力F与角速度ω可以表示为：

F＝mω

ω＝2πn (6)

式中，m为柱状试样的质量，n为转速，R为旋转半径；

燃气热冲击下高速旋转材料失效的临界角速度与温度分布、材料的弹性常数、热膨胀系数以及几何尺寸等有关，表示为：

式中，ω为旋转的角速度，T为涂层表面的温度场分布，α

根据(5)(6)(7)式，得：

则根据离心力来估算剪切应力，即：

式中，S为柱状试样的横截面的面积。

由于是圆柱试样，则横截面的面积S为：

S＝πr

式中，r为圆柱试样的曲率半径。

由(9)(10)得：

将(11)式代入(4)式得功能函数：

若Z>0表示材料或结构处于安全状态，而Z<0则被认为结构处于失效状态。当Z＝0时，结构处于临界状态。

所述功能函数是用于区分热障涂层安全和失效区域的边界；若功能函数大于零，则热障涂层处于安全状态；若功能函数小于零，则热障涂层处于失效状态；若功能函数等于零，则热障涂层处于临界状态；

所述步骤三具体为：

建立了失效准则后，需要确定涂层材料参数和载荷(热力化耦合服役环境、微结构、性能的各参量)的随机统计特征，以涂层厚度为例进行说明；

首先获取热障涂层截面的SEM图，根据标尺等距测量出不同位置处的涂层厚度值，得到大量热障涂层厚度测量值；

其次，把所有的测量结果从小到大排成顺序x

式中Φ

然后计算标准正态分布分位数q

最后，把数对(q

由式(12)的失效准则可以看出，燃气热冲击下高速旋转热障涂层失效的概率分析中有三类变量：材料参数(如断裂韧性、热膨胀系数、泊松比和密度)、外部载荷(如温度、转速)以及几何参数(如涂层厚度、曲率半径)。其中，分布函数类型近似为正态分布、威布尔分布，是因为生产与科学实验中很多随机变量的概率分布都可以近似地用正态分布来描述，威布尔分布是根据最弱环节模型或串联模型得到的，能充分反映材料缺陷和应力集中源对材料疲劳寿命的影响。

所述步骤四具体为：

由步骤二知功能函数为g(X)，则随机变量X的联合概率密度函数为f

结构的失效概率写成：

式中，I[g

因此，根据式(14)，失效概率的预测值为：

I[g

根据上述公式(14-17)建立基于Monte-Carlo法的计算模型，具体计算流程如下：

(1)假设有n个随机变量X＝[X

(2)每组数据代表一个样本，将每个样本值代入功能函数(12)中进行判断，若Z<0，则失效次数N

(3)当模拟次数≥N时，模型计算结束。如果此时的失效域为F＝{x:g(x)≤0}，并且落入失效域的样本数量为N

所述步骤五具体为：

采用ANSYS、COMSOL有限元软件的计算流体动力学对涡轮叶片的服役环境进行模拟，得到叶片表面的温度分布；然后将叶片表面每个节点的温度数据提取出来；

将不同的温度数据和功能函数(12)中其他随机变量输入到步骤四的循环模型中，即可算出所有节点的失效概率；

最后将所有节点的失效概率采用Tecplot等可视化软件做成可靠性分布云图，由此得到了叶片表面失效概率的分布。

所述步骤六具体为：

通过拟合失效概率和各参数的二次函数，计算各参数的敏感性因子；

假设得到了TBCs失效概率P

Pf＝g(x) (18)

则敏感性因子的表达式可如下所示：

其中，x表示某个随机参数，x

用加权最小二乘法拟合出失效概率和各随机参数均值的二次函数，然后根据式(19)计算出各随机参数的平均值敏感性因子。

所述步骤七具体为：

将步骤六计算出来的敏感性因子采用origin等作图软件作成柱状图，每个随机变量的灵敏度用平均值来描述；

随机变量在正方向上的平均灵敏度在正意义上影响失效概率。失效概率将随着正平均灵敏度的增加而增加。相反，如果随机变量的平均灵敏度为负并增加，则失效概率将降低。本发明的有益效果：

本发明从概率统计的角度提出了一种基于蒙特卡罗的高速旋转热障涂层可靠性评估方法，将其应用于热障涂层材料的可靠性分析领域，计算出燃气热冲击下高速旋转热障涂层发生失效的概率，并分析出各随机参数对热障涂层失效的敏感程度，能更加准确的评估热障涂层的可靠性。该方法为航空发动机上的安全应用提供可靠的依据，因此在航空发动机热障涂层可靠性评估领域具有重大的意义。

另外，本发明提出的可靠性评估方法考虑了涡轮叶片热障涂层高速旋转的真实服役工况，还考虑了热障涂层微结构影响和内部力、热、化耦合失效的本质机制。本发明在提高可靠性评价精度的同时，还能保证可靠性评价的效率，为热障涂层的设计和安全应用提供了保障。

附图说明

图1为本发明基于蒙特卡罗的高速旋转热障涂层可靠性评估方法的流程图。

图2为燃气热冲击下高速旋转EB-PVD热障涂层微观失效形貌的SEM图。

图3为燃气热冲击下高速旋转EB-PVD热障涂层的失效机理示意图。

图4为蒙特卡罗计算方法的具体步骤示意图。

图5为不同随机参数值的高速旋转热障涂层失效概率图。

图6为高速旋转涡轮叶片热障涂层表面失效概率分布及实验结果图。

图7为高速旋转涡轮叶片热障涂层的灵敏度分析图。

具体实施方式

下面结合附图对本发明作进一步详细说明。

本发明提供如下技术方案：

如图1所示，一种基于蒙特卡罗的高速旋转热障涂层可靠性评估方法，包括以下几个步骤：

步骤一(S1)：分析燃气热冲击下高速旋转热障涂层的失效行为，获得热力化耦合主要破坏模式。

热障涂层的服役环境极其恶劣，承受着高温高速旋转、燃气热冲击、冲蚀、CMAS腐蚀等热、力、化学及三者耦合的载荷作用。在这些服役环境中，高速旋转涡轮叶片与燃气冲击交互作用是造成热障涂层剥落的关键因素。图2是燃气热冲击下高速旋转EB-PVD热障涂层微观失效形貌的SEM图。高速旋转热障涂层陶瓷层与粘结层的结合界面较为完整，并未发现明显的内氧化现象。热生长氧化物层(TGO)生长缓慢且致密，其厚度最大时仅为1.35μm，这并不构成热障涂层剥落的薄弱环节。此时，陶瓷层内部的柱状晶发生严重破损甚至断裂，陶瓷层在整体上呈现出减薄剥落的趋势，且不同高度位置的断口较为平整。因此，与静态热障涂层相比，高速旋转热障涂层的微观失效位置发生在陶瓷层内部，而不是粘结层和热生长氧化物层的位置。

图3为燃气热冲击下高速旋转EB-PVD热障涂层的失效机理示意图。高速旋转过程的离心载荷F

目前，根据是否考虑高速旋转，寿命实验平台主要分为热力化耦合静态试验模拟装置和动态试验模拟装置。

热力化耦合静态试验模拟装置主要应用在导向叶片等静止件的考核，目前大多数的研究集中在导向叶片，燃气热冲击下的失效模式主要为TGO的生长造成界面开裂。

热力化耦合动态试验模拟装置主要应用在工作叶片的考核，考虑了高速旋转状态，燃气热冲击下高速旋转的失效模式主要为陶瓷层柱状晶的断裂，造成破坏的因素主要为离心载荷，其次是高速气流的冲刷。

步骤二(S2)：根据主要破坏模式建立失效准则，得到功能函数或极限状态函数。

评估材料(或结构)的可靠性需要建立有针对性的失效准则，而失效准则的数学定义对于评估材料或结构的可靠性评估至关重要。这可以概括为一个函数公式，该公式被命名为极限状态函数或功能函数：

Z＝g(X)＝R(X)-S(X) (1)

其中X＝X

根据步骤一对燃气热冲击下高速旋转热障涂层失效模式的分析，当能量释放率G等于陶瓷层的断裂韧性Γ

Z＝Γ

在简单的应变条件，刚性基底上涂层从自由边界开裂时的能量释放率G满足：

式中，h为热障涂层的厚度，E

因此，由(2)(3)式得：

燃气热冲击下，离心力F与角速度ω可以表示为：

F＝mω

ω＝2πn (6)

式中，m为柱状试样的质量，n为转速，R为旋转半径。

燃气热冲击下高速旋转材料失效的临界角速度与温度分布、材料的弹性常数、热膨胀系数以及几何尺寸等有关，可以表示为：

式中，ω为旋转的角速度，T为涂层表面的温度场分布，α

根据(5)(6)(7)式，得：

则可以根据离心力来估算剪切应力，即：

式中，S为柱状试样的横截面的面积。

由于是圆柱试样，则横截面的面积S为：

S＝πr

式中，r为圆柱试样的曲率半径。

由(9)(10)得：

将(11)式代入(4)式得功能函数(极限状态函数)：

若Z>0表示材料或结构处于安全状态，而Z<0则被认为结构处于失效状态。当Z＝0时，结构处于临界状态。

步骤三(S3)：确定热力化耦合服役环境、微结构、性能等各参量的随机统计特性。

由式(12)的失效准则可以看出，燃气热冲击下高速旋转热障涂层失效的概率分析中有三类变量：材料参数(如断裂韧性、热膨胀系数、泊松比和密度)、外部载荷(如温度、转速)以及几何参数(如涂层厚度、曲率半径)。其中，分布函数类型近似为正态分布、威布尔分布，是因为生产与科学实验中很多随机变量的概率分布都可以近似地用正态分布来描述，威布尔分布是根据最弱环节模型或串联模型得到的，能充分反映材料缺陷和应力集中源对材料疲劳寿命的影响。

各参数的随机特性见表1。

表1各参数的统计特性

步骤四(S4)：利用蒙特卡罗法程序进行热障涂层可靠性计算，获得叶片上每点的失效概率。

蒙特卡罗(Monte-Carlo)方法利用每个随机变量的概率分布中的随机样本来计算响应。失效概率是通过失效频率除以样本总数来确定的。具体地，按照f

结构的失效概率可以写成：

式中，I[g

因此，根据式(14)，失效概率的预测值为：

I[g

根据上述Monte-Carlo的原理，建立基于Monte-Carlo法的计算模型。具体计算过程如下：

(1)假设有n个随机变量X＝[X

(2)每组数据代表一个样本。假设失效域为F＝{x:g(x)≤0}，并且落入失效域的样本数量为N

据上述蒙特卡罗方法的原理及运算步骤，其编程的算法流程图如图4所示。按照该算法可以算出热障涂层的失效概率。在温度T

另外，分析了燃气热冲击下高速旋转热障涂层的可靠性，预测了不同参数分布条件下的失效概率以及失效概率与各参数变化的关系。具体而言，预测了断裂韧性Γ

可以看出，随着涂层热膨胀系数、涂层表面温度、转速、涂层厚度和涡轮叶片质量的增加，热障涂层的失效概率随之增加。而热障涂层的失效概率随着断裂韧性、涂层密度和曲率半径的增大而较小。这些结论与实验及文献的结论一致。

步骤五(S5)：采用Tecplot等可视化软件将计算结果做成可靠性分布云图。

为进一步确定涡轮叶片热障涂层的失效位置，并同时考虑热障涂层表面温度分布的不均匀性，需要获得热障涂层表面的温度分布。而热障涂层表面的实际温度分布很难通过实验获得。

因此，采用ANSYS等软件的计算流体动力学对涡轮叶片的服役环境进行模拟，得到叶片表面的温度分布。然后将叶片表面每个节点的温度数据提取出来。

将温度数据作为平均值和其它数据输入到Monte Carlo模型中，即可算出该节点上的失效概率。自动重复图4所示的步骤，换句话说，在程序中加一个循环，即可逐个计算出所有节点的失效概率。

最后将所有节点的失效概率采用Tecplot等可视化软件做成可靠性分布云图，由此获得了叶片表面失效概率的分布，如图6所示。

可以发现，热障涂层左侧面的最上端失效概率值较大，说明涂层在该位置发生了失效，与其实验结果相一致。

步骤六(S6)：拟合失效概率和各参量的二次函数，计算各参量的敏感性因子。

为了得到各参数变化时对可靠性的影响程度指标，我们通过拟合失效概率和各参数的二次函数，来计算各参数的敏感性因子。以此提炼影响热障涂层剥落失效的关键因素。

假设得到了TBCs失效概率P

P

则敏感性因子的表达式可如下所示：

其中，x表示某个随机参数，x

用加权最小二乘法拟合出失效概率和各随机参数均值的二次函数，如图5中的式子所示。然后根据式(18)计算出各随机参数的平均值敏感性因子，如图7中的数字所示。

步骤七(S7)：根据敏感性分析，获得高速旋转热障涂层的关键影响因素。

将步骤六计算出来的敏感性因子采用origin等作图软件作成柱状图，如图7所示。每个随机变量的灵敏度用平均值来描述。

随机变量在正方向上的平均灵敏度在正意义上影响失效概率。失效概率将随着正平均灵敏度的增加而增加。相反，如果随机变量的平均灵敏度为负并增加，则失效概率将降低。

也就是说，随机变量Γ

这表明燃气热冲击下高速旋转TBCs的可靠性与曲率半径成正比，与转速和温度成反比，这与实验结果非常一致。

另外，n和m的灵敏度最大，其次是α

在服役环境参数(如转速n)不可控的情况下，可以通过调整涂层制备中可控的敏感性参数(如涡轮叶片质量m、涂层热膨胀系数α