基于Copula函数的多功能机械系统竞争失效分析方法

刘敬一, 庄新臣, 张玉刚, 庞欢

1.长安大学 汽车学院, 陕西 西安 710064; 2.清华大学 机械工程学院, 北京 100091;3.西北工业大学 航空学院, 陕西 西安 710072

随着机械设计技术的发展和加工制造水平的提高，机械产品通常具有多个功能，与之相关的失效形式也更为多样化，尤其对应用于航空航天、兵器装备和重型运载工具等领域的机械系统而言，恶劣的工作条件和极端的载荷环境更易导致多失效模式的发生。对于多失效模式的机械系统，其各个失效模式之间往往是相互作用的，这种相互作用体现于：①失效相关，多个失效模式之间并非独立，而是因共用组件、同受载荷等原因彼此作用；②相互竞争，某一种失效模式的出现会引起机构工作状态的变化甚至导致机构系统停止工作，进而造成其他失效模式不会再有机会在同等条件下发生[1]。准确描述多失效模式之间相关且竞争的相互作用关系，成为对机械产品进行可靠性评估、寿命预测和维护策略制定的关键。

关于相关竞争失效的理论和建模方法方面，已有大量的文献开展了相关研究。罗格斯大学的Pham教授较早地把竞争失效的概念推广应用于可靠性工程领域，提出了可应用于工程上的竞争失效模型，认为对于受突发失效和退化失效影响的系统，任一种失效模式的出现都会导致系统的破坏，从而使另一种失效模式没有机会再次发生，也就是说2种失效模式以相互竞争的方式综合决定了系统的可靠性[2-5]。随后罗格斯大学Coit教授分别考虑变失效阈值、系统多部件、冲击分类和变退化速率等[6-9]因素，对相关竞争失效的建模方法进行了扩展，以适用于不同工程背景机械系统的研究。同时，实际工程应用中机械系统的相关竞争失效研究也逐渐受到关注。航空航天方面，李伟、王华伟、郭庆和刘晓娟等[10-14]分析了航空发动机的竞争失效，针对可靠性分析、可靠性评估和寿命预测等内容开展研究；车辆工程方面，闫书法等[15]对某型号重型车辆的传动机构进行了竞争失效的建模，并利用试验对所建立的模型进行验证和应用。Yu等[16]提出可用于分析多态退化可修系统的可靠性模型，考虑系统状态可退化的情况对其可靠性表现进行评估，并把相应模型应用于汽车车架结构的疲劳失效分析；船舶及水下运载工具领域，Bocchetti等[17]考虑磨损失效和热应力失效，研究了某型轮船燃气轮机的竞争失效现象，并基于此对该机械系统进行了可靠性评估，鄢伟安[18]则研究了水下鱼雷翼面展开机构的竞争失效；兵器装备领域，王浩伟、蔡忠义、孙丽和逯程等分别研究了导弹贮存[19]、弹翼机构[20]、起爆电容[21]和大型装备[22]的竞争失效现象及维护保障方法。

综合现有研究可知，当前研究更多地关注因载荷过大所导致的结构破坏失效模式，却很少考虑“结构并未破坏但因退化因素累积导致功能无法完成”所造成的机械系统功能失效。在如今的材料科学技术和加工制造水平下，机械系统发生结构破坏失效的比例正逐步降低，更多的失效形式体现于运动卡滞、运动精度不足或构件磨损等形式的退化型功能失效；同时，当前研究往往假设失效模式相互独立或使用相关系数对失效模式之间的相关性进行描述，独立假设并不符合工程实际，而通过相关系数无法获取相关情况下多失效模式的表征量的联合概率密度函数，难以进行可靠性的计算求解。因此，本文建立了面向多功能机械系统的建立竞争失效分析模型，针对处于突发型结构破坏失效和多个退化型功能失效综合作用下的机械系统，研究各个失效模式的表征方法和退化规律；基于Copula函数定量化失效模式之间的相关性，获取多失效模式功能表征量的联合概率密度函数及演化规律，形成可靠性评估方法；最后利用本文方法针对某型飞机舱门锁机械系统进行了竞争失效分析，证明了本文所提方法的正确性和工程实用性。

