近年来，许多****对运动机构的可靠性及灵敏度进行了大量的研究。孙中超等[1]研究了销轴间隙的不确定性对连杆机构运动精度的全局灵敏度，结果表明，销轴间隙均值的变化对重要度影响较小，标准差的变化对重要度影响较大。张娟和王艳艳[2]以某型飞机起落架收放机构为研究对象，研究了各输入变量对起落架收放机构时变可靠性的影响大小和作用机理。国志刚等[3]给出了单个铰链与多个铰链的磨损可靠性分析方法和计算模型，并采用曲柄滑块机构验证了模型与方法的有效性和可行性。成凯等[4]提出了随机输入变量和模糊输入变量的全局灵敏度新指标，并采用Kriging代理模型提高了混合不确定性全局灵敏度指标的计算效率，通过算例验证了方法的高效性和准确性。唐成虎等[5]研究了飞机缝翼机构发生性能退化后的参数灵敏度，并采用神经网络代理模型，有效提高了机构灵敏度分析效率。游令非等[6]将混合不确定性引入机构时变问题中，建立了机构产品的模糊-随机时变可靠性模型，采用四连杆机构验证了方法的可行性，具有较高的计算精度。
本文在ADAMS软件中建立了飞机舱门泄压阀机构多体动力学仿真模型，研究在销轴处发生磨损后，对泄压阀运动定位精度的影响。建立了销轴磨损量和泄压阀定位角度差值的Kriging代理模型，并主动地挑选符合学习准则的点加入到样本池中以提高代理模型精度。基于构建好的代理模型，采用Monte Carlo方法计算失效概率和全局灵敏度指标，计算结果可以为工程实际提供指导作用。
1 泄压阀机构多体动力学仿真分析 泄压阀机构作为飞机舱门系统的主要组成部分，在舱门打开阶段对于平衡飞机内外气压具有重要作用。本文使用ADAMS软件建立了泄压阀机构的多体动力学仿真模型。按照泄压阀的运动特征添加约束、驱动、接触力，并建立必要的零部件模型。在泄压阀机构发生磨损的部位施加平面副进行约束，保证其运动不偏离回转平面。所建立的仿真模型能够准确考虑各个销轴由于长期磨损造成的销轴间隙增大，导致运动精度降低的情况。
图 1为建立的泄压阀机构仿真模型及考虑发生磨损的5个销轴ri(i=1, 2, 3, 4, 5)所在位置。在内手柄处提供驱动，通过各个构件传递运动关系，最终使泄压阀打开，完成既定功能。图 2中，3条曲线分别表示在未发生磨损、磨损量为0.1 mm、磨损量为0.2 mm时泄压阀打开角度随时间的变化。可知，发生磨损后会导致泄压阀运动精度降低，进而造成舱门部分功能丧失。
 |
| 图 1 飞机舱门泄压阀机构示意图 Fig. 1 Schematic diagram of aircraft door's relief valve mechanism |
| 图选项 |
 |
| 图 2 不同磨损量下泄压阀的打开角度 Fig. 2 Opening angles of relief valve under different wear amounts |
| 图选项 |
2 泄压阀机构可靠性与灵敏度分析 2.1 基于Archard理论的销轴磨损量计算 英国****Archard[7]提出的磨损模型广泛应用在工程当中。因此，本文使用Archard黏着磨损模型公式对飞机舱门泄压阀机构各个销轴磨损量进行分析计算。Archard磨损计算表达式为
 | (1) |
式中：V为磨损体积；K为磨损因数；P为接触面的法向压力；L为相对滑移距离；H为材料布氏硬度，通常取H≈3σs，σs为材料屈服强度。
在单次磨损过程中，销轴相对滑移距离可由式(2)得出：
 | (2) |
单个销轴的磨损体积为
 | (3) |
式中：b为销轴高度；r为销轴半径；φ为销轴相对旋转角度；n为销轴磨损次数。
联合式(1)~式(3)，求得销轴半径的磨损量Δr为
 | (4) |
