CN102663176A - 针对高可靠机械产品的主动可靠性分析评价方法 - Google Patents

Abstract

本发明提出一种针对高可靠机械产品的主动可靠性分析评价方法，首先建立机械产品的极限状态函数，利用蒙特卡洛仿真方法产生符合参数分布的样本点，利用Kriging元模型模拟极限状态函数，通过主动更新DOE，提高模拟精度并确定最优抽样半径，构造重要抽样密度函数，由重要抽样密度函数产生随机样本点，不断更新DOE，然后抽取服从卡方分布的半径随机样本点，确定机械产品的失效概率和失效变异系数，由失效变异系数来限定更新是否结束，最后得到机械产品的可靠度。本发明方法具有效率高、健壮性好和仿真精度高的优点，实现了对空间机构在轨磨损寿命与可靠度综合设计分析，对于航天装备复杂机构的可靠性设计技术工程化具有重要意义。

Description

针对高可靠机械产品的主动可靠性分析评价方法

技术领域

本发明是一种针对高可靠机械产品，Kriging(克里金方法)与半径外重要抽样混合的主动可靠性分析评价方法，属于机械结构可靠性技术研究领域。

背景技术

随着机械产品向着高精度、轻型和高效能的方向发展，高可靠的机械产品正成为当今的主流，而目前机械产品可靠性设计、评价和验证面临的首要困难就是数据缺乏和相关高可靠定量设计和验证技术的空白，另外机械产品构成复杂且功能多样，各功能之间相互耦合作用，造成机械产品失效的极限状态函数高度非线性，因此目前高可靠机械产品的可靠性分析难点在于极限状态函数隐式、高度非线性特点以及小失效概率事件。

一次二阶矩(FORM)和二次二阶矩(SORM)是目前分析机械结构可靠性经典的方法，但是这些方法需要计算极限状态函数关于基本变量的偏导数，而通常情况下极限状态函数是隐式，无法直接微分或解析微分。针对隐式极限状态函数，普遍采用蒙特卡洛仿真方法(MCS)，但是MCS方法对于复杂机械产品小概率事件需要通过大量抽样，尤其是协同FEA(有限元)仿真时，计算时间长，工程上难以进行可靠度计算。因此一些减少方差技术如重要度抽样、线抽样、拉丁超立方抽样、径向抽样以及定向抽样方法应用于结构可靠性分析当中，文献1(J Luc Schueremans，Dionys VanGemert.Benefit of splines and neural networks in simulation based structural reliabilityanalysis[J].Structural & Safety，2005，V27(3)：246-261)将蒙特卡洛仿真与元模型(Metamodel)结合分析复杂结构系统的可靠度，并将基于直接抽样、重要度抽样方法与元模型(低阶多项式、样条函数、神经网络)相结合，在结构可靠度计算效率等方面进行了比较，结果表明，该方法有效的提高了仿真效率，但是该方法的缺陷在于需要大量样本训练神经网络和样条函数。

重要抽样方法将抽样密度函数的抽样中心移到验算点时抽样效率最高。而对于隐式极限状态函数，验算点未知，因此需要通过解析的算法确定设计验算点。目前随机响应面法(SRSM)发展成为模拟隐式极限状态函数的重要研究方向，但是对于高度非线性极限状态函数和小失效概率事件，SRSM有一定的局限。而Kriging方法在回归多项式的基础上考虑模型误差的空间相关性，对于形状复杂、非线性程度较高的极限状态函数可靠度计算具有适用性，然而Kriging方法相关参数优化问题至今没有得到很好解决，在实用中基于Kriging的结构可靠度仿真计算精度并不稳定，限制了它在结构可靠性分析领域的进一步发展，因此需要利用Kriging局部预测值与局部方差数据构造主动学习函数，利用学习函数更新DOE(实验设计)，以判断Kriging模拟极限状态函数的精度，减少可靠度计算误差。

综上所述，如何针对高维、非线性隐式极限状态函数以及小概率事件，进行可靠性分析与评价是目前研究难点。

发明内容

本发明针对高维小失效概率事件以及极限状态函数为隐式、高度非线性特点，提出了一种Kriging元模型与半径外重要抽样(Monte Carlo Radius-Outside ImportanceSampling，MCROIS)混合的针对高可靠机械产品的主动可靠性分析评价方法。

本发明提出一种针对高可靠机械产品的主动可靠性分析评价方法，具体分以下步骤：

步骤一：确定机械产品的可靠性重要件和关键失效模式。

步骤二：利用应力-强度干涉模型，建立机械产品的极限状态函数G(x)，x表示极限状态函数中的随机向量。

步骤三：确定极限状态函数中随机向量x中各参数的随机统计特性，包括参数的分布类型和均值。

步骤四：根据步骤三得到的参数的分布类型，利用蒙特卡洛仿真方法，在样本空间Ω中产生服从步骤三得到的参数的分布类型的随机样本点，样本量为N_mcs。

