小样本下考虑磨损退化的机构时变可靠性分析

王佳聪， 高 鹏， 常筠袖， 田思祺， 崔永海

（1.辽宁石油化工大学 机械工程学院，辽宁 抚顺 113001； 2.中国石油锦州石化公司，辽宁 锦州 121001；3.抚顺新钢铁有限责任公司，辽宁 抚顺 113001）

机构作为机械产品的重要组成部分，主要用来传递运动或改变运动形式。由于制造、装配误差等非累积性因素会引起机构在运转中产生累积性误差，进而导致机构运动精度逐渐下降，因此对在产品设计后与投入使用前进行可靠性分析十分必要。在分析各种参数对机械机构运转影响时，考虑经济成本、时间成本等原因，学者们往往采用动力学仿真软件等来模拟机构的运转情况。I.KHEMILI 等[1]利用ADAMS 软件对曲柄滑块刚柔耦合机构进行了仿真，考虑间隙故障，将模拟结果与从自己的试验台获得的实验数据进行比较，证明了仿真结果的准确性。Z.F.BAI 等[2]基于Archard 磨损模型，通过使用ADAMS 软件对四杆机构完成了动力学仿真，给出了带回转间隙连接的多体系统的动力学响应。A.GUMMER 等[3]提出了一种利用多体仿真工具RecurDyn 对带有不完美转动关节的机构进行仿真的方法，开发了一个计算给定恢复系数下模型参数的专用程序。樊智敏等[4]应用ADAMS 软件建立了双渐开线齿轮结构，综合考虑齿面摩擦与齿轮啮合刚度二因素,对比分析了齿轴分解参数对齿轮性能的影响。在进行可靠度计算时，经常需要大量数据样本提供支撑。因此，学者们又利用一些可以减小样本量但能保持一定精度的计算方法近似地得到了机构可靠度。C.G.BUCHER 等[5]最早提出利用平均值和标准差的基本统计信息来提高效率和精度，使用响应面法（RSM）评估可靠性，以此提高计算效率。S.K.YOO 等[6]采用传统响应面法（RSM）和连续标准响应面法（SRSM）进行了基于双向渐进结构优化（BESO）的可靠性拓扑优化（RBTO）。张春宜等[7-8]提出了极值响应面法，通过以动态参数的随机样本值作为输入、机构应力和变形的最大值为输出响应，分析了柔性机构可靠性。D.Q.ZHANG等[9]在基于时间的可靠性分析中采用响应面法，计算了随时间变化的故障概率和可靠性指标，并用于估计随机载荷作用下非确定性结构的时变可靠性。A.T.ESHGHI 等[10]提出自适应改进响应面法（AIRSM），并对其与降维法（DR）的可靠性进行了比较。

本文为了在节省分析时间、减少样本数量的前提下保证计算精准度，考虑机构旋转副磨损退化故障这一累积性因素，从仿真到可靠度计算给出了一完整可靠性分析模型。在仿真过程中，此模型使用ADAMS 与MATLAB 软件的联合仿真得出一定数量的样本；提出了ARIMA 与RSM 相结合的方法（ARI-RSM 方法），并通过此方法实现样本预测；通过应力-强度干涉模型建立极限状态方程，计算时变可靠度。最后，以曲柄滑块为例，验证了方法的可行性。

1 仿真模型

1.1 磨损量的计算

在机构运行过程中，含有间隙的铰链会导致销轴与轴套间发生相对运动或接触，此时会产生磨损，逐渐引起铰链接触面的形状和尺寸的改变。随着销轴与轴套磨损的加剧，铰链间隙逐渐变大，并呈现出非常规形状特征，影响铰链连接的多体动力学特性，降低机构的运行精度或导致机构发生故障，进而影响机构的使用寿命。因此，对磨损量进行精准计算及预测非常有必要[11]。磨损量的计算学者们多数采用三种方法，即经验公式、克拉盖尔斯公式及Archard 磨损计算公式（Archard 公式，下同）[12-13]。其中，经验公式在使用前需要通过大量实验计算出材料的磨损率与一些待定参数；克拉盖尔斯公式虽有一定的通用性，但计算繁琐；Archard 公式简单且参考资料较多，计算精度较高。因此，本文采用Archard 公式计算磨损量[13]。Archard 公式的表达式为：