发生磨损的销轴材料分别为Q275钢和Q215钢，查阅资料[8]可知，其σs分别为275 MPa和215 MPa。在计算销轴半径的磨损量时，还需要确定的参数主要是接触面的法向压力P。本文采用如下方法获得较为准确的法向压力。
当未发生磨损时，通过所建立的多体动力学仿真模型仿真测得各个销轴接触处的接触力P0i(i=1, 2, 3，4, 5)，随后以P0i计算飞机舱门泄压阀机构磨损n次过后的磨损量，将该磨损量反映在模型中得到新的正压力Pni。当n取值越小时，越接近实际磨损时的法向压力值。
图 3为随着磨损次数的增加，5个销轴处的法向压力变化情况。
 |
| 图 3 不同销轴处的法向压力随磨损次数变化的曲线 Fig. 3 Change of normal pressure with wear times for different pins |
| 图选项 |
图 4为销轴发生磨损后磨损量随磨损次数的变化曲线。可以看出，在不同磨损次数下，不同销轴的磨损情况有较大差异。本文将5个销轴ri(i=1, 2, 3, 4, 5)处的磨损量Δri(i=1, 2, 3, 4, 5)视为随机变量，记为Xi(i=1, 2, 3, 4, 5)，假设磨损量服从正态分布。随机变量的均值随着磨损次数而不断变化，变异系数取0.1。表 1为30 000次磨损时输入变量的概率分布信息。
 |
| 图 4 不同销轴处的磨损量随磨损次数变化的曲线 Fig. 4 Change of wear amounts with wear times for different pins |
| 图选项 |
表 1 30 000次磨损时输入变量概率分布信息 Table 1 Probability distribution information of input variables at 30 000 times of wear
| 变量 | 分布类型 | 均值/mm | 变异系数 |
| X1 | 正态分布 | 0.255 | 0.1 |
| X2 | 正态分布 | 0.065 | 0.1 |
| X3 | 正态分布 | 0.074 | 0.1 |
| X4 | 正态分布 | 0.101 | 0.1 |
| X5 | 正态分布 | 0.117 | 0.1 |
表选项
2.2 泄压阀机构的可靠性模型 如2.1节所述，销轴磨损量X=[X1, X2, X3, X4, X5]作为输入变量，记泄压阀允许的角度误差阈值为Δ0，发生磨损后的泄压阀定位角度为Y，未发生磨损时的定位角度为Y*(Y*=16.87°)，二者之差为Δ=|Y-Y*|。首先定义销轴磨损量X与泄压阀定位角度Y的响应函数关系为
 | (5) |
建立反映泄压阀机构定位精度的功能函数为
 | (6) |
当角度差Δ超过阈值Δ0时即为失效，失效域定义为
 | (7) |
失效概率为
 | (8) |
式中：
 | (9) |
其中：IF(X)为失效域的指示函数；D为变量空间；E[·]为数学期望算子。
根据可靠性理论中的Monte Carlo方法[9]，可以采用失效域指示函数的数学期望来计算泄压阀定位角度的失效概率，即
 | (10) |
式中：N为总抽样次数；Nf为落入失效域的样本个数。
2.3 全局灵敏度指标 可靠性分析可以计算出机构的失效概率或可靠度，灵敏度分析则研究输入变量的分布参数或不确定性对机构输出性能统计特征的影响[10]。灵敏度分析包括局部灵敏度分析和全局灵敏度分析，本文采用全局灵敏度指标来衡量销轴磨损量Xi(i=1, 2, 3, 4, 5)作为随机变量时对泄压阀机构定位精度Y及功能函数值g(X)的影响。
2.3.1 基于方差的全局灵敏度 基于方差的重要性分析目的是得到各输入变量对响应量方差的贡献，采用Sobol[11]提出的模型分解方法，可以将响应函数Y=h(X)分解如下:
 | (11) |
式中：h0为响应函数h(X)的均值。
通过将响应量方差V(Y)分解为各项方差之后，可得
 | (12) |
根据式(12)，响应量方差可以分解为每个输入变量的方差贡献。其中，Vi是输入变量Xi对响应量一阶方差的贡献，表达式如下：
 | (13) |