步骤五：从N_mcs个样本点中任意选取N个构造初始的实验设计(Design ofexperiment，DOE)。

步骤六：根据当前的DOE，构造初始Kriging元模型，模拟极限状态函数，得到极限状态函数的预测值i＝1，2，…，N_mcs；x_i表示N_mcs中第i个样本点。

步骤七：首先，通过极限状态函数的预测值

和方差

构造主动学习函数L(x_i)：

$L\left(x_i\right)=(\hat{G}\left(x_i\right)-a)[2\mathrm{\Φ}\left(\frac{a-\hat{G}\left(x_i\right)}{{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}}\right)-\mathrm{\Φ}\left(\frac{(a-\mathrm{\ζ})-\hat{G}\left(x_i\right)}{{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}}\right)-\mathrm{\Φ}\left(\frac{(a+\mathrm{\ζ})-\hat{G}\left(x_i\right)}{{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}}\right)]$

$-{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}[2\mathrm{\φ}\left(\frac{a-\hat{G}\left(x_i\right)}{{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}}\right)-\mathrm{\φ}\left(\frac{(a-\mathrm{\ζ})-\hat{G}\left(x_i\right)}{{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}}\right)-\mathrm{\φ}\left(\frac{(a+\mathrm{\ζ})-\hat{G}\left(x_i\right)}{{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}}\right)]---\left(1\right)$

$+[\mathrm{\Φ}\left(\frac{a-\hat{G}\left(x_i\right)}{{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}}\right)-\mathrm{\Φ}\left(\frac{(a-\mathrm{\ζ})-\hat{G}\left(x_i\right)}{{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}}\right)]$

其中，a表示极限状态函数阈值，Φ(.)为标准正态累积分布函数，φ(.)为标准正态概率密度函数，中间变量

然后，对N_mcs个样本点都计算主动学习函数值，找到其中最大的函数值：max(L(x_i))。

步骤八：判断学习法则max(L(x_i))≤0.001是否成立：若不成立，将max(L(x_i))对应的样本点加入当前的DOE中，即N＝N+1，更新当前DOE中的样本点，然后转步骤六执行；若成立，当前所得到的max(L(x_i))所对应的样本点就是最优样本点x^*，然后进入步骤九执行。

步骤九：将步骤八中所得到最优样本点x^*作为初始设计验算点，将

的样本值带入式(2)确定初始抽样半径r₁，

然后确定抽样区域(r₁，r_k)，其中r_k＝r₁+3或者r_k＝r₁+4；

$r_1={\mathrm{\β}}_1=\sqrt{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\left(x_{1,j}^{*}\right)}^2}---\left(2\right)$

其中，r₁表示由验算点

所构造的抽样半径，

表示设计验算点

的第j个分量，n表示样本点

的维数，β₁表示标准正态空间中设计验算点

到原点的距离。

步骤十：构造半径外重要抽样密度函数。将原始坐标系点x_i＝[x_i，1，x_i，2，…x_i，n]转化为标准正态空间的随机变量x′_i＝[x′_i，1，x′_i，2，…x′_i，n]，构造半径外重要抽样密度函数：

$\left\{\begin{array}{cc}{f}_{\mathrm{tr}}\left(x_i^{\′}\right)=K_{\mathrm{tr}}\·f\left(x_i^{\′}\right),& \left|\right|x_i^{\′}\left|\right|>r_1\\ 0,& \left|\right|x_i^{\′}\left|\right|\≤r_1\end{array}\right.---\left(3\right)$

其中，中间参数

其中

步骤十一：由步骤十的半径外重要抽样密度函数产生一个随机样本点，加入到当前的DOE中，更新DOE，利用Kriging元模型，计算新的设计验算点

步骤十二：首先，构造新的抽样半径：根据设计验算点

的样本值构造抽样半径

然后，判断误差Δr＝r₂-r₁＜ε是否成立，ε为设置的误差精度，若成立，则将当前得到的r₂作为最优的抽样半径

继续执行下一步；若不成立，则转步骤十一执行，直到满足误差精度要求。

步骤十三：抽取半径随机样本点。抽取服从卡方分布χ²(n)的半径随机样本点R，判断样本点R是否落在抽样区域

内，ξ≤3，如是，将该样本点记为R_j，j＝1，2，…，N_MCROIS，N_MCROIS表示要取得半径样本点的总个数。若不是，则重新进行抽取，直到样本点的数量达到N_MCROIS。

步骤十四：抽取落在抽样区域

的随机样本Y′_j。将极限状态函数的随机变量x正态化，即

其中u_x和σ_x分别为随机向量x的均值和标准差，设随机变量的第j个样本x′_j＝[x′_1.j，x′_2.j，…，x′_n.j]，令

进而获得位于抽样区域的随机变量样本Y′_j：Y′_j＝[y′_j，1，y′_j，2，……，y′_j，n]＝R_j·a_j。

将随机样本点Y′_j带入Kriging元模型计算极限状态函数的预测值

并进行判断得到二元函数

的值：如果

否则

步骤十五：根据式(4)确定机械产品的失效概率，根据式(5)确定机械产品的失效变异系数，具体如下：

${\hat{p}}_f=(1-{\mathrm{\χ}}_n^2\left({\mathrm{\β}}^2\right))\frac{1}{N_{\mathrm{MCROIS}}}\underset{j=1}{\overset{N_{\mathrm{MCROIS}}}{\mathrm{\Σ}}}\left\{I\right[\hat{G}\left(Y_j^{\′}\right)\left]\right\}---\left(4\right)$

${\mathrm{COV}}_{{\hat{p}}_f}=\sqrt{\frac{1-{\hat{P}}_f}{{\hat{p}}_f\·N_{\mathrm{MCROIS}}}}---\left(5\right)$

其中， $\mathrm{\β}=r_n^{*};$

然后，设ρ_max为设定的变异系数上限值，判断

是否成立，若成立，则本方法结束，得到机械产品的可靠度，若不成立，转步骤十三执行。

本发明方法的优点和积极效果在于：

1)采用Kriging元模型来模拟隐式极限状态函数，通过主动的学习方法，不断更新DOE，减小Kriging模拟极限状态函数误差，提高了机械产品可靠度仿真精度；

