一种航空发动机典型机构多构件失效相关的可靠性仿真方法

1.本发明属于航空发动机可靠性分析与评估技术领域，具体涉及一种航空发动机典型机构多构件失效相关的可靠性仿真方法。


背景技术：

2.航空发动机机构构件数量多、各构件装配关系复杂，在恶劣的服役条件下，由于载荷环境、结构设计公差、材料性能参数、间隙配合公差等因素的影响，极易发生多构件失效，严重影响了机构的可靠性和安全性。此外，由于材料参数、载荷参数和模型参数的多重不确定性导致构件寿命呈现较大的分散性。目前对于可靠性理论方法研究较多，但针对航空发动机典型机构的可靠性仿真方法研究较少。
3.为了获得具有多失效性能函数的最可能失效点(mpp)或最可能失效域(mpfrs)，基于metropolis-hastings(m-h)算法的马尔可夫链蒙特卡罗(mcmc)方法被广泛研究。然而，现有方法大多通过数值分析或工程经验确定马尔科夫链的初始状态，具有较大的主观性。对于复杂的高非线性工程问题，主观确定初始状态是非常困难的。同时，传统的重要抽样方法通常以整个重要抽样样本作为候选样本，较大的样本量严重影响了计算效率。


技术实现要素：

4.为了提高多失效可靠性分析的计算精度和效率，本发明提出了一种确定马尔科夫链初始状态的理论方法和候选样本的缩减策略，通过结合mcmc、主动学习kriging、自适应核密度估计和自适应重要抽样，可以获得一定可靠度下的疲劳寿命和失效风险，用于指导航空发动机机构的可靠性优化设计。因此，为克服现有技术的不足，本发明提供一种航空发动机典型机构多构件失效相关的可靠性仿真方法，提高可靠性分析的计算精度和计算效率。
5.本发明采用的技术方案是：一种航空发动机典型机构多构件失效相关的可靠性仿真方法，包括：
6.步骤s1、根据航空发动机典型机构多构件失效的功能函数，确定随机变量及分布类型；
7.步骤s2、对每个随机变量在μ
±
σ范围内利用拉丁超立方抽样生成少量样本ns×
n，其中，ns为样本数量，n为随机变量维数；并对样本点进行刚柔耦合动力学协同仿真，分别得到各失效构件的疲劳寿命响应值；其中，μ为均值，σ为标准差；
8.步骤s3、基于步骤s1的功能函数，以随机变量作为输入变量，最小疲劳寿命作为输出响应，构建多失效的初始代理模型；
9.步骤s4、对每个随机变量进行蒙特卡洛随机抽样，生成n
×
n个随机样本，n为蒙特卡洛样本数量，n为随机变量维数；以n个样本为候选样本，识别关键样本点并更新初始代理模型；
10.步骤s5、利用步骤s4所更新的代理模型在n个候选样本中选择失效域的样本，并将
其作为马尔科夫链的初始状态；
11.步骤s6、基于马尔科夫链初始状态，利用改进的m-h算法确定失效域中的更多样本，进而识别功能函数中所有的最可能失效域；
12.步骤s7、基于步骤s6所确定的可能失效域中的样本建立自适应核密度估计函数，并将其作为拟最优重要抽样函数，生成n
is
×
n个重要样本，其中n
is
为重要样本数量，n为随机变量维数；
13.步骤s8、为减少步骤s7中候选重要样本的数量，采用自适应截断功能函数可能区域的方法确定自适应缩减样本池；
14.步骤s9、将步骤s8所确定的缩减样本池中的重要样本作为代理模型更新的候选样本，再次识别关键样本点并更新代理模型，根据更新的代理模型进一步确定步骤s8的自适应缩减样本池，如此反复直至代理模型接近真实极限状态函数lsf；
15.步骤s10、利用步骤s9所得到的最终更新的代理模型进行可靠性评估，计算失效概率和变异系数，如果变异系数大于0.05，则增加重要样本数量n
is
以扩大重要样本池，回到步骤s7，执行相应的更新过程，直到再次满足停止条件。
16.进一步地，所述步骤s1中，多失效功能函数g(x)具体表示为：
17.g(x)＝min(g1(x),g2(x),
…
,g
p
(x))
18.其中，g
p
表示第p个失效功能函数；p为功能函数的个数。
19.进一步地，所述步骤s3中，代理模型采用kriging代理模型，初始kriging代理模型g0(x)具体表示为：
20.g0(x)＝f(x)
t
β z(x)
21.其中，x＝[x1,x2,
…
,xm]表示输入变量；m为输入参数的维数；f
t
(x)＝[f1(x),f2(x),
…
,fm(x)]为回归基函数；β＝[β1,β2,
…
,βm]为回归系数；z(x)为局部随机偏差，服从高斯分布(0,σ2)。
[0022]
进一步地，所述步骤s4中，第k次更新的代理模型gk(x)具体表示为：
[0023][0024]
其中，d0为初始训练样本集；dk为第k次增加的训练样本集；g0为初始代理模型；gk为第k次更新的代理模型；f是单位列向量；f
t
(x)＝[f1(x),f2(x),
…
,fm(x)]为回归基函数；为回归系数；g为函数真实响应值；r表示预测点与候选样本点之间的相关性，r
t
为r的转置向量；r表示候选样本点集合中任意两个候选样本点的相关性。
[0025]
进一步地，所述步骤s6中，m-h算法为：
[0026][0027]
其中，r为候选状态x
*
的条件概率密度分布与马尔科夫链前一状态的条件概率密度分布的比值；h(x)为最优重要抽样密度函数；f(x)为输入变量x的联合概率密度函数；if(
·
)为失效指示函数，如果g(x)≤0，值为1，否则值为0；x
(j-1)
为第(j-1)个马尔科夫链状态；x
(j)
为第j个马尔科夫链状态,j＝1,2,...,m；x
*
为备选状态；如果r≥1，x
(j)
＝x
*
；如果r《1,以r的概率接受x
*
为第j个状态x
(j)
，以1-r的概率接受x
(j-1)
为第j个状态x
(j)
。
[0028]
进一步地，所述步骤s7中，自适应核密度函数具体表示为：
[0029][0030]
其中，
[0031][0032][0033]
其中，ω为窗宽参数，λj为局部带宽系数；m为构建密度函数的样本数；n为变量空间的维数；md为不同样本的个数(md≤m)；f(x)为输入变量x的联合概率密度函数；k(
·
)为核概率密度函数；x＝{x
(1)
,x
(2)
,...,x
(m)
}为样本点集；x
(j)
和x
(k)
分别为x的第j个和第k个变量；c为样本点集x＝{x
(1)
,x
(2)
,...,x
(m)
}的协方差矩阵，具体表示为：
[0034][0035]
其中，是样本的平均值。
[0036]
进一步地，步骤s10所述的失效概率和变异系数的计算公式具体为：
[0037][0038]
其中，为失效概率；xi为第i个样本点；n为候选样本数量；f
x
(x)为x的联合概率密度函数；h
x
(x)为重要抽样密度函数；if(
·
)为失效指示函数，如果功能函数g(x)≤0，值为1，否则值为0。
[0039]
所述失效概率的估计方差和变异系数计算公式为：
[0040][0041][0042]
本发明的有益效果：
[0043]
本发明通过结合主动学习kriging(ak)、马尔可夫链蒙特卡罗(mcmc)、自适应核密度估计和重要抽样的优点，提出了一种航空发动机典型机构多构件失效相关的可靠性仿真方法，以解决航空发动机典型机构可靠性分析中的高非线性、多失效区域、失效相关性和小失效概率问题。该方法首先采用ak模型确定马尔可夫链的初始状态，识别所有的最可能失
效域；然后利用失效域中的样本构造自适应核密度估计函数，并生成一系列重要样本；进一步采用自适应缩减样本池策略，选择缩减样本池中的重要样本作为候选重要样本，通过学习函数识别关键样本点继续更新kriging模型和候选样本池；最后基于最终更新的代理模型，计算失效概率。该方法可以通过少量的实际性能函数调用次数获得较高精度的可靠度结果。与目前的相关方法相比，该方法具有三个计算优势：(1)ak模型指导的mcmc程序可以精确地确定马尔可夫链的初始状态，并识别出所有的最可能失效区域。(2)自适应缩减样本池策略使代理模型在更新时仅需要较少的候选样本和较小的候选区域，极大提高了计算效率。(3)利用高精度的代理模型技术代替真实的有限元或动力学仿真模拟，极大降低了计算成本，提高了计算效率。结合代理模型技术，本发明所提出的多构件失效相关的可靠性仿真方法适用于具有高度非线性、多失效区域、失效相关性和小失效概率的复杂工程可靠性分析问题，具有较强的工程应用价值。
附图说明
[0044]
图1是本发明的一种航空发动机典型机构多构件失效相关的可靠性仿真方法。
具体实施方式
[0045]
为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。此外，下面所描述的本发明各个实施方式中所涉及到的技术特征只要彼此之间未构成冲突就可以相互组合。
[0046]
如图1所示，本发明一种航空发动机典型机构多构件失效相关的可靠性仿真方法，以静叶调节机构为例对发明内容进行阐释，包括以下步骤：
[0047]
步骤s1、根据静叶调节机构多构件失效的功能函数，确定随机变量及分布类型。
[0048]
静叶调节机构在工作过程中的多构件疲劳失效是其最主要的失效模式，将多构件疲劳失效作为功能函数。通过灵敏度分析，确定影响多构件功能函数的随机变量，主要包括材料参数(弹性模量e，材料密度ρ，泊松比u)、温度参数(气体温度t，导热系数λ，膨胀系数α)、载荷参数(气动力矩τ，摩擦系数mu)和寿命模型参数(常数a、b、s0)。各随机变量分布类型通过数据收集和假设检验得到。为简化计算复杂度，本实施例考虑了两个摇臂构件的疲劳失效，其功能函数可表示为：
[0049][0050]
s.t.g1(x)＝[n]-n1(x)
[0051]
g2(x)＝[n]-n2(x)
[0052]
其中，x＝{t,λ,α,τ,e,ρ,u,mu,γ,β,s0}；n1(x),n2(x)为两个构件的疲劳寿命；[n]为许用寿命；g1，g2分别表示第1个构件和第2个构件的失效功能函数。
[0053]
步骤s2、对步骤s1的11个关于材料参数的随机变量依其分散性在μ
±
σ(其中μ为均值，σ为标准差)范围内利用拉丁超立方抽样生成少量随机样本ns×
n，其中，ns为样本数量，n为随机变量维数；并根据整个机构的刚柔耦合动力学协同仿真得到危险构件的铰接力极值
载荷，考虑温度对机构的影响，将危险构件作为柔性体，加载极值载荷和热-固耦合分析，基于三参数应力寿命预测模型，获得各构件的疲劳寿命。
[0054]
三参数应力寿命预测模型为：
[0055]
lgn＝lga-blg(s
t-s0)
[0056]
式中，s
t
为加载应力幅值；n为疲劳寿命；a、b、s0为材料常数，可以由材料疲劳数据得到。
[0057]
步骤s3、基于步骤s1的功能函数以材料参数的随机变量作为输入变量，两个构件的最小疲劳寿命作为输出响应，并将步骤s2所述的ns个样本所对应的输入变量和输出响应组成的样本集作为初始训练样本集，构建两个构件失效的初始代理模型。
[0058]
代理模型采用kriging模型，初始kriging模型g0(x)表示为：
[0059]
g0(x)＝f(x)
t
β z(x)
[0060]
其中，x＝[x1,x2,
…
,xm]表示输入变量；m为输入参数的维数；f
t
(x)＝[f1(x),f2(x),
…
,fm(x)]为回归基函数；β＝[β1,β2,
…
,βm]为回归系数；z(x)为局部随机偏差，服从高斯分布(0,σ2)。
[0061]
步骤s4、对每个随机变量进行蒙特卡洛随机抽样，生成n
×
n个随机样本，n为蒙特卡洛样本数量，n为随机变量维数；以n个样本为候选样本，识别危险样本点并更新初始代理模型。具体方法为：估计所有候选样本的功能函数预测值基于预测值和标准差计算样本的u学习函数值；若最小u学习函数值小于2，则将该最小u学习函数对应的样本点作为最佳样本，进一步计算最佳样本点下两个构件的刚柔耦合协同仿真疲劳寿命，选择最小寿命作为疲劳寿命响应值；将最佳样本及其疲劳寿命响应值加入到初始训练样本集中更新代理模型。
[0062]
第k次更新的代理模型gk(x)表示为：
[0063][0064]
其中，d0为初始训练样本集；dk为第k次增加的训练样本集；g0为初始代理模型；gk为第k次更新的代理模型；f是单位列向量；f
t
(x)＝[f1(x),f2(x),
…
,fm(x)]为回归基函数；为回归系数；g为函数真实响应值；r表示预测点与候选样本点之间的相关性，r
t
为r的转置向量；r表示候选样本点集合中任意两个候选样本点的相关性。
[0065]
学习函数u为：
[0066][0067]
其中，为功能函数的预测值；为功能函数所对应的标准差。
[0068]
步骤s5、利用步骤s4所更新的代理模型在n个候选样本中选择失效域的样本，即g(xi)《0(i＝1,2,...,n)的样本，并将其作为马尔科夫链的初始状态。
[0069]
步骤s6、基于步骤s5所确定的马尔科夫链初始状态，利用改进的m-h算法确定失效域中的更多样本{x
(1)
,x
(2)
,...,x
(m)
}，m为失效域中样本数，从而获得功能函数中所有的最可能失效域。
[0070]
m-h算法表示为：
[0071][0072]
其中，r为候选状态x
*
的条件概率密度分布与马尔科夫链前一状态的条件概率密度分布的比值；h(x)为最优重要抽样密度函数；f
x
(x)为输入变量x的联合概率密度函数；if(
·
)为失效指示函数，如果g(x)≤0，值为1，否则值为0；x
(j-1)
为第(j-1)个马尔科夫链状态；x
(j)
为第j个马尔科夫链状态,j＝1,2,...,m；x
*
为备选状态；如果r≥1，x
(j)
＝x
*
；如果r《1,以r的概率接受x
*
为第j个状态x
(j)
，以1-r的概率接受x
(j-1)
为第j个状态x
(j)
。
[0073]
步骤s7、基于步骤s6所确定的失效域中的样本建立自适应核密度估计函数，并将其作为拟最优重要抽样函数，生成n
is
×
n个重要样本，其中n
is
为重要样本数量，n为随机变量维数。
[0074]
自适应核密度函数为：
[0075][0076]
其中：
[0077][0078][0079]
其中，ω为窗宽参数，λj为局部带宽系数；m为构建密度函数的样本数；n为变量空间的维数；md为不同样本的个数(md≤m)；f(x)为输入变量x的联合概率密度函数；k(
·
)为核概率密度函数；x＝{x
(1)
,x
(2)
,...,x
(m)
}为样本点集；x
(j)
和x
(k)
分别为x的第j个和第k个变量；c为样本点集x＝{x
(1)
,x
(2)
,...,x
(m)
}的协方差矩阵，具体表示为：
[0080][0081]
其中，是样本的平均值。
[0082]
步骤s8、为减少步骤s7中候选重要样本的数量，提高计算效率，采用自适应截断功能函数可能区域的方法确定自适应缩减样本池。
[0083]
多失效功能函数g(x)＝min(g1(x),g2(x))在输入(x)空间和g空间的可能区域和可定义为：
[0084][0085]
[0086]
其中，τ是一个2元素向量，其所有元素都等于阈值；ω
g-τ
与ω
g τ
分别表示lsf性能函数的g-τ与g τ区域。
[0087]
步骤s9、将步骤s8所确定的缩减样本池中的重要样本作为代理模型更新的候选样本，再次更新kriging代理模型，根据更新的代理模型进一步确定步骤s8的自适应缩减样本池，如此反复直至多失效功能函数满足u函数停止准则，此时所更新的代理模型接近真实极限状态函数lsf。
[0088]
步骤s10、利用步骤s9所得到的最终更新的代理模型进行可靠性评估，计算失效概率和变异系数，一般认为系统失效概率估计的变异系数在0.05左右为可接受，如果变异系数过高(大于0.05)，则增加重要样本量n
is
以扩大重要样本池，回到步骤s7，执行相应的更新过程，直到再次满足停止条件。
[0089]
失效概率为：
[0090][0091]
其中，为失效概率；xi为第i个样本点；n为候选样本数量；f
x
(x)为x的联合概率密度函数；h
x
(x)为重要抽样密度函数；if(
·
)为失效指示函数，如果功能函数g(x)≤0，值为1，否则值为0。
[0092]
失效概率的估计方差和变异系数为：
[0093][0094][0095]
提供以上实施例仅仅是为了描述本发明的目的，而非要限制本发明的范围。对于本领域的技术人员来说，本发明可以有各种更改和变化。凡是在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的权利要求范围之内。