由此可知，Vi能够反映出在固定Xi后，响应量方差的平均变化情况，Vi能够较好地反映输入变量Xi对响应量方差的影响。通常只考虑一阶方差的影响时，式(12)可以简化为
 | (14) |
一阶方差贡献Vi一般被称为主贡献，由此定义基于方差的全局灵敏度指标δi：
 | (15) |
δi能够衡量输入变量Xi对于模型响应量变异性的影响。在泄压阀机构中，即反映第i个销轴磨损量对泄压阀定位角度Y的方差的影响程度。
2.3.2 基于概率分布的全局灵敏度 除方差外，概率分布函数(包括概率密度函数和累积分布函数)也是描述随机性的重要统计特征量[12-13]。对于响应函数Y=h(X1, X2, …, Xn)记其无条件概率密度函数和分布函数分别为fY(Y)和FY(Y)，当输入变量Xi取其实现值Xi*时，可以得到Y的条件概率密度函数和分布函数fY|Xi(Y)和FY|Xi(Y)。
Borgonovo[14]提出了全局灵敏度用来衡量输入变量对响应量概率密度函数(PDF)的影响。当响应量Y从-∞到+∞变化时，输入变量Xi取其实现值Xi*时对响应量分布密度的影响可以用式(16)积分值s(Xi)来衡量：
 | (16) |
当Xi按照其分布规律取所有可能实现值时，Xi对响应量概率密度累积影响的平均值可由式(16)的期望值，即EXi[s(Xi)]来描述，由此定义输入变量Xi基于响应量概率密度分布的矩独立全局灵敏度指标ηi如下：
 | (17) |
ηi从概率密度函数出发，衡量了泄压阀机构中第i个销轴磨损量对于泄压阀定位角度Y概率特征的影响。
2.3.3 基于失效概率的全局灵敏度 对于可靠性问题而言，研究输入变量对于失效概率或可靠性的影响程度大小更接近可靠性设计所需的信息。基于失效概率的全局灵敏度指标主要用来衡量输入变量在其整个分布域内变化时对失效概率的平均影响[15]。记基于失效概率的全局灵敏度指标为
 | (18) |
通过求解εi，可以确定第i个销轴磨损量对泄压阀可靠性模型失效概率的影响程度。输入变量的εi灵敏度指标越大，说明该变量对可靠性的影响就越大。
2.4 基于Kriging代理模型的可靠性与灵敏度分析 可靠性和灵敏度计算依赖大量的样本数据，尤其对于计算成本高的模型，采用代理模型能够有效地提高计算效率。目前，工程领域常用的代理模型有响应面、神经网络和Kriging代理模型等。Kriging代理模型作为一种估计方差最小的无偏估计模型，具有全局近似与局部随机误差相结合的特点，能够提供很好的全局预测[16]。
Kriging代理模型可以近似表达为一个随机分布函数和一个多项式之和，如下：
 | (19) |
式中：gK(x)为未知的Kriging代理模型；f(x)=[f1(x), f2(x), …, fp(x)]T为随机向量x的基函数，提供了设计空间内的全局近似模型；β=[β1, β2, …, βp]为回归函数待定系数，其值可通过已知的响应值估计得到，p为基函数的个数；z(x)为一随机过程，是在全局模拟的基础上创建的期望为0且方差为σ2的局部偏差，其协方差方程为
 | (20) |
其中：R(xi, xj)为样本池中任意2个样本点xi、xj的相关函数，Kriging代理模型中常用高斯相关函数，其形式为
 | (21) |
式中：θk为未知的相关参数，一般采用最大似然估计得到：
 | (22) |
回归系数
 | (23) |
 | (24) |
式中：1为元素为1的k×1向量；y为训练样本响应值。
对于任意的未知点x，通过Kriging代理模型预测该点函数值gK(x)服从高斯分布，即gK(x)~N(μgK(x), σ2gK(x))，对于均值和方差计算可以采用MATLAB工具箱DACE[17]来实现。
Kriging代理模型对于预测点的估计为无偏估计，其预测误差σ2gK(x)也通常被称为均方误差(MSE)。由于其能够反映预测值的稳健性，可将其作为学习函数来更新泄压阀Kriging代理模型，提高代理精度，即
 | (25) |