1 机械系统多功能失效描述

功能失效指的是“机械系统在部件尚未发生明显破坏时因功能无法完成而导致的失效”，机械产品的特定功能可以通过一个正确描述功能现象的指标进行表达，该指标定义为“功能表征量”。常见的功能表征量包括运动精度类、驱动力/阻力类和运动时间类等[23]。随着机械系统本身的复杂化和综合化，通常一个机械产品具有多个功能。对于大型机械系统而言，具有精密度高、作用载荷大、功能复杂等特点，需要进行合理的可靠性建模分析。机械产品的失效模式可以概括为突发型结构破坏失效和功能无法完成导致的退化型功能失效。

理想情况下机械系统的各个部件有足够的能力承受工作载荷。但是由于随机性的存在，工作载荷本身并非定值，而是一个随机变量。且在复杂工作环境下的机械系统，一些偶发因素也会引发额外载荷，例如外界物体对机械系统的碰撞，机械系统本身的振动等，都会造成量值较大的载荷作用于各个部件，这种明显大于工作载荷的偶发载荷，可看作是一种冲击，当冲击的量值超过部件的承受能力时，部件发生断裂，进而引起机械系统突发失效。另一方面，部件参数的退化会引起功能表征量的退化，进而发生功能失效，例如磨损导致摩擦阻力增大而卡滞，运动副间隙过大引起运动精度不足等。机械系统的功能表征量具有明显的退化特性，且各个功能指标之间具有相关性，需要综合退化过程和相关性对机械系统的功能失效进行分析。

图1 多功能机械系统竞争失效分析流程

任一失效模式的发生都会导致机械系统工作特性发生改变，进而其他失效模式不再有机会于同等条件下发生，即各个失效模式处于相互竞争的关系。对处于竞争失效模式下的机械系统进行可靠性分析时，首先需要对机械系统的运动特性进行分析，确定失效模式的具体表现形式(含突发型失效和退化型功能失效)；基于载荷随机特性和机械系统部件抵抗外界冲击的能力，确定所承载荷分布和阈值，进行突发型失效建模；根据机械系统各个功能属性定义功能表征量，基于运动几何学、运动学和动力学确定功能表征量及其影响因素之间的数学函数关系。对于难以直接表征的功能函数关系，可利用合适的代理模型构造每个功能表征量的功能函数，得到每个功能的累积分布函数；因功能之间的相关性，引入Copula函数构造多功能失效联合概率分布函数，进而得到考虑失效相关性的多失效模式机械系统的可靠度或失效概率。基于Copula函数的多功能机械系统竞争失效分析的完整流程如图1所示。

2 多失效模式机械系统失效建模方法

2.1 突发型失效建模方法

机械系统工作的过程中会承受相应的工作载荷，根据对多功能机械系统失效的描述，一些偶然因素会引起量值较大的冲击载荷，如重型车辆在恶劣路面行驶过程中的意外冲击，飞机飞行过程中结构颤振导致的冲击载荷等，从而造成部件的破坏失效。根据定义，这种冲击载荷在到来时间和每次载荷的量值大小上都具有不确定性。对于每一次突发冲击载荷而言，一旦超过机械系统构件所能承受的最大阈值，系统就会因构件破坏而失效，故机械系统不发生突发失效的概率表示为公式(1)。

式中：Wi为第i次突发载荷的大小；D0为系统可以承受最大载荷的阈值。突发载荷的到来时刻符合速率为λ的泊松过程，有如下公式成立

式中，N(t)和N(t+τ)分别表示截止时刻t和时刻t+τ突发载荷的发生次数。当t为0时，根据突发载荷的物理意义有N(t)=0，则

即

考虑到时刻t为止外界冲击载荷次数的所有可能情形，机构不发生突发失效的概率表示为公式(3)。

2.2 退化型功能失效建模方法