2)在得到最优抽样半径的基础上，构造重要抽样密度函数，使样本点更多落在抽样半径确定的球区域附近，加速失效概率计算的收敛，提高仿真效率，相对于MCS方法，仿真时间缩短，效率提高，适用于工程的实际应用；

3)本发明的方法适用于小概率事件的机械产品可靠性分析评价，能够解决机械产品由于数据缺乏和相关高可靠性小概率事件的分析评价难题。

附图说明

图1为本发明的主动可靠性分析评价方法的整体步骤流程图；

图2为半径外重要抽样的原理图；

图3为基于Kriging的MCROIS确定最优抽样半径的过程示意图；

图4为空间结构锁对接和解锁的过程示意图；

图5为空间结构锁动力学仿真结果图：a为主被动锁钩接触碰撞力随时间变化的曲线示意图；b为锁钩工作面中点位移随时间变化的曲线示意图；

图6为本发明实例中磨损失效的极限状态空间的示意图；

图7为Kriging与RSM模型对比的示意图：(a)为采用RSM模型，锁钩最大接触应力与弹性模量的关系示意图；(b)为采用Kriging元模型，锁钩最大接触应力与弹性模量的关系示意图；(c)为采用RSM模型，锁钩最大接触应力与密度的关系示意图；(d)为采用Kriging元模型，锁钩最大接触应力与密度的关系示意图；

图8为采用MCS仿真得到的可靠性指标与失效概率变异系数的变化曲线图；

图9为采用本发明方法仿真得到的可靠性指标与失效概率变异系数的变化曲线图；

图10为本发明实施例中可靠性参数的重要度示意图。

具体实施方式

下面结合附图和实施例子对本发明做进一步说明。

本发明方法利用Kriging元模型模拟隐式极限状态函数，然后通过主动学习迭代算法，计算最优点(最接近设计验算点的样本点)，更新DOE，提高Kriging元模型的模拟精度，在此基础上，利用Kriging元模型确定最优抽样半径，构造半径外重要抽样密度函数，使样本点更多落在抽样半径确定的球区域附近，加速失效概率计算的收敛，从而有效的解决了高维、小失效概率事件以及隐式、非线性极限状态函数的可靠性分析评价难题。

如图1所示，本发明提出的主动可靠性分析评价方法，包括以下几个步骤：

步骤一：确定机械产品的失效模式。通过系统故障模式分析，如FMECA(FailureMode，Effects and Criticality Analysis，故障模式、影响和危害性分析)、FTA(Fault TreeAnalysis，故障树分析法)，确定机械产品可靠性重要件和关键失效模式。

步骤二：建立机械产品的极限状态函数G(x)，x表示极限状态函数中的随机向量。根据步骤一确定的可靠性重要件和关键失效模式，利用“应力-强度”干涉模型，建立失效模式对应的极限状态函数。

步骤三：确定极限状态函数中随机向量中各参数的随机统计特性。根据产品试验数据和材料数据手册，确定影响机械产品可靠性参数的随机特性，如分布类型，均值和标准差等。

步骤四：产生随机样本点。根据参数的分布类型，利用MCS在样本空间Ω产生服从参数分布类型的样本点N_mcs。

步骤五：建立初始的DOE。在样本空间N_mcs个样本点中任意选出N个样本点构造DOE。DOE会通过下面的主动学习函数逐步更新。样本点N的选取取决于可靠度大小，可靠度越高，需要的样本量就越大，也可以首先预估一下失效概率