所用加点收敛准则定义如下：
 | (26) |
式中：X(i)表示根据式(25)第i次所挑选的样本点。当Cr < 5×10-5时，满足收敛准则，即Kriging代理模型停止加点。该学习函数能够挑选出预测误差大的点加入到样本池中，所用收敛准则能减小在加点过程中2次预测产生的差异，从而进一步提高代理模型的全局精度，为计算全局灵敏度构建精度更高的代理模型。
针对飞机舱门泄压阀机构发生磨损时的可靠性和灵敏度分析，本文采用主动学习Kriging代理模型构建泄压阀输入变量与输出响应量的隐式关系。基于构建好的功能函数，采用Monte Carlo方法计算泄压阀机构的失效概率，并计算不同的全局灵敏度分析指标，得到泄压阀输入变量对输出响应量的重要性排序，分析流程如图 5所示，具体步骤如下：
 |
| 图 5 泄压阀机构可靠性和灵敏度分析流程 Fig. 5 Flowchart of reliability and sensitivity analysis for relief valve mechanism |
| 图选项 |
步骤1 ??创建总样本池S。采用Sobol序列抽样方法得到分布在整个样本空间的总样本池S。
步骤2 ??构建初始训练样本池S0。从总样本池S中随机选择N0=500组作为初始训练样本，调用ADAMS得到500组响应值组成初始训练样本池S0。
步骤3 ??训练Kriging代理模型。从总样本池S中根据MSE学习函数选择新的样本点更新训练样本池S0，不断拟合代理模型，直至满足收敛条件。
步骤4?? 失效概率和灵敏度计算。基于构建好的Kriging代理模型，采用Monte Carlo方法计算不同磨损次数下的失效概率和3种全局灵敏度指标。
3 结果分析 对于复杂的工程模型，可靠性和灵敏度计算需要大量的样本值，构建代理模型是提高计算效率的有效途径。图 6为不同磨损次数下(1×104, …, 6×104次)所构建代理模型的加点收敛过程。由于未加点时假设σgK(X)=1，所以在加点初期，衡量代理精度的指标值Cr变化较大。当加点数量超过10时，建模精度开始进入稳定期。当加点个数超过40时，逐渐满足收敛条件，完成代理模型的构建。
 |
| 图 6 代理模型加点收敛过程 Fig. 6 Convergence progress of surrogate model adding points |
| 图选项 |
图 7给出了不同磨损次数下，随着角度阈值的改变，泄压阀机构失效概率的变化规律。随着角度阈值Δ0的增大，泄压阀机构的失效概率逐渐减小，表明机构规定角度阈值越大，即对精度要求越低，机构可靠性越高。
 |
| 图 7 不同磨损次数下失效概率随角度阈值的变化规律 Fig. 7 Rule of failure probability varying with threshold under different wear times |
| 图选项 |
通过采用不同的全局灵敏度分析方法，其结果能够反映出泄压阀发生磨损退化后，输入变量分布参数的改变对输出响应量的影响。泄压阀机构从1×104次到6×104次磨损过程中，各个销轴发生磨损形式的性能退化后，求得全局灵敏度指标分别如图 8和图 9所示。
 |
| 图 8 δi随销轴磨损次数的变化规律 Fig. 8 Rule of δi varying with pins wear times |
| 图选项 |
 |