机构系统的特定功能可以使用具体的功能表征量进行定量化表征，根据指标判定机构功能是否正常。系统不发生功能失效的概率表示为公式(4)。

R(t)=P(|FS-F0|<δ)(4)

式中：F0为功能指标的理想值；FS为功能指标的实际取值；δ为指标理想值和实际值之间的允许误差。上述表达式的意义就是如果实际功能指标和功能指标的理想值之间的误差过大，则机械系统发生功能失效。

机械功能的实现受到多种因素的影响，包括各个部件特性、载荷、连接关系等。功能是否可以完成受到各个部件的共同作用影响。由于各部件在工作中伴随着不可避免的退化作用(例如磨损、老化、腐蚀等)，进而导致机构的功能也会发生退化。功能表征量与其影响因素存在函数关系，如公式(5)所示。

FS=f(C,L)(5)

式中：C表示机构部件特性和连接关系；L为载荷。函数表达式f的具体形式，需要根据具体机械系统的实际功能相应确定。机械系统的某些参数会随着工作循环次数和时间发生无法避免的退化(如磨损、老化)，进而通过功能函数，把这种退化传递给功能表征量，因此需要对参数的退化进行建模分析。漂移布朗运动是一种典型的用来对退化过程建模的方法。典型的漂移布朗运动的表达式为

X(t)=σB(t)+μt+X0(6)

式中：X(t)表示t时刻的退化量；σ为波动参数；B(t)为标准布朗运动；μ为漂移参数；X0表示退化量的初始值。通常情况下机械系统的退化是加速的，例如磨损的加剧会导致磨损接触面变得更加粗糙，增大接触面之间的摩擦力，从而进一步加速磨损；疲劳裂纹的出现使裂纹附近出现应力集中，进而加速裂纹的扩展。为了表示这种加速现象，在典型的漂移布朗运动中引入指数时间项，即把时间项t变为tq，其中q表示加速系数。此时漂移布朗运动的表达式变为公式(7)。

X(t)=σB(t)+μtq+X0(7)

当参数σ,μ,q和退化量的初始值X0为常量时，退化量X(t)～N(μtq+X0,σ2)。利用布朗运动对退化运动建模有几点优势[24]：①漂移参数的大小代表着退化过程中整体速率的高低，而引入指数的时间项则有效描述了工程中机械产品退化过程中的加速现象；②波动参数和标准布朗运动的乘积项代表着退化过程中退化量的随机性，即相同时间段内的退化量并非常值，而是一个有着波动向量的不确定值；③正态分布具有比较高效的计算效率和精确的计算值。

3 基于Copula函数的多功能机械系统可靠度计算方法

对于多功能机械系统而言，分析其可靠性的核心在于建立多功能表征量的联合分布函数(CDF)，但是复杂的相关关系，使得常规构造CDF的方法存在难度。而Copula函数则可以很好地描述相关性，并建立多变量的联合分布函数。

“Copula”一词起源于拉丁词汇“link”和“tie”，为一组边缘分布符合标准均匀分布变量的联合分布函数，表示为公式(8)。