和失效概率变异系数

则初始的DOE的样本量N可根据条件：来确定。

步骤六：计算Kriging响应面模型。利用当前所构建的DOE，建立Kriging元模型，模拟隐式极限状态函数，得到极限状态函数的预测值

x_i表示N_mcs中第i个样本点，是属于当前DOE中的样本点。

步骤七：通过主动学习函数识别最优点。通过Kriging元模型模拟得到的极限状态函数预测值

和方差

构造如式(1)所示的主动学习函数，并计算样本空间N_mcs个点的函数值，通过主动学习函数识别最优点x^*，所述的最优点x^*指最接近设计验算点的样本点，并利用学习准则判别是否停止迭代。

$L\left(x_i\right)=(\hat{G}\left(x_i\right)-a)[2\mathrm{\Φ}\left(\frac{a-\hat{G}\left(x_i\right)}{{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}}\right)-\mathrm{\Φ}\left(\frac{(a-\mathrm{\ζ})-\hat{G}\left(x_i\right)}{{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}}\right)-\mathrm{\Φ}\left(\frac{(a+\mathrm{\ζ})-\hat{G}\left(x_i\right)}{{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}}\right)]$

$-{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}[2\mathrm{\φ}\left(\frac{a-\hat{G}\left(x_i\right)}{{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}}\right)-\mathrm{\φ}\left(\frac{(a-\mathrm{\ζ})-\hat{G}\left(x_i\right)}{{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}}\right)-\mathrm{\φ}\left(\frac{(a+\mathrm{\ζ})-\hat{G}\left(x_i\right)}{{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}}\right)]---\left(1\right)$

$+[\mathrm{\Φ}\left(\frac{a-\hat{G}\left(x_i\right)}{{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}}\right)-\mathrm{\Φ}\left(\frac{(a-\mathrm{\ζ})-\hat{G}\left(x_i\right)}{{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}}\right)]$

通常情况下，极限状态函数阈值a设置为0。

找出N_mcs个样本点的主动学习函数值的最大值，即寻找：max(L(x_i))，i＝1，2，...，N_mcs。

步骤八：根据学习判别准则更新DOE。本发明实施例中设定的学习法则为：max(L(x_i))≤0.001，当max(L(x_i))≤0.001不成立时，更新DOE，具体是，将max(L(x_i))对应的样本点加入当前的DOE中，即N＝N+1，更新当前DOE，然后转步骤六执行，利用更新的DOE计算Kriging元模型；当max(L(x_i))≤0.001成立时，当前所得到的max(L(x_i))所对应的样本点就是最优样本点x^*，然后进入步骤九执行。

步骤九：构造初始抽样半径。将步骤八中所得到最优样本点x^*作为初始设计验算点

n表示每个样本点的维数，

表示设计验算点

的第j个分量(j＝1，2，…，n)。利用当前构造的Kriging元模型在初始设计验算点

构造初始抽样半径r₁，具体将当前设计验算点的各参数值带入式(2)得到。

$r_1={\mathrm{\β}}_1=\sqrt{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\left(x_{1,j}^{*}\right)}^2}---\left(2\right)$

其中，β₁表示标准正态空间中设计验算点

到原点的距离。

确定初始抽样区域(r₁，r_k)，其中，r_k＝r₁+3或者r_k＝r₁+4。

步骤十：构造半径外重要抽样密度函数。将原始坐标系点x_i＝[x_i，1，x_i，2，…x_i，n]转化为标准正态空间的随机变量x′_i＝[x′_i，1，x′_i，2，…x′_i，n]，构造如式(3)所示的半径外重要抽样密度函数：