| 图 9 ηi随销轴磨损次数的变化规律 Fig. 9 Rule of ηi varying with pins wear times |
| 图选项 |
基于方差的全局灵敏度指标衡量输入量对于响应量方差的影响，图 8即为不同销轴δi随磨损次数的变化。从1×104次到6×104次磨损过程中，X1基于方差的灵敏度指标一直保持最高，X2的灵敏度指标最低。在泄压阀磨损初期(1×104次至2×104次)，X1灵敏度逐渐升高，由2×104次到3×104次过程中，其余销轴的灵敏度缓慢上升，X1灵敏度逐渐下降，在磨损4×104次时，各个销轴的灵敏度变化达到极值，在随后的4×104次到6×104次磨损过程中，基于方差的全局灵敏度指标排序逐渐保持稳定。如图 9所示，ηi变化趋势与δi大致相同，ηi从响应量概率密度函数角度出发，衡量了输入变量的随机取值对响应量分布的影响。在不同磨损次数下，ηi中X3与X4的排序与δi略有不同，这也反映了仅从单一角度衡量输入变量对响应量的影响程度是不全面的。目前的研究结果表明，尚没有任何一种指标能够全面地反映输入不确定性对响应量不确定性的影响，应当根据分析问题的侧重点来选择合适的指标。
不同磨损次数下对应的基于失效概率的全局灵敏度指标εi的求解依赖于泄压阀机构功能函数落入失效域，即失效域指示函数为1时才能求出正确解。由图 7可知，若给定阈值为1时，失效概率的变化范围仅在3×104次附近波动。所以在给定阈值时，仅计算泄压阀在3×104次至3.5×104次时εi的变化情况，计算结果如图 10所示。基于失效概率的全局灵敏度指标εi的变化规律，与前2种指标排序类似，磨损程度最大的X1指标值一直保持最高，X2指标值最低。X4与X5随着磨损次数的增加指标值逐渐相同，X3指标值先减小后增大。图 11为3×104次磨损时，基于失效概率的全局灵敏度指标εi随所设定的角度阈值的变化规律。当阈值从0.8变化至1.1时，5个销轴磨损量基于失效概率的全局灵敏度指标εi经历了先增大后减小的过程，但是在排序上并没有实质性的变化。在不同的阈值下，X1的指标值都保持最高，X5次之，X2最低，X3与X4大致相同。
 |
| 图 10 给定阈值下εi随销轴磨损次数的变化规律 Fig. 10 Rule of εi varying with pins wear times under fixed threshold |
| 图选项 |
 |
