0 引言
榫连接结构是航空发动机叶片与轮盘连接的主要形式,在航空发动机飞行剖面转速工况下承受巨大的交变离心载荷,使得微动疲劳损伤成为其主要失效机制。而航空发动机全生命周期存在几何、载荷[1]、材料、监测误差等多种不确定性因素[2],这些因素导致发动机在试验、服役过程中性能偏离设计值,带来性能不达标甚至引发故障。因此针对榫连接结构微动疲劳寿命开展不确定性分析具有重要的意义。
蒙特卡洛法基于大数定律,适用于各种极限状态方程形式、变量维数、变量分布特征,广泛应用于航空发动机结构强度可靠性评估[3]。崔海涛等[4]以弹性模量、摩擦系数和寿命预测模型参数为随机变量,通过104次蒙特卡洛模拟,得到了不同可靠度下的燕尾榫结构微动疲劳寿命;同时,这一方法被用于燕尾榫结构高周微动疲劳寿命可靠性计算,然后基于Miner疲劳累积损伤准则,建立了高低周微动疲劳寿命对数正态分布模型[5]。此外,蒙特卡洛法也应用于航空发动机其他关键部件的研究,例如涡轮盘持久寿命分析[6],且蒙特卡洛抽样次数达到106。因此,对于高安全性要求的航空发动机限寿件失效风险评估问题,其失效属于小概率事件,传统蒙特卡洛方法在模拟时存在计算成本高、效率不足的问题。为此,工程上发展了基于抽样的传统蒙特卡洛模拟方法及其改进算法[7],与传统蒙特卡洛法计算时间成本相比,重要性抽样方法能够将时间降低约1个数量级。
目前,国内关于航空发动机榫连接结构微动疲劳寿命不确定性的研究相对较少。为此,本文中针对航空发动机榫连接结构微动疲劳寿命不确定性分析这一工程问题,提出了改进蒙特卡洛抽样方法,其应用贝塔分布抽样和双峰正态分布抽样对样本进行抽样,然后基于重要抽样思想推导得到变量原始分布特征的失效概率。结合非线性算例分析了分布参数、缩放参数和平移参数对改进蒙特卡洛抽样方法失效概率评估性能的影响。进一步地,针对航空发动机风扇叶片榫连接结构微动疲劳寿命不确定性分析问题,采用双峰正态分布抽样完成了失效概率计算及全局灵敏度分析。
1 改进蒙特卡洛抽样
蒙特卡洛法利用计算机生成服从变量空间随机分布的抽样样本,然后进行随机模拟试验,从而根据试验结果获得系统失效概率。设X=(X1,X2,…,Xn)为系统变量参数,其分布概率密度函数记为fX(X),功能函数记为g(x),则系统失效概率Pf可计算为
Pf=
If[g(x)]fX(x)dx=E{If[g(x)]}
(1)
其中If[g(x)]为功能函数g(x)的指示函数,计算如下:
(2)
蒙特卡洛法借助计算机产生服从分布fX(x)的样本xi,然后模拟系统特性,则式(1)系统失效概率Pf的估计值为
(3)
从式(3)可知产生服从变量分布的随机样本为蒙特卡洛法以及改进的数字模拟方法的基础[8]。式(3)所示的简单随机抽样广泛应用于工程可靠性分析。但是在航空发动机零部件参数不确定性分析中,设计点变量通常为对应分布均值[1,10],属于小概率失效问题,简单随机抽样存在效率低、样本量巨大的问题。针对这一问题,研究人员提出了重要抽样法,其通过采用重要抽样密度函数代替变量分布密度函数,使得落入失效域的样本概率增加,从而提高失效概率预估精度和收敛速度[9]。在重要抽样法中,式(1)可计算为
(4)
式(4)中:hX(x)为重要抽样密度函数,此时失效概率估计值为:
(5)
式(5)中:xi不再服从变量分布概率密度函数fX(x),而是根据重要抽样密度函数hX(x)随机生成。
对于不同的抽样方法和次数,开展轮蒙特卡洛抽样试验,功能函数失效概率期望
及其标准差
为
(6)
1.1 贝塔分布抽样
贝塔分布是一组定义在(0,1)区间的连续概率分布,其概率密度函数计算如下
(7)
式中:α、β为贝塔分布的2个参数,B(α,β)计算为
(8)
α、β均小于1时呈现中心概率密度小、边缘概率密度大的特点,均大于1时正相反,而均为1时则转化为均匀分布。
针对简单随机抽样在失效区域样本分布不足问题,将贝塔分布函数构造为重要抽样密度函数。引入缩放参数kB,将贝塔分布区间其扩增到目标变量区间,使抽样样本xi能够囊括失效区域,计算如下
(9)
式中:
为根据贝塔分布生成的样本,σ为抽样样本标准差,此时重要抽样密度函数可写为
(10)
1.2 双峰正态分布抽样
双峰正态分布由2个正态分布叠加得到,其表达式为
(11)
式中:μ1、μ2、σ1、σ2分别为2个正态分布的均值和标准差,r为第一个正态分布所占的比重,本文取σ1=σ2=σ,r=0.5。抽样密度如图1所示。为了增加分布样本在失效区域的抽样概率,引入平移参数kD,使得
图1 双峰正态分布曲线
Fig.1 Bimodal normal distribution curve
(12)
式中:μ为设计点,则式(11)为重要抽样密度函数hX(x)。
1.3 失效概率评估流程
失效概率评估流程包括固定抽样次数和依据收敛准则2种流程。前者按照抽样次数N直接进行随机模拟试验,然后根据试验结果评估失效概率。相比于简单随机抽样,重要抽样法在计算失效概率时考虑了抽样密度函数的影响。
按照收敛准则的失效概率评估流程如图2所示,收敛判据为前后两次预测失效概率的相对误差连续5次小于收敛指标。初始抽样次数为10,样本池样本为空集,按照抽样次数进行蒙特卡洛模拟,计算样本池的预测失效概率,当预测失效概率满足收敛指标时收敛次数加1;当预测失效概率不满足收敛指标时,收敛次数归零,然后按照抽样次数进行蒙特卡洛模拟,并将模拟样本放入样本池中,在样本池数目达到抽样次数的10倍时,则应更改抽样次数为样本池数目,如此往复,直至连续5次满足收敛指标后输出预测失效概率。
图2 依据收敛准则的失效概率评估流程
Fig.2 The failure probability evaluation process with convergence criteria
2 抽样性能分析
2.1 算例1
计算案例采用西北工业大学吕震宙等[8]提出的非线性功能函数,如下所示:
(13)
式中:X1、X2相互独立,且X1服从正态分布N(10,25),X2服从正态分布N(9.9,25),属于小概率失效。
针对上述功能函数,采用简单随机抽样、贝塔分布抽样和双峰正态分布抽样3种抽样方法,设置蒙特卡洛抽样次数分别为[10,102,103,104,105,106,107,108],开展100轮蒙特卡洛试验。
2.1.1 简单随机抽样
按照简单随机抽样方法,不同抽样次数下的功能函数失效概率期望及其标准差如图3所示,预测失效概率
在抽样次数达到106后趋于稳定,失效概率标准差随着抽样次数的增加呈现幂级衰减,因此在抽样次数比较少的情况下得到的预测失效概率波动也比较大。
图3 简单随机抽样预测结果
Fig.3 The estimated results based on simple sampling
2.1.2 贝塔分布抽样
影响贝塔抽样密度函数及其失效概率评估的因素包括分布参数α、β以及缩放参数kB,通过随机模拟,分析这些参数对失效概率预估性能的影响。
1) 贝塔分布参数。设置缩放参数kB为8,对应样本分布范围为(μ-4σ,μ+4σ)。分布参数α=β且分别取值为0.75、1.00、1.25、1.50、1.75。如图4(a)所示,随着抽样次数的增加,不同分布参数的失效概率预估值均趋于稳定。在α=β=1.50时波动迅速减小并稳定在收敛值上,此时样本能够囊括绝大多数失效区域,且中心区域仍占有较多的样本,使得失效概率评估收敛较快。预估失效概率标准差如图4(b)所示,α=β=1.50情况下预估失效概率标准差在不同抽样次数下均处于较低水平。
图4 不同分布参数下贝塔分布抽样预测结果
Fig.4 The estimated results based on beta sampling with different distribution parameters
2) 缩放参数。设置分布参数α=β=1.50,缩放参数kB分别取值6、8、10、12、14,不同抽样次数下的单轮预测失效概率如图5(a)所示。在kB取值为6时,由于模拟样本分布范围为(μ-3σ,μ+3σ),使得超出区域并没有样本,导致预估失效概率值低于准确值。如图5(b)所示,由于缩放参数增加导致其在概率密度较大的失效区域分布的样本量反而减少,使得预测失效概率标准差随缩放参数增加而增大,在取值为10时可以达到较好的收敛效果。
图5 不同缩放参数下贝塔分布抽样预测结果
Fig.5 The estimated results based on beta sampling with different scaling parameters
2.1.3 双峰正态分布抽样
影响双峰正态分布抽样密度函数的因素为平移参数kD。如图6(a)所示,不同平移参数单轮抽样结果均在106时趋于稳定,图6(b)所示的标准差表明kD在取值为2.0时可以得到较好的收敛速度,此时2个正态分布的均值分别为μ-2σ和μ+2σ。预估失效概率的标准差并非随着平移参数kD的增加而单调递增或单调递减,这是因为当平移参数kD较小时其抽样落入失效区域的样本偏少,其预测失效概率标准差与简单随机抽样相似;而当平移参数kD较大时其抽样样本落入失效区域的概率增加,但是落入概率密度较大的失效区域的样本量减少,导致其标准差也会增加。
图6 不同平移参数下双峰正态分布抽样预测结果
Fig.6 The estimated results based on bimodal normal sampling with different translation parameters
2.1.4 失效概率评估性能分析
设置贝塔分布参数α=β=1.50、kB=10,双峰正态分布参数kD=2.0,进一步对比三种抽样方法对失效概率预测性能的影响。1 000次随机抽样下不同抽样方式样本点的分布如图7所示,简单随机抽样只有极少数样本分布于非失效区域(黑色曲线左侧),贝塔分布抽样和双峰正态分布抽样样本分布较为均匀,从而能够提高失效概率预测精度。
图7 不同抽样方式样本分布
Fig.7 The sample distribution of different sampling methods
如图8(a)所示,在抽样次数为10时,由于双峰正态分布和贝塔分布抽样样本点数目在失效和非失效区域分布量级一致,两者给出了与理论值量级一致的预测失效概率。另一方面,图8(b)为预测失效概率标准差相对于简单随机抽样的值,可以看出双峰正态分布失效概率预测标准差大约为简单随机抽样的一半甚至更少,证明其在收敛性能上要优于另外2种抽样方式。
图8 特定参数下不同抽样方法预测性能对比
Fig.8 The performance of different sampling methods with certain parameters
按照图2流程,设置收敛指标δ为5%,进行10轮蒙特卡洛试验,不同抽样方式达到收敛所需的抽样次数如表1所示。取简单随机抽样108次预测失效概率期望值0.581 2%作为理论值,双峰正态分布抽样最为接近理论计算结果,其平均抽样次数不到简单随机抽样的1/10;而且简单随机抽样标准差最大,双峰正态分布抽样最小。
表1 不同抽样方法算例1预测结果对比
Table 1 The estimation of example 1 with different sampling methods
抽样方式平均抽样次数(105)平均预测失效概率/%标准差/%简单随机18.2900.589 90.018 8贝塔分布4.890 00.579 00.011 8双峰正态分布1.390 00.580 90.008 7
2.2 算例2
九盒段结构由64个杆元件和42个板元件构成,材料为铝合金,其2个主要失效模式的功能函数为[8]:
(14)
式中:R68、R77、R68为单元强度,外载荷P均值为120 kg,所有参数均满足正态分布,如表2所示。
表2 算例2计算参数
Table 2 The computation parameters of example 2
参数均值标准差分布R68/kg83.510.02正态分布R77/kg83.510.02正态分布R78/kg83.510.02正态分布P/kg12030正态分布
类似地,设置收敛指标δ为5%,进行10轮蒙特卡洛试验,其中贝塔分布参数α、β设置为1.5,缩放参数kB=10,双峰正态分布平移参数kD=2.0,3种抽样方法计算得到的结果如表3所示,双峰正态分布抽样达到收敛的平均次数为7.19×105,小于简单随机抽样和贝塔分布抽样,验证了改进蒙特卡洛抽样的准确性和高效性。
表3 不同抽样方法算例2预测结果对比
Table 3 The estimation of example 2 with different sampling methods
抽样方式平均抽样次数(105)平均预测失效概率/%标准差/%简单随机39.600 00.239 50.003 6贝塔分布14.300 00.237 90.005 6双峰正态分布7.190 00.237 20.004 2
3 微动疲劳寿命不确定性分析
3.1 有限元模型
参考某型号风扇叶片、轮盘提取榫连接结构,简化为平面应力模型,如图9所示,榫头顶端施加风扇叶片离心载荷产生的压强P,榫槽两侧施加周向约束,榫槽底部施加全约束,设置接触区域摩擦系数为0.3。榫头、榫槽均采用TC4钛合金材料,如表4所示,参数均值来源于《中国航空材料手册》[11]TC4棒材室温下的低周疲劳数据,考虑随机因素为压强P和应变疲劳参数σ′f,假设两者均服从正态分布。
表4 榫连接结构计算参数
Table 4 The computation parameters of tenon structure
参数均值标准差分布压强P/MPa20020正态分布密度/(kg·m-3)4 620——弹性模量/GPa109——泊松比0.34——应变疲劳参数σ′f/MPa1 56420正态分布b-0.07——ε′f/%2.69——c-0.96——
图9 榫连接结构有限元模型
Fig.9 The finite element model of tenon structure
取表4参数均值,计算结果如图10所示,接触区域应力分布与涡轮叶片榫连接结构[12]、枞树形榫连接结构[13]有限元计算结果相近,最大Mises应力均位于接触区边缘。
图10 榫连接结构应力分布
Fig.10 The stress distribution of tenon structure
3.2 微动疲劳寿命计算
SWT损伤参量在钛合金燕尾榫连接结构微动疲劳预测方面具有较好的效果,该参量综合考虑了临界平面上的循环正应变幅εn,a和最大正应力σn,max为
SWT=σn,maxεn,a
(15)
其中临界平面为SWT值最大的平面,对应微动疲劳寿命为:
(16)
提取榫连接结构接触区域的应力和应变信息,计算各个节点的SWT损伤参量,如图11所示,在滑动区域两侧边缘存在2个SWT极值,与榫连接结构试验件计算结果一致[12]。将最大SWT值代入式(16)得到微动疲劳寿命为91 720,满足10 000周循环的寿命要求。
图11 接触区域损伤参量
Fig.11 The damage parameter in the contact zone
3.3 微动疲劳寿命失效概率
考虑拉伸压强P和应变疲劳参数σ′f两个随机因素(表4),基于双峰正态分布抽样进行蒙特卡洛法模拟,计算榫头微动疲劳寿命失效概率。设定功能函数为
g(X)=Nf(X)-104
(17)
双峰正态分布抽样结果如图12所示,拉伸压强P的增加和材料应变疲劳参数σ′f的减少,都将导致榫连接结构微动疲劳寿命下降,甚至失效。对于同一个失效的样本点,相比于增加材料应变疲劳参数σ′f,改变更少的拉伸压强P就能满足寿命要求,如图12两个箭头所示。因此,降低转速波动更能有效降低结构失效概率。将抽样结果代入式(5)计算得到失效概率为0.046 4%。
图12 榫连接结构抽样样本分布
Fig.12 The sample distribution of tenon structure
3.4 失效概率影响因素分析
考虑拉伸压强P和应变疲劳参数σ′f的不确定性对榫连接结构失效概率的影响时,采用Cui等[14]提出的失效概率全局灵敏度指标,其表达式为
(18)
式中:Pf|Xi为变量Xi不变时的条件失效概率。对于本文中提出的改进蒙特卡洛抽样方法,当积分步长dxi取极小值时,可认为变量Xi服从均匀分布
然后按照式(5)计算得到Pf|Xi。
根据式(18)计算得到榫连接结构微动疲劳寿命失效概率影响因素的全局灵敏度指标如表5所示。拉伸压强P的全局灵敏度指标最大,说明其对榫连接结构微动疲劳寿命失效概率影响最为显著。这是因为拉伸压强P增大时将同时增加接触面的摩擦力和正压力,使得循环正应变幅εn,a和最大正应力σn,max增大,导致微动疲劳损伤参量SWT值迅速增加,且在同等标准差下拉伸压强P相对增加值更大,使得榫连接结构微动疲劳寿命迅速下降。因此,拉伸压强P,即转速波动对榫连接结构微动疲劳寿命影响不可忽略。
表5 全局灵敏度指标
Table 5 Global sensitivity index
因素拉伸压强(P)应变疲劳参数(σ′f)全局灵敏度0.045 80.022 8
4 结论
针对蒙特卡洛法简单随机抽样在评估失效概率效率不足这一问题,本文中提出了基于贝塔抽样和双峰正态分布抽样的蒙特卡洛模拟方法,在不引入大量计算复杂度的前提下提高了不确定性评估效率。得到如下结论:
1) 算例分析结果表明:影响贝塔分布抽样预测失效概率性能的因素包括分布参数α、β和缩放参数kB,缩放参数过小时将会导致失效区域分布样本不足。
2) 算例分析结果表明:影响双峰正态分布抽样预测失效概率性能的因素为平移参数kD,平移参数过小时其失效概率预测性能将会趋于简单随机抽样,过大将会导致落在概率密度较大的失效区域的样本量减少。给定算例下双峰正态分布抽样在预测失效概率精度和收敛性方面均表现更好。
3) 基于双峰正态分布抽样完成了风扇榫连接结构微动疲劳寿命不确定性分析,计算得到失效概率为0.046 4%。基于全局灵敏度指标分析表明拉伸压强P对榫连接结构微动疲劳寿命失效概率影响最为显著,在设计中应给予考虑。
[1] 张屹尚,陈健.航空发动机轮盘结构可靠性全局灵敏度分析[J].兵器装备工程学报,2021,42(4):244-249.
ZHANG Yishang,CHEN Jian.Global sensitivity analysis for turbo-engine disk structures with uncertainties based on failure probability[J].Journal of Ordnance Equipment Engineering,2021,42(4):244-249.
[2] 郑新前,王钧莹,黄维娜,等.航空发动机不确定性设计体系探讨[J].航空学报,2022,42(12):27099-27107.
ZHENG Xinqian,WANG Junying,HUANG Weina,et al.Discussion on uncertainty-based design system for aeroengines[J].Acta Aeronautica et Astronautica Sinica,2022,42(12):27099-27107.
[3] GOVINDARAJAN N,KHOSROW R,UDO N.Fatigue life estimation of aero engine mount structure using Monte Carlo simulation[J].International Journal of Fatigue,2016,83(4):53-58.
[4] 崔海涛,汪震,温卫东,等.微动疲劳寿命可靠性分析方法[J].航空动力学报,2009,24(6):1279-1283.
CUI Haitao,WANG Zhen,WEN Weidong,et al.Reliability analysis method of fretting fatigue life[J].Journal of Aerospace Power,2009,24(6):1279-1283.
[5] 汪震.燕尾榫结构微动疲劳寿命可靠性分析研究[D].南京:南京航空航天大学,2008.
WANG Zhen.Research on fretting fatigue life reliability of dovetail joint[D].Nanjing:Nanjing University of Aeronautics and Astronautics,2008.
[6] 牟园伟,陆山.涡轮盘持久及低周疲劳寿命可靠性评估[J].燃气涡轮试验与研究,2015,28(3):13-18.
MOU Yuanwei,LU Shan.Reliability prediction on creep rupture and LCF life for a turbine disk[J].Gas Turbine Experiment and Research,2015,28(3):13-18.
[7] 严宇飞,李方义.基于蒙特卡洛重要抽样法的结构时变可靠性分析[J].重庆理工大学学报(自然科学),2021,35(6):104-112.
YAN Yufei,LI Fangyi.Time-variant reliability analysis of structures based on important sampling Monte Carlo method[J].Journal of Chongqing University of Technology (Natural Science),2021,35(6):104-112.
[8] 吕震宙,宋述芳,李璐祎,等.结构/机构可靠性设计基础[M].西安:西北工业大学出版社,2019.
LYU Zhenzhou,SONG Shufang,LI Luyi,et al.Reliability design basis of structural mechanism[M].Xi’an:Northwestern Polytechnical University Press,2019.
[9] 赵国藩.工程结构可靠性理论与应用[M].大连:大连理工大学出版社,1996.
ZHAO Guofan.Reliability theory and its applications for engineering structures[M].Dalian:Dalian University of Technology Press,1996.
[10] 柯爽.航空发动机涡轮轴的疲劳寿命预测与疲劳可靠性分析[D].成都:电子科技大学,2020.
KE Shuang.Fatigue life prediction and fatigue reliability analysis of aero-engine turboshaft[D].Chengdu:University of Electronic Science and Technology of China,2020.
[11] 《中国航空材料手册》编辑委员会.中国航空材料手册[M].北京:中国标准出版社,2001.
Editorial board of China aeronautical materials handbook.China aeronautical materials handbook[M].Beijing:Standards Press of China,2001.
[12] ZHANG R,MENG X H,SUN K,et al.An investigation of high and room temperature fretting fatigue of DD6-FGH96 dovetail joint in aero-engine:Experiment and numerical analysis[J].International Journal of Fatigue,2022,154(2):106537-106545.
[13] 魏大盛,王延荣.枞树形榫连结构接触应力的有限元分析及建模研究[J].燃气涡轮试验与研究,2010,23(2):5-9.
WEI Dasheng,WANG Yanrong.On the fem modeling and contact analysis of fir tree attachment in blade disk assemblies[J].Gas Turbine Experiment and Research,2010,23(2):5-9.
[14] CUI L J,LU Z Z,ZHAO X P.Moment-independent importance measure of basic random variable and its probability density evolution solution[J].Science China Technological Sciences,2010,53(4):1138-1145.