$\left\{\begin{array}{cc}{f}_{\mathrm{tr}}\left(x_i^{\′}\right)=K_{\mathrm{tr}}\·f\left(x_i^{\′}\right),& \left|\right|x_i^{\′}\left|\right|>r_1\\ 0,& \left|\right|x_i^{\′}\left|\right|\≤r_1\end{array}\right.---\left(3\right)$

其中，中间参数

其中

步骤十一：确定新的设计验算点。根据步骤十的半径外重要抽样密度函数生成一个随机样本点，加入到当前的DOE中，更新DOE，利用Kriging元模型，计算新的设计验算点

抽样原理如图2所示，在n维欧式空间Rⁿ中，以原点为中心，半径为β的n维球将空间划分为两部分，半径外区域为失效域，半径内区域为安全域，半径外重要抽样就是要确定设计验算点，及MPP(Most Probable Point)，从而确定最优抽样半径

也就是步骤十三中的

步骤十二：重新构造抽样半径。根据步骤十一得到的设计验算点

构造抽样半径

表示设计验算点

的第j个分量，并判断误差Δr＝r₂-r₁是否满足判别条件：Δr＝r₂-r₁＜ε，其中ε表示误差精度，ε一般取值为0.01，误差精度选取太小，会产生无穷迭代的情况，直接选取过大，会影响仿真精度。若不满判别条件，则转步骤十一继续更新DOE，重新进行判断，直到满足精度要求；若满足判别条件，则将当前得到的r₂作为最优的抽样半径

最优抽样半径确定过程如图3所示，通过迭代，寻找设计验算点MPP，得到最优抽样半径

也就是步骤十三中的

步骤十三：抽取半径随机样本点。根据最优抽样半径

确定抽样区域

其中，ξ≤3，抽取服从卡方分布χ²(n)的半径随机样本点R，并判断所产生的随机样本是否落入抽样区域

若是，则将该样本点记为R_j，j＝1，2，…，N_MCROIS，否则，重新抽取，直到半径随机样本点的数量达到N_MCROIS。

步骤十四：通过Kriging计算样本点函数值。通过Kriging元模型计算落在抽样区域内的样本点Y′_j的极限状态函数的预测值，如果

否则

是一个二元函数，用于在步骤十五中计算变异系数的值。

步骤十五：计算失效概率和变异系数。

当失效概率变异系数过大时，转步骤十三执行，更新步骤十三的随机样本点，即重新生成N_MCROIS个半径随机样本点，重新计算失效概率和变异系数，直到变异系数满足要求，即

成立，其中ρ_max为设定的变异系数上限值，通常情况下ρ_max≤0.1。

实施例

下面通过一个具体的对空间结构锁微动磨损的可靠性分析评价的工程实例，来介绍本发明的主动可靠性分析方法。

1)空间结构锁微动磨损的失效模式以及机理模型。

在两航天器对接和分离过程中，主动锁钩和被动锁钩在接触面上会因出现周期性的小幅相对运动而导致接触面产生磨损微粒，这种微粒难以从接触表面排除，则会在锁钩之间接触面处发生微动磨损。如图4所示，为空间结构锁对接和解锁的过程。

锁钩的微动运动形式比较复杂，其微动形式为切向微动，很多因素都会影响微动磨损，如相对位移振幅、接触压力、环境中的磨粒、润滑特性、循环次数、接触表面的状态、环境条件等。本发明实施例中设主被动锁钩锁紧时，接触面间法向压力为P_H，磨损量W与主被动锁钩接触表面的法向压力P_H以及主被动锁钩的相对滑动位移振幅S成正比，与锁钩材料的硬度H成反比，则t时刻主被动锁钩的磨损量可以表示为：W(t)＝K.P_H.S(t)/H，其中K为磨损比例系数，假设锁钩间接触区域为一个矩形摩擦副，磨损体积(也就是磨损量)为锁钩摩擦副宽度b，长度L和磨损深度Δh的乘积，即：W＝b.L.Δh，进而得到磨损深度Δh为：

Δh＝K.P_H.S(t)/(H.b.L)

单套结构锁N次对接分离的磨损深度为：

Δh(N)＝K.P_H.S.N/(H.b.L)＝K.P.S.N/H，其中P＝P_H/(b.L)，为接触应力。

2)空间结构锁微动模式失效机理的参数进行仿真。

锁钩的法向接触应力P_H和相对位移振幅S是影响空间结构锁微动磨损的两个重要参数，通过单套结构锁的ADAMS(Automatic Dynamic Analysis of MechanicalSystems，机械系统动力学自动分析)仿真来确定上述两个参数。

通过ADAMS动力学仿真来模拟锁钩的碰撞过程并获得接触压力，碰撞接触压力随着时间变化曲线以及主动锁钩工作面中点的水平位移曲线如图5所示。

仿真结果表明，两锁钩从初始位置1.75秒开始发生碰撞，仿真得到的最大碰撞接触压力为30270.97N；锁钩工作面在拉紧过程中，由于主被动钩的咬合，产生了微小的相对水平位移，主被动锁钩之间最大相对位移差为15.076-14.6326＝0.4434mm，则可以确定微动磨损中的最大相对位移振幅为0.4434mm。因此，通过单套结构锁动力学仿真得到了影响锁钩之间微动磨损的法向接触应力P_H和相对位移振幅S两个关键参数。

3)建立空间结构锁微动模式对应的极限状态函数。

由微动磨损机理模型，应用“应力-强度”干涉模型，建立单套结构锁基于微动磨损失效的隐式非线性极限状态函数：

G＝Δh_max-Δh₀，Δh₀＝N.K.f(e，q，u，d，T).S/H

其中，Δh₀表示锁钩实际磨损量(mm)；Δh_max为锁钩的许用磨损量(mm)；f(e，q，u，d，T)为锁钩碰撞法向接触应力(Pa)；S为锁钩之间相对位移振幅(mm)；e，q，u，d，T分别为锁钩材料的弹性模量、泊松比、剪切模量、密度以及工作温度。