| 图 11 3×104次磨损时εi随角度阈值的变化规律 Fig. 11 Rule of εi varying with angle threshold under 3×104 wear times |
| 图选项 |
4 结论 1) 在ADAMS中建立了某型飞机舱门泄压阀机构的多体动力学仿真模型，基于Archard磨损理论计算得到了5个销轴处磨损量随磨损次数的变化规律。
2) 考虑磨损量的不确定性，建立了泄压阀机构的可靠性模型，结果表明，在不同磨损次数下，阈值大小对于失效概率的影响有较大差异。计算得到3种全局灵敏度指标，从不同角度分析输入变量对响应量的影响程度。对于磨损量较大的销轴，其全局灵敏度指标值也较大，不同指标随磨损次数的变化规律并不相同。
3) 为了提高求解效率，本文采用了Kriging代理模型，只需对原始模型进行少量样本点选取即可精确代理，大大减少了计算量。
参考文献
| [1] | 孙中超, 喻天翔, 崔为民, 等. 销轴间隙对连杆机构运动精度的重要性测度分析[J]. 中国机械工程, 2014, 25(21): 2874-2879. SUN Z C, YU T X, CUI W M, et al. Importance measure of revolute joint clearance about motion accuracy of linkage mechanism[J]. China Mechanical Engineering, 2014, 25(21): 2874-2879. DOI:10.3969/j.issn.1004-132X.2014.21.007 (in Chinese) |
| [2] | 张娟, 王艳艳. 起落架收放机构时变可靠度与时变全局可靠性灵敏度分析[J]. 机械科学与技术, 2018, 37(3): 358-363. ZHANG J, WANG Y Y. Research on time-variant reliability and time-varying global reliability sensitivity of aircraft landing gear mechanism[J]. Mechanical Science and Technology for Aerospace Engineering, 2018, 37(3): 358-363. (in Chinese) |
| [3] | 国志刚, 冯蕴雯, 冯元生. 铰链磨损可靠性分析及计算方法[J]. 西北工业大学学报, 2006, 24(5): 644-648. GUO Z G, FENG Y W, FENG Y S. Wear reliability of linkage under periodic loading[J]. Journal of Northwestern Polytechnical University, 2006, 24(5): 644-648. DOI:10.3969/j.issn.1000-2758.2006.05.024 (in Chinese) |
| [4] | 成凯, 吕震宙, 石岩. 基于可能性矩的混合不确定性全局灵敏度分析[J]. 北京航空航天大学学报, 2017, 43(8): 1705-1712. CHENG K, LU Z Z, SHI Y. Global sensitivity analysis under mixed uncertainty based on possibilistic moments[J]. Journal of Beijing University of Aeronautics and Astronautics, 2017, 43(8): 1705-1712. (in Chinese) |
| [5] | 唐成虎, 周长聪, 侯伟, 等. 考虑性能退化的飞机典型部件灵敏度分析[J]. 西安交通大学学报, 2019, 53(4): 158-166. TANG C H, ZHOU C C, HOU W, et al. A sensitivity analysis of typical aircraft components with performance degradation[J]. Journal of Xi'an Jiaotong University, 2019, 53(4): 158-166. (in Chinese) |
| [6] | 游令非, 张建国, 翟浩, 等. 模糊-随机混合参数的机构运动可靠度计算方法[J]. 北京航空航天大学学报, 2019, 45(4): 714-721. YOU L F, ZHANG J G, ZHAI H, et al. Computation method on motional reliability of mechanism under mixed parameters with fuzziness and randomness[J]. Journal of Beijing University of Aeronautics and Astronautics, 2019, 45(4): 714-721. (in Chinese) |
| [7] | ARCHARD J F. Contact and rubbing of flat surfaces[J]. Journal of Applied Physics, 1953, 24(8): 981-988. DOI:10.1063/1.1721448 |
| [8] | 沈宁福. 新编金属材料手册[M]. 北京: 科学出版社, 2003. SHEN N F. New handbook of metal materials[M]. Beijing: Science Press, 2003. (in Chinese) |
| [9] | 吕震宙, 宋述芳, 李洪双, 等. 结构机构可靠性及可靠性灵敏度分析[M]. 北京: 科学出版社, 2009. LU Z Z, SONG S F, LI H S, et al. Reliability and reliability sensitivity analysis of structural mechanism[M]. Beijing: Science Press, 2009. (in Chinese) |
| [10] | 吕震宙, 李璐祎, 宋述芳, 等. 不确定性结构系统的重要性分析理论与求解方法[M]. 北京: 科学出版社, 2015. LU Z Z, LI L Y, SONG S F, et al. Importance analysis theory and solution method of uncertain structural system[M]. Beijing: Science Press, 2015. (in Chinese) |
| [11] | SOBOL I M. Global sensitivity indices for nonlinear mathematical models and their Monte Carlo estimates[J]. Mathematics and Computers in Simulation, 2001, 55(1): 271-280. |
| [12] | ZHOU C C, TANG C H, LIU F C, et al. Regional moment-independent sensitivity analysis with its application in engineering[J]. Chinese Journal of Aeronautics, 2017, 30(3): 1031-1042. DOI:10.1016/j.cja.2017.04.006 |
| [13] | ZHOU C C, ZHANG Z, LIU F C. Sensitivity analysis for probabilistic anti-resonance design for aeronautical hydraulic pipelines[J]. Chinese Journal of Aeronautics, 2019, 32(4): 948-953. DOI:10.1016/j.cja.2019.01.001 |
| [14] | BORGONOVO E. A new uncertainty importance measure[J]. Reliability Engineering and System Safety, 2007, 92(6): 771-784. DOI:10.1016/j.ress.2006.04.015 |
| [15] | LI L Y, LU Z Z, FENG J, et al. Moment-independent importance measure of basic variable and its state dependent parameter solution[J]. Structural Safety, 2012, 38(5): 40-47. |
| [16] | 韩忠华. Kriging模型及代理优化算法研究进展[J]. 航空学报, 2016, 37(11): 3197-3225. HAN Z H. Kriging surrogate model and its application to design optimization: A review of recent progress[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(11): 3197-3225. (in Chinese) |
| [17] | LAURENT L, RICHR R L, SOULIER B, et al. An overview of gradient-enhanced metamodels with applications[J]. Archives of Computional Methods in Engineering, 2019, 26(1): 61-106. DOI:10.1007/s11831-017-9226-3 |