C(u1,u2,…,un)=P(U1

(8)

式中，Ui为标准正态分布,即Ui～U(0,1),i=1,2,…,n。根据Skler定理[25],边缘分布为F1(x1),F2(x2),…,Fn(xn)的一维分布变量x1,x2,…,xn的n维联合分布函数F表示为公式(9)。

F(x1,x2,…,xn)=C(F1(x1),F2(x2),…,Fn(xn))

(9)

式中，C为对应的Copula函数。如果边缘分布F1(x1),F2(x2),…,Fn(xn)都是连续函数,则存在唯一Copula函数C满足公式(10)。

(10)

式中，u1=F1(x1),u2=F2(x2),…,un=Fn(xn)。

阿基米德Copula作为一类重要的函数族,因构造形式简单、计算结果准确而被经常使用。最典型的3种Copula函数分别为Clayton函数、Frank函数和Gumbel函数。3种函数形式如(以二元为例)表1所示。

表1 典型阿基米德Copula函数

当考虑竞争的因素时,如果系统需要保持可靠,则包括突发失效和功能失效在内的任何失效模式都不允许出现。综合公式(2)、公式(4)和公式(9),可得机械系统的可靠度表达式。

(11)

式中：FSi和F0i分别表示第i个功能失效模式指标的实际值和理想值；δi为第i个功能失效模式许用误差,其中i=1,2,…m；N(t)表示速率为λd的泊松过程。基于公式(11),即可利用Copula函数计算得到多失效模式相关机械系统可靠度及其演化规律。

4 案例:某飞机舱门锁机构

4.1 锁机构组成及工作原理分析

本节针对某型飞机舱门锁机构进行竞争失效分析。根据外场记录数据了解到该型号锁机构存在3种失效模式,失效模式一描述为因外界突发载荷造成的部件破坏,该失效模式属突发型失效;失效模式二描述为因摩擦因数恶化和弹簧刚度演化引起开锁过程阻力过大导致的功能失效,该失效模式属于退化型功能失效;失效模式三描述为关锁过程中总阻力过大而无法锁闭导致的功能失效,该失效模式属于退化型功能失效。

锁机构的组成结构如图2所示,主要包括作动筒、活塞、锁钩以及进行运动传递的3个连杆。作动

图2 锁机构组成

图3 锁机构工作状态

筒为驱动源,通过活塞带动各个连杆按照既定规律运动;锁钩为锁机构运动传递的目标构件;锁环固定在飞机机身,通过与锁钩的相对位置关系实现锁机构的关锁和开锁功能。锁机构的工作状态如图3所示。

开锁过程详细描述为,锁机构处于锁闭状态,在液压力的作用下锁钩克服锁环负载和锁机构的摩擦阻力,从锁钩锁定位置运动至锁钩最大旋转角度位置,这个阶段需要最大的液压驱动力克服运动过程中的阻力;之后从锁钩最大旋转位置运动到锁钩锁环脱离接触的位置,此时作用在锁钩上的载荷成为运动的动力而不是阻力;最后锁环锁钩脱离接触运动到打开位置,此时锁机构空载,液压力只需克服较小的重力。关锁过程详细描述为,作动筒驱动活塞和一系列连杆带动锁钩从初始位置开始运动到接触锁环,在锁钩锁环接触之前锁钩空载,整个机构的阻力主要来源于各个构件的重力;之后锁钩接触锁环继续运动到锁钩旋转到最大角度位置,此阶段锁钩开始接受来自锁环传递的负载,且负载的作用极大地增加了锁机构运动的阻力,需要提供较大的液压驱动力;最后从锁钩的最大旋转角度回复到锁死位置。这个过程锁环传递给锁钩的拉力变成了驱动力,液压力所需要克服的运动阻力值也较小。锁闭状态时,锁机构始终承受锁环位置传递的作用力,以及飞机运行过程中振动等原因引起的突发冲击载荷。

4.2 锁机构竞争失效建模

4.2.1 锁机构突发失效建模

锁机构处于锁闭状态时,锁钩与锁环紧密连接,受到机身传递来的突发冲击载荷作用,这类载荷的作用时间和量值都具有随机性。以Wi表示载荷的量值大小,H表示锁机构构件所能承受的临界值,根据公式(2),锁机构不发生突发失效的概率表示为

P(Wi

(12)

式中，Wi是独立同分布(independent identically distributed)的变量,表示每次出现突发载荷的量值。使用泊松过程表示突发载荷的随机出现时刻性,则机构不发生突发失效的概率表示为

(13)

式中，N(t)表示泊松过程,用于表示突发冲击载荷出现时刻的不确定性。不发生突发失效要求没有突发冲击载荷出现(即式中i取值为0)或出现的载荷量值皆小于阈值H。

4.2.2 锁机构退化型功能失效建模

根据对锁机构工作状态的分析,失效模式二和失效模式三发生于开/关锁过程中的驱动力无法克服运动阻力之时。开闭锁过程中的运动阻力都受到锁钩处摩擦因数、锁钩与锁环接触的角度、锁钩处的最大接触力和弹簧的弹性系数影响,且摩擦因数和弹簧刚度具有明显的随时间演化的特征,故此二类失效模式属于具有退化特征的功能失效,功能表征量分别为“开锁过程中的最大阻力”和“关锁过程中的最大阻力”,分别使用Fopen和Fclose表示。因此锁机构不发生开锁功能失效和关锁功能失效的概率分别表示为公式(14)和公式(15)。

Ropen(t)=Fopen(t)

(14)

Rclose(t)=Fclose(t)

(15)

式中：F(t)表示开锁或关锁过程中的最大阻力会因退化的存在随使用时间发生变化；FM表示作动筒所能提供的最大驱动力。为了提高计算效率,使用二阶响应面构造功能表征量及其影响因素之间的数学关系,步骤如下:

Step1 依据运动原理搭建参数化机构动力学仿真模型;

Step2 按照各输入参数的分布(见表1～3)进行拉丁超立方抽样,将抽样值代入仿真模型,得到成组的输入输出;

Step3 利用Step2中得到的结果构造响应面并进行精度验证。

构造得到的锁机构开关锁过程功能表征量(即所需驱动力)的响应面函数分别为

Fopen(t)=-4 416.10+1 382.69x1+

23.93x1x2+0.32x1x3-0.02x1x4-0.002 4x2x3-

0.000 39x2x4-1.06×10-6x3x4

Fclose(t)=4 557 030.76-33 176.20x1+

115.15x1x2+0.33x1x3+5.07x1x4-

0.009 1x2x3+0.093x2x4+0.000 23x3x4

式中，xi(i=1,2,3,4)为相互独立的影响因素,详细含义及取值见表2。当所需驱动力大于作动筒所能提供的最大驱动力时,就发生相应阶段的功能失效。根据公式(11)可得,该锁机构的可靠度表达式进一步推导为公式(16)。

(16)

表2 参数含义及取值

4.3 计算结果

综合公式(12)～(16)及锁机构开关锁过程功能表征量的响应面函数,将各参数取值(见表2)代入,即可得到计算结果。

图4为开锁过程所需驱动力的概率密度分布演化示意图,图5为闭锁过程所需驱动力的概率密度分布演化示意图。从图中可以看出,随着开闭锁循环次数的增加,开锁过程和闭锁过程所需驱动力的均值逐步上升,且分散性逐渐增大。

图4 开锁过程所需驱动力分布演化图 图5 闭锁过程所需驱动力分布演化图图6 可靠度变化曲线

根据公式(16)可以对锁机构可靠度进行计算,锁机构可靠度的演化结果如图6所示。图中横坐标为使用锁机构开闭锁循环次数表示的等效时间,纵坐标为可靠度。分别将基于3种类型Copula函数、独立假设情况与蒙特卡洛方法计算所得结果相比,各个方法的均方根误差(root mean square error,RMSE)如表3所示。根据计算结果可以得到,锁机构的可靠度随着使用而逐渐发生退化,以蒙特卡洛方法计算结果为对照,相比于独立假设,本文所提基于Copula函数的计算结果误差更小,且其中Gumbel函数形式误差最小,更适合于本案例的研究。

表3 各类方法所得结果误差

5 结 论

本文针对处于竞争失效模式作用下的多功能机械系统，基于Copula函数建立了可靠性分析方法。Copula函数可以有效地建立多功能相关情况下功能失效的联合概率分布函数，进而结合突发失效，计算得到机械系统的总可靠度，并利用所提方法，对某飞机舱门锁机构进行了可靠性分析。根据分析和计算的结果，可以得到如下结论：

1) 飞机舱门锁机构开闭锁所需驱动力随时间发生演化，且分散性逐渐变大；

2) 本文所提基于Copula函数的考虑失效模式相关性的可靠性分析结果更接近蒙特卡洛方法所得结果，证明本文所提方法的正确性；

3) 针对本文所提案例，Gumbel形式的Copula函数计算结果的RMSE为6.19×10-4，误差最小，更适合于对此案例进行分析。