式中：V为体积磨损量，mm³；S为相对滑动距离，mm；Fn为碰撞产生的法向作用力，N；K为磨损系数；H为材料硬度，HB。S和Fn均由仿真过程得出。对于所需的磨损深度，将所得体积磨损量与孔轴接触面积相除即可得到。

ADAMS 软件采用IMPACT 函数表达构件间的碰撞，应用非线性弹簧阻尼器模型计算碰撞力[11]。该模型在接触深度为负数时处于失效状态，在接触深度为正数时处于激活状态，可表示为：

式中：k为刚度系数，N/mm；d为接触深度，mm；STEP 为阶跃函数；dmax为最大接触深度，mm；Cmax为最大阻尼系数，m/s；δ为碰撞深度，mm；n为指数。

刚度系数k的表达式见式（3）。

式中：R为构件的曲率半径，mm；v为构件的泊松比；E为构件的弹性模量，Pa；σ为常数，可由泊松比和弹性模量计算得到；下标i、j 分别表示不同杆件。

在仿真过程中，通过测量可得两构件的转动角位移，再通过式（4）即可得相对滑动距离。

式中：ϕ为两构件转动角位移之差，（°）。

1.2 联合仿真

ADAMS 是一个多体动力学仿真软件，其求解器采用多刚体系统动力学理论中的拉格朗日方程，对机械系统的运动学和动力学都具有很好的分析能力。但是，此软件由于存在局限性，并不能在仿真过程中模拟构件磨损退化现象，需要每次仿真后导出实验数据，计算单次磨损量，依据计算结果更改仿真模型中构件尺寸后再次进行仿真，对一定数量的样本机构进行可靠性分析时，需要大量的时间及人力，并不能满足实际需求。

MATLAB 软件在非线性动态系统的建模和仿真中具有强大的计算能力，被广泛应用于各种科学研究与工程设计中。虽然可以通过动力学方程建立机构模型并应用此软件进行计算，但消耗大量的时间。MATLAB 中的Simulink 模块能够方便搭建控制系统，并与ADAMS 之间存在交互接口，因此本文模型的仿真过程由ADAMS 与MATLAB 共同完成。

本文模型以ADAMS 为基础。 首先，在ADAMS 中建立机构模型，根据计算磨损量公式中所用到的数据，建立测量模块，将测量结果添加到输出单元中；将考虑随机变量的机构尺寸添加为设计变量，软件的Adams Controls 模块与MATLAB软件的Simulink 模块相互连接；在MATLAB 中编写程序，设定循环次数，生成随机变量样本，通过Simulink 模块调用ADAMS，动力学分析过程由ADAMS 软件完成；通过得到的测量结果，在MATLAB 中计算磨损量。

2 可靠度计算模型

2.1 ARIMA 与RSM

由于磨损退化，机构的轴孔半径随时间逐渐缩小，因此可以作为时间序列问题进行分析。本文采用差分自回归移动平均模型（ARIMA）对轴孔半径的变化进行预测[14]。ARIMA 结合了自回归模型（AR）、移动平均模型（MA）和差分法。通过差分法可将非平稳时间序列转化为平稳时间序列，其具体模型为：

式中：△为差分符号；b为差分次数；N为误差估计；μ为常数；φ为自回归参数；ε为误差；p为自回归模型阶数；q为移动平均模型阶数；θ为移动平均相关参数；j为与回溯白噪声数量相关的参数；r为与时间步长相关的参数。p和q的值通过自相关函数和偏自相关函数确定。

响应面法（RSM）是利用多项式函数来近似逼近隐式功能函数的一种方法，因其具有计算准确性与便捷性的特点，在可靠性分析过程中得到了广泛的应用[15]。为保证计算量与计算精度，本文的机构功能函数由不含交叉项的二次多项式来代替，其表达式为：

式中：ġ(x)为响应近似值；xi为随机变量；a、bi、ci均为待定系数；m为随机变量个数。

2.2 计算模型

本文建立了ARIMA 与RSM 相结合的方法，即ARI-RSM 方法，得出了机构预测响应值。由式（5）及式（6）可得：

机构时变误差函数：

式中：t为机构运转时间；y为样本结束时间；T为预测结束时间；X为样本随机变量；△d Ni为预测随机变量；G(X,t)为机构理想响应值；g(X,t)为样本的实际响应值为预测响应值。根据应力强度干涉模型[16]，可得机构时变极限状态函数：

机构时变可靠度（R(t)）的表达式为：

此可靠度计算模型的基本思路为：