当G＞0对应安全区域，G＜0对应失效区域，如图6所示，其中f_{Δhmax，Δh0}表示Δ_hmax和Δh₀的联合概率密度函数。

4)确定极限状态函数中各参数的随机统计特性。

将以上影响微动磨损可靠性的相关参数随机化，假设所有的参数相互独立。根据材料设计手册和实验统计，得到随机变量的均值、变异系数以及分布类型，如表1所示。

表1参数随机统计特性

序号	随机变量	均值	变异系数	分布类型
					1	K	5.6×10-8	0.055	正态分布
2	e	139181MPa	0.042	正态分布
					3	q	0.3	0.032	正态分布
4	d	2.7kg/mm3	0.05	正态分布
					5	u	462810MPa	0.01	正态分布
6	T	50℃	0.02	正态分布

5)利用MCS在样本空间Ω产生服从表1所示参数分布类型的随机样本，样本数量为N_mcs个。从N_mcs个样本中任意选取N个样本构造初始DOE。根据初始DOE建立Kriging元模型，模拟隐式极限状态函数，得到极限状态函数的预测值。

根据具体实施中的步骤六～八，通过主动的学习方法不断更新DOE，使得DOE中包含更多的最优点，部分DOE结果如表2所示。

表2结构锁部分DOE结果

弹性模量.e	泊松比q	剪切模量u	密度d	温度T	响应值f
						104385.825	0.18647	482758.71	2.118	49.57	146.8749
106785.5	0.17198	578512.5	2.072	37.5	150.3368
						109185.17	0.25888	450840.78	3.003	40.95	150.3368
111584.85	0.21905	434881.81	2.351	45.26	156.7841
						113984.52	0.25164	387004.91	3.189	56.47	159.884
116384.2	0.22991	379025.43	2.863	60.78	163.4412
						118783.87	0.20095	355086.98	2.258	58.19	167.0338
121183.54	0.20457	506697.16	2.723	59.05	170.3816
						123583.22	0.25526	498717.67	2.63	39.22	173.3086
125982.89	0.24802	522656.12	2.909	38.36	176.7502
						128382.57	0.21543	410943.36	2.816	50.43	180.4164
130782.24	0.1756	347107.5	2.956	62.5	184.0943
						133181.91	0.1575	426902.33	2.491	44.4	187.6017
135581.59	0.2444	363066.47	2.444	40.09	190.2559
						137981.26	0.20819	474779.22	3.142	57.33	193.9672
140380.94	0.21181	530635.6	2.025	61.64	197.3094

表2中的响应值f表示锁钩最大接触应力。

利用DOE得到的样本生成二次多项式响应面模型与Kriging模型对比如图7所示。图7中，纵坐标f表示锁钩最大接触应力，单位是牛顿，横坐标e表示弹性模量，单位是MPa，横坐标d表示密度，单位是kg/mm³。从对比可以明显看出，锁钩微动磨损失效的极限状态函数非线性程度非常高，RSM表示二次多项式形式的响应面法，采用RSM模型很难真实的反映极限函数的非线性程度，这样会造成计算精度上的误差，而本发明采用Kriging元模型建立一个与已知信息相关的插值模型来模拟结构的响应，通过主动的学习算法使得Kriging元模型对非线性系统的模拟更加精确，特别在后期半径外重要抽样仿真中，能够更加准确的模拟抽样半径，提高仿真效率和计算精度。

6)结构锁微动磨损可靠性仿真分析评价。

空间结构锁微动磨损失效的极限状态函强非线性且为高维、小失效概率事件，简单MCS方法需要大量抽样，仿真效率低，计算量大。

在Kriging模拟非线性极限状态函数基础上，根据实施步骤九～十五，将DOE、Kriging与MCROIS相结合的主动可靠性仿真方法得到最优抽样半径β，在β球区域附近进行抽样，并与简单MCS抽样方法进行仿真对比，通过以下三方面来验证仿真效率、健壮性和计算精度：

(1)仿真效率：通过MCS仿真抽样次数对比，可以验证主动可靠性仿真方法的效率；

(2)健壮性：针对隐式、非线性和多重失效点的极限状态函数失效概率仿真计算，验证主动可靠性仿真方法的健壮性；

(3)计算精度：通过分析仿真结果误差，来验证主动可靠性仿真方法的计算精度。

图8和图9分别为采用MCS方法与采用本发明方法主动可靠性仿真(Kriging&MCROIS)得到可靠性指标和失效概率变异系数随仿真次数变化的曲线图，表3列出两种方法的抽样次数以及结果对比。

表3分析结果对比

方法	抽样次数	失效概率	变异系数	计算误差
					MCS	1000	0.243	0.05627	11.03％
本发明方法	100	0.258	0.05363	10.51％

根据MCS方法和本发明的主动可靠性仿真方法得到的失效概率和变异系数，可以计算给定置信水平下的最大误差：