（Ⅰ）将MALAB 中通过仿真计算得到的孔轴半径样本作为RSM 的输入、对应的机构响应值作为RSM 的输出，计算待定系数，得出二次多项式；

（Ⅱ）通过ARIMA 对存在磨损退化故障的孔轴半径样本进行分析，得出后续磨损过程中的预测结果；

（Ⅲ）将步骤（Ⅱ）中的预测结果输入至由步骤（Ⅰ）得到的二次多项式，得出估计机构响应值，完成样本扩增；

（Ⅳ）根据应力强度干涉模型，构建极限状态方程，计算机构时变可靠度。

3 可靠性分析模型检验

可靠性分析模型检验分为两部分，一部分为仿真检验，另一部分为可靠度计算检验。

曲柄滑块机构作为典型的连杆机构，其应用范围广泛，且在运行过程中具有销轴磨损量大的特点，故本文选取的实例模型为曲柄滑块机构。曲柄滑块机构如图1 所示，其参数见表1。此次仿真实验做简化处理，只考虑销轴A 和销轴B 的磨损退化情况。

表1 曲柄滑块机构参数

图1 曲柄滑块机构

3.1 仿真结果验证与分析

设定曲柄旋转一周为1 个磨损周期，以周期为单位对机构参数重构，使用ADAMS 及ADAMSMATLAB（ABADAMS 和MATLAB 联合仿真）对比了曲柄滑块机构在10 个周期内销轴A、B 半径的变化，以及第1 周期滑块质心位移的变化和整体仿真时间的变化，结果如图2—4 所示。

图2 销轴A 及销轴B 半径的变化曲线

由图2、图3 可以看出，在2 种仿真方式下，2 个销轴的半径变化及滑块位移变化完全吻合，说明在联合仿真过程中并没有造成数据损失及计算误差。

图3 第1 周期滑块质心位移变化曲线

由图4 可以看出，这一仿真方法具有很好的便捷性，除了能保证精度，还能大大缩短仿真时长，可用于可靠度计算的前期样本准备。

图4 整体仿真时间变化曲线

3.2 可靠度计算验证与分析

考虑到在这一机构模型中曲柄旋转1 周产生的磨损量小于5×10－6mm，实验设计曲柄旋转2 000 周为1 个磨损周期，共进行10 000 次的蒙特卡洛仿真实验来模拟机构运转50 周期。其中，样本中的一部分作为本文模型的样本，用于预测扩展样本数量，其余作为模型预测结果的检验。以销轴A、B的半径作为式（8）的随机变量[17]，滑块质心最远端位移作为式（8）的响应值，位移误差阈值设定为0.18 mm。

以蒙特卡洛方法的计算结果作为准确值，采用5个样本与预测数量的关系（指的是模型样本与扩展样本数量之比），将本文模型得到的计算结果与蒙特卡洛法进行了对比验证[18]，结果见图5。图5 中，括号内数值表示模型样本与扩展样本数量之比。

图5 机构可靠度对比结果

在5 个数量关系下，样本数据完全相同，所以只对比本文ARI-RSM 模型预测结果的准确性。由图5 可以看出，在数量关系为1∶4 的情况下，计算结果与准确值有一定的出入，尤其是在磨损周期为30～50 时，与准确值有较大差距，但也能近似地看出机构可靠度变化趋势；在数量关系为2∶3、1∶1 的情况下，计算结果与准确值在较长时间内分保持贴合，且数量关系为2∶3 时贴合周期更长，误差能够保持在一定范围内；在数量关系为3∶2 及4∶1 时，两条曲线基本重合，误差非常小，可以近似代替准确值使用。综合考虑，在样本数量与预测数量之比为2∶3 时最为实用，既能保证运算精度，又能大大缩短所用时长。

4 结 论

1）本文以小样本、高效率及高精度为目的，考虑机构磨损退化故障，给出了一套完整的机构时变可靠性分析方法。此方法包含仿真与可靠度计算两部分：利用ADAMS 与MTALAB 联合仿真缩短样本获取时间；应用ARIMA 与RSM 方法，建立ARI-RSM 方法，扩充样本容量，对机构后续运转情况进行预测并计算了可靠度。

2）以曲柄滑块机构为实例，通过对比分析验证了本文方法的准确性与实用性，并确定了较为合适的样本与预测数量之比（2∶3）。

3）本方法具有一定普适性，可应用于包括其他故障的机构时变可靠性分析。