其中γ为置信水平，设定置信水平为95％，COV_pf表示变异系数，则MCS和主动可靠性仿真的最大计算误差分别为：

$E_{\mathrm{CMC}}^{\mathrm{pf}}=1.96{\mathrm{COV}}_{\mathrm{pf}}\≈11.02\%$

$E_{\mathrm{MCROIS}}^{\mathrm{pf}}=1.96{\mathrm{COV}}_{\mathrm{pf}}\≈10.51\%$

从计算结果对比中可以看出，在置信水平95％下，MCS与本发明主动可靠性仿真计算误差小于12％，相对误差为6.173％。

分析结论：

(1)在结构锁微动磨损失效极限状态函数为隐式、高维以及非线性情况下，本发明的Kriging与MCROIS混合的主动可靠性仿真方法计算得到的微动磨损可靠度结果与MCS接近，在计算误差基本一致的情况下，显然本发明的主动可靠性仿真方法抽样次数要比MCS少，仿真效率大大提高，因此通过上述分析可以验证本文论述主动可靠性仿真方法的高效率、计算精度和健壮性；

(2)由可靠性参数重要度结果，如图10所示，很明显可以看出，锁钩的工作温度和材料密度为影响锁钩微动磨损可靠性的最重要影响因素，因此在设计改进中应重点考虑以上两个设计参数。

(3)上述计算得到空间结构锁对接20次后，结构锁微动磨损可靠度为0.742，同理利用主动可靠性仿真方法分别计算结构锁对接不同对接次数下的磨损可靠度结果。因此可以确定空间结构锁在轨磨损寿命与磨损可靠度之间的关系，从而可以确定其微动磨损可靠寿命。

Claims (5)

1.一种针对高可靠机械产品的主动可靠性分析评价方法，其特征在于，包括如下步骤：

步骤一：确定机械产品的可靠性重要件和关键失效模式；

步骤二：利用应力-强度干涉模型，建立机械产品的极限状态函数G(x)，x表示极限状态函数中的随机向量；

步骤三：确定极限状态函数中随机向量x中各参数的随机统计特性，包括参数的分布类型和均值；

步骤四：利用蒙特卡洛仿真方法，在样本空间Ω产生服从步骤三得到的参数分布类型的随机样本点，样本点的数目为N_mcs；

步骤五：从N_mcs个样本点中任意选取N个样本点构造初始的DOE，DOE表示实验设计；

步骤六：根据当前的DOE建立Kriging元模型，模拟极限状态函数，得到极限状态函数的预测值

i＝1，2，…，N_mcs，x_i表示N_mcs中第i个样本点；

步骤七：首先，通过极限状态函数的预测值

和方差

构造主动学习函数L(x_i)：

$L\left(x_i\right)=(\hat{G}\left(x_i\right)-a)[2\mathrm{\Φ}\left(\frac{a-\hat{G}\left(x_i\right)}{{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}}\right)-\mathrm{\Φ}\left(\frac{(a-\mathrm{\ζ})-\hat{G}\left(x_i\right)}{{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}}\right)-\mathrm{\Φ}\left(\frac{(a+\mathrm{\ζ})-\hat{G}\left(x_i\right)}{{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}}\right)]$

$-{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}[2\mathrm{\φ}\left(\frac{a-\hat{G}\left(x_i\right)}{{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}}\right)-\mathrm{\φ}\left(\frac{(a-\mathrm{\ζ})-\hat{G}\left(x_i\right)}{{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}}\right)-\mathrm{\φ}\left(\frac{(a+\mathrm{\ζ})-\hat{G}\left(x_i\right)}{{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}}\right)]---\left(1\right)$

$+[\mathrm{\Φ}\left(\frac{a-\hat{G}\left(x_i\right)}{{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}}\right)-\mathrm{\Φ}\left(\frac{(a-\mathrm{\ζ})-\hat{G}\left(x_i\right)}{{\mathrm{\σ}}_{\hat{G}\left(x_i\right)}}\right)]$

其中，a表示极限状态函数阈值，Φ(.)为标准正态累积分布函数，φ(.)为标准正态概率密度函数，中间变量

然后，对N_mcs个样本点分别计算主动学习函数值，找到其中最大的函数值：max(L(x_i))，i＝1，2，...，N_mcs；

步骤八：判断学习法则max(L(x_i))≤0.001是否成立：若不成立，将max(L(x_i))对应的样本点中加入当前的DOE中，更新DOE，然后转步骤六执行；若成立，当前所得到的max(L(x_i))所对应的样本点就是最优样本点x^*，进入步骤九执行；

步骤九：将步骤八中所得到最优样本点x^*作为初始设计验算点

将

的样本值带入式(2)确定初始抽样半径r₁，

然后确定抽样区域(r₁，r_k)，其中，r_k＝r₁+3或者r_k＝r₁+4；

$r_1={\mathrm{\β}}_1=\sqrt{\underset{j=1}{\overset{n}{\mathrm{\Σ}}}{\left(x_{1,j}^{*}\right)}^2}---\left(2\right)$

其中，

表示设计验算点

均第j个分量，n表示样本点的维数，β₁表示标准正态空间中设计验算点

到原点的距离；

步骤十：首先，将原始坐标系的样本点x_i＝[x_i，1，x_i，2，…x_i，n]转化为标准正态空间的随机变量x′_i＝[x′_i，1，x′_i，2，…x′_i，n]；然后，构造半径外重要抽样密度函数：

$\left\{\begin{array}{cc}{f}_{\mathrm{tr}}\left(x_i^{\′}\right)=K_{\mathrm{tr}}\·f\left(x_i^{\′}\right),& \left|\right|x_i^{\′}\left|\right|>r_1\\ 0,& \left|\right|x_i^{\′}\left|\right|\≤r_1\end{array}\right.---\left(3\right)$

其中，中间变量

β₁的卡方分布值

步骤十一：由步骤十的半径外重要抽样密度函数产生一个随机样本点，加入到当前的DOE中，更新DOE，利用Kriging元模型，计算新的设计验算点

步骤十二：首先，根据设计验算点

的样本值构造抽样半径

然后，判断误差Δr＝r₂-r₁＜ε是否成立，ε为设置的误差精度，若成立，则将当前得到的r₂作为最优的抽样半径继续执行下一步；若不成立，则转步骤十一执行；

步骤十三：抽取服从卡方分布χ²(n)的半径随机样本点R，判断样本点R是否落在抽样区域

内，ξ≤3，如是，将该样本点记为R_j，j＝1，2，…，N_MCROIS，N_MCROIS表示要取得半径随机样本点的总个数，若不是，则重新进行抽取，直到样本点的数量达到N_MCROIS；

步骤十四：首先，抽取落在抽样区域

的随机样本Y′_j，具体是：

将极限状态函数的随机变量x正态化，即

其中u_x和σ_x分别为随机向量x的均值和标准差，设随机变量的第j个样本x′_j＝[x′_1.j，x′_2.j，…，x′_n.j]，设参数

则随机样本Y′_j＝[y′_j，1，y′_j，2，……，y′_j，n]＝R_j·a_j；

然后将随机样本点Y′_j带入Kriging元模型计算极限状态函数的预测值

并进行判断得到二元函数

的值：如果

否则

步骤十五：首先，根据式(4)确定机械产品的失效概率

根据式(5)确定机械产品的失效变异系数：

${\hat{p}}_f=(1-{\mathrm{\χ}}_n^2\left({\mathrm{\β}}^2\right))\frac{1}{N_{\mathrm{MCROIS}}}\underset{j=1}{\overset{N_{\mathrm{MCROIS}}}{\mathrm{\Σ}}}\left\{I\right[\hat{G}\left(Y_j^{\′}\right)\left]\right\}---\left(4\right)$

${\mathrm{COV}}_{{\hat{p}}_f}=\sqrt{\frac{1-{\hat{P}}_f}{{\hat{p}}_f\·N_{\mathrm{MCROIS}}}}---\left(5\right)$

其中， $\mathrm{\β}=r_h^{*};$

然后，设ρ_max为设定的变异系数上限值，判断是否成立，若成立，则本方法结束，得到机械产品的可靠度，若不成立，转步骤十三执行。

2.根据权利要求1所述的一种针对高可靠机械产品的主动可靠性分析评价方法，其特征在于，步骤七中所述的极限状态函数阈值a设置为0。

3.根据权利要求1所述的一种针对高可靠机械产品的主动可靠性分析评价方法，其特征在于，步骤十二中所述的误差精度ε取值为0.01。

4.根据权利要求1所述的一种针对高可靠机械产品的主动可靠性分析评价方法，其特征在于，步骤十五中所述的变异系数上限值ρ_max≤0.1。

5.根据权利要求1所述的一种针对高可靠性机械产品的主动可靠性仿真分析评价方法，其特征在于，所述的机械产品为空间结构锁，在对空间结构锁微动磨损的可靠性进行分析评价时，建立的极限状态函数为：G＝Δh_max-Δh₀，Δh₀＝N.K.f(e，q，u，d，T).S/H；其中，Δh₀表示锁钩实际磨损量，单位是毫米；Δh_max为锁钩的许用磨损量，单位是毫米；f(e，q，u，d，T)为锁钩碰撞法向接触应力，单位是帕；S为锁钩之间相对位移振幅，单位是毫米；e，q，u，d，T分别为锁钩材料的弹性模量、泊松比、剪切模量、密度以及工作温度；影响空间结构锁微动磨损的可靠性的参数有弹性模量、泊松比、剪切模量、密度、工作温度以及磨损比例系数K。

Images