考虑竞争失效与多源不确定性的扑翼飞行器驱动机构可靠性分析

Reliability Analysis for Driving Mechanisms of Flapping-wing Aircrafts Considering Competitive Failure and Multi-source Uncertainty

摘要

针对仿生扑翼飞行器驱动机构因铰链磨损导致寿命短的缺陷，提出含间隙铰链的扑翼驱动机构动态可靠性评估方法。基于扑翼飞行器驱动机构的工作原理对机构的失效模式进行分析，确立了磨损与运动不对称性为扑翼飞行器驱动机构的关键失效模式，融合Hertz接触力模型与Archard磨损理论建立了非均匀磨损动力学模型。通过离散控制点绘制样条曲线重构磨损轮廓，实现了多周次磨损动态更新仿真。在结合Kriging和蒙特卡罗模拟的主动学习方法（AK-MCS）基础上考虑竞争失效模式，将其改进为Improved-AKMCS，以量化铰链磨损与摇臂运动不对称的竞争失效机制。结果表明：机构前65 000周次运行时以磨损失效为主，至68 000周次因运动不对称失效加剧，总可靠度由0.9989降至0.0165。灵敏度分析证实了弹性模量、磨损率及铰链初始间隙对可靠性具有显著相关性。

Abstract

Aiming at the lifespan bottlenecks of the bionic flapping-wing aircraft drive mechanisms caused by hinge wear， a dynamic reliability assessment method was proposed considering clearance hinges. Based on the working principle of the flapping-wing aircraft drive mechanisms， the failure modes of the mechanisms were analyzed， and wear and motion asymmetry were established as the key failure modes. A non-uniform wear dynamics model was developed by integrating the Hertz contact force model and the Archard wear theory. The wear profile was reconstructed by drawing spline curves through discrete control points， enabling dynamic update simulation for multiple cycles of wear. Building upon the active learning method combining Kriging and Monte Carlo simulation（AK-MCS） and considering competitive failure modes， improved-AKMCS was used to quantify the competitive failure mechanism between hinge wear and rocker motion asymmetry. The results indicate that wear failure dominates the first 65 000 operation cycles of the mechanisms； until 68 000 cycles， the failure due to motion asymmetry intensifies， and the overall reliability drops from 0.9989 to 0.0165. Sensitivity analysis confirms that the elastic modulus， wear coefficient， and initial hinge clearance have significant correlations with reliability.

Graphical abstract

关键词

扑翼驱动机构 / 失效模式 / 铰链磨损 / 机构可靠性 / 非均匀磨损

Key words

flapping-wing drive mechanism / failure mode / hinge wear / mechanism reliability / non-uniform wear

0 引言

在我国航空航天技术快速发展的背景下，仿生微型扑翼飞行器因兼具可巡航、可悬停的能力及尺寸小、机动性高、隐蔽性强的特点，在军民领域应用前景广阔^［1］，然而仿生微型扑翼飞行普遍存在的寿命短缺陷严重制约了技术发展。现有机械式驱动系统^［2］因含大量旋转铰链，加工装配中易产生初始缺陷、误差及磨损，进而诱发铰链间隙，从而导致运动轨迹偏移与精度下降，加速系统失效，故亟需开展磨损工况下可靠性退化机理研究。

铰链作为扑翼机构磨损核心部位，其研究聚焦非均匀磨损机制与动态响应。FLORES^［3］基于Archard准则建立了四连杆动力学模型，揭示了间隙运动副磨损空间的非均匀性；ZHUANG等^［4］通过自适应有限元耦合贝叶斯推断，实现了含润滑铰链磨损系数的动态更新；ZHU等^［5］提出的接触刚度算法量化了磨损对曲柄滑块机构动态响应的耦合影响。

针对实际机构多间隙铰链磨损预测需求，当前研究实现了三重突破：ERKAYA等^［6］将无质量连杆实验与轴套连续接触理论融合，实现了双铰链磨损深度精准预测；喻天翔等^［7］基于运动监测数据的综合预测框架，通过贝叶斯理论动态更新了磨损系数分布；GENG等^［8］采用B样条曲线拟合非均匀磨损轮廓，建立了多体动力学预测模型；BAI等^［9］量化揭示了间隙尺寸与铰链数量对机构振动特性的显著影响规律。

为控制机构间隙引发的运动精度衰退，可靠性评估领域发展出多尺度解决方案：庞欢等^［10］综合加工误差、变形与磨损构建体系，验证了多失效模式下系统的可靠性；ZHANG等^［11］基于径向基神经网络建立了参数-位姿映射，将蒙特卡罗与混合学习结合用来评估机器人定位可靠性；WU等^［12］提出了概率阶跃模型使理论实验误差小于5%；黄鹏^［13］将微分运动学与鞍点近似法融合，解析了混合不确定性的精度可靠性；SUN等^［14］建立了随机-认识双重不确定性的运动学精度通用框架；周长聪等^［15］通过Kriging代理与蒙特卡罗方法量化了销轴磨损对失效概率的敏感性影响。

近年来，基于Kriging代理模型的主动学习策略已成为可靠性分析新范式：结合Kriging和蒙特卡罗模拟的主动学习方法（AK-MCS）^［16］基于U学习函数自适应采样，将计算量降至传统蒙特卡罗模拟（MCS）的20%，但未考虑时变退化；用于时变可靠性分析的Kriging辅助随机过程离散化技术（K-TRPD）^［17］引入随机过程离散化处理时变可靠性，在曲柄滑块机构中实现了误差小于5%结果，但其线性假设难以描述非均匀磨损；梯度下降的自适应Kriging方法（AK-GD）^［18］通过梯度下降优化采样点选择，高维问题中样本调用量减少40%，但在多失效模式竞争中易早收敛；并行Kriging^［19］针对小失效概率问题开发了高效建模策略，在航空航天结构案例中使效率提高了30%，但缺乏动态轮廓更新机制。

现有研究主要存在以下三点不足：磨损建模单一性、失效模式孤立性、不确定性处理局限性。针对上述局限性，本文以扑翼驱动机构为研究对象，在AK-MCS方法的基础上考虑竞争失效模式，将其改进为Improved-AKMCS，建立考虑关键失效模式的扑翼驱动机构可靠性分析方法。同时，考虑构件尺寸误差与材料参数误差的影响，将其作为随机因素，建立机构可靠性功能函数，通过计算得到机构的可靠度。此外，还对机构进行了可靠度灵敏度分析，分析随机因素对扑翼驱动机构可靠性的影响。

1 扑翼驱动机构功能分析

1.1 扑翼驱动机构组成与工作原理

扑翼驱动机构作为扑翼飞行器的关键组成部分，与齿轮配合运转能够实现将电机电能转化为两侧扑翼做往复运动的机械能，是维持扑翼飞行器飞行状态不可或缺的一部分，其结构示意见图1。

本文研究的扑翼飞行器驱动机构设计和制造工艺相对简单，整个驱动机构主要包括三部分：齿轮减速器、机架和扑动机构齿轮减速器。该机构能够将由电机驱动的高转速运动转化为适合带动扑动机构的低转速运动。扑动机构的主要组成部件有曲柄轴、连杆和摇臂等。扑动机构通过固定在低速轴大齿轮上的曲柄轴和多个旋转铰链将低速轴大齿轮、连杆和摇臂组合起来形成单曲柄双摇杆机构，从而将齿轮的旋转运动转化为摇臂上下往复的扑动运动。扑翼驱动机构的机架作为一个整体，主要起固定、支撑驱动电机、齿轮减速器和扑动机构的作用，是整个机构的基础。

1.2 磨损机理与影响因素分析

扑翼驱动机构旋转铰链的失效行为与黏着磨损机制具有显著相关性。通过对接触面黏着特性等关键参数的研究，可有效建立铰链磨损的量化模型，为机构可靠性建模提供理论支撑。

2 含间隙铰链的扑翼驱动机构动力学仿真

2.1 间隙铰链建模

为了更合理地表述旋转铰链的间隙，本文引入图2所示的间隙矢量模型。假设旋转铰链由轴和轴套组成，轴的中心点为O_{\mathrm{1}}、半径为R_{\mathrm{1}}，轴套中心点为O_{\mathrm{2}}、半径为R_{\mathrm{2}}，则铰链间隙可表示为

c=R_{\mathrm{2}}-R_{\mathrm{1}}（1）

铰链间隙矢量的起点为O_{\mathrm{1}}，终点为O_{\mathrm{2}}，间隙矢量的长度被限制在以O_{\mathrm{2}}为中心、以c为半径的间隙圆内。间隙矢量\mathit{e}可由全局坐标系下的轴位置矢量{\mathit{r}}_{o{o}_{\mathrm{1}}}和轴套位置矢量{\mathit{r}}_{o{o}_{\mathrm{2}}}之差表示：

\mathit{e}={\mathit{r}}_{O{O}_{\mathrm{1}}}-{\mathit{r}}_{O{O}_{\mathrm{2}}}（2）

间隙矢量的单位长度可定义为

e=\sqrt[]{{\mathit{e}}^{\mathrm{T}}\mathit{e}}（3）

由式（1）和式（3）可以得到轴与轴套接触、碰撞引起变形的计算公式：

\delta =e-c（4）

其中，\delta为铰链元素之间的法向穿透深度，\delta的数值大小可作为轴与轴套运动状态的参考。当\delta >\mathrm{0}时，轴与轴套接触并发生弹性变形，称为碰撞状态（图3a）；当\delta =\mathrm{0}时，轴与轴套开始发生接触或开始分离（图3b）；当时，轴与轴套不产生接触，为自由状态（图3c）。

2.2 扑翼驱动机构动力学建模

2.2.1 法向接触力模型

铰链法向接触力的计算模型为Hertz模型，其公式如下：

F_{\mathrm{N}}=K_{\mathrm{n}}{\delta }^n（5）

式中：F_{\mathrm{N}}为铰链的法向接触力；K_{\mathrm{n}}为广义刚性系数；n为力系数，取决于接触元素的材料和接触区域的几何特性，对于金属接触，一般取1.5。

广义刚性系数由下式计算：

K_{\mathrm{n}}=\frac{\mathrm{4}}{\mathrm{3}\mathrm{\pi }(h_{\mathrm{1}}+h_{\mathrm{2}})}\sqrt[]{\frac{R_{\mathrm{1}}{R}_{\mathrm{2}}}{R_{\mathrm{1}}-R_{\mathrm{2}}}}（6）

h_l=\frac{\mathrm{1}-{\nu }_l^{\mathrm{2}}}{E_l}l=\mathrm{1,2}

式中：h_{\mathrm{1}}、h_{\mathrm{2}}分别为铰链轴和轴套的材料参数；{\nu }_l为材料的泊松比；E_l为材料的弹性模量。

2.2.2 切向接触力模型

铰链的切向接触力F_{\mathrm{T}}可看作是法向接触力的函数：

F_{\mathrm{T}}=\mu F_{\mathrm{N}}（7）

式中：\mu为摩擦因数。

2.3 旋转铰链磨损计算

扑翼驱动机构旋转铰链的磨损更符合黏着磨损机理。Archard模型同时考虑了滑移速度、垂直载荷、材料特性等参数对磨损量的影响，是最接近黏着磨损机理的模型之一，其计算公式如下：

h=K_{\mathrm{H}}\frac{F_{\mathrm{N}}S}{A_{\mathrm{a}}}（8）

式中：h为磨损深度；K_{\mathrm{H}}为磨损率，它与铰链元素的材料属性相关；S为铰链元素之间的相对滑移距离；A_{\mathrm{a}}为铰链元素的接触面积。

根据扑翼驱动机构的工作原理及运动规律，各铰链的实际磨损应为非均匀磨损。为了更加准确地描述铰链的非均匀磨损特性，本文采用一种基于样条曲线的非均匀磨损计算方法。该方法使用N个控制点均匀划分轴套内轮廓，并利用样条曲线将控制点连接起来，构成接触碰撞中铰链的轴套内轮廓，如图4a所示。相比于将轴套划分为多个离散区域分别与轴建立接触碰撞模型的传统方法，本文所采用的方法保持了轴承轮廓的平滑连续性，使得计算出的接触力更加精确，从而获得更加精确的铰链磨损特性。

设机构仿真总时长为T₀，将其离散为m个时间区间。对于第j（j=1，2，…，m）个时间间隔[t_j,t_j+\mathrm{\Delta }t]和第i(i=\mathrm{1,2},\dots ,N)个控制点，它在当前时间间隔内的磨损深度计算公式为

h_{j,i}=K_{\mathrm{H}}\frac{F_{j,i}{S}_{j,i}}{A_{j,i}}（9）

式中：F_j，i 为第j个时间间隔下第i个控制点所受的法向接触力；S_j，i 为第j个时间间隔下第i个控制点处铰链元素之间的相对滑移距离；A_j，i 为第j个时间间隔下第i个控制点处铰链元素的接触面积。

图4b示意了铰链轴套在运行若干周次后的非均匀磨损曲线，\stackrel{̑}{AB}表示轴套发生非均匀磨损后轴和轴套的接触曲线，由于轴套的实际磨损非常轻微，这里将磨损深度进行放大演示。假设每个控制点的累积磨损深度为h_{t,i}(i=\mathrm{1,2},\dots ,N)，则每个控制点在全局坐标系（GCS）中的坐标为

{\mathit{P}}_i=(\frac{\mathrm{2}\mathrm{\pi }}{N}(i-\mathrm{1}),R_{\mathrm{2}}+\mathrm{\Delta }{h}_{t,i})（10）

式中：t为运行时间周次；i为控制点个数；\mathrm{\Delta }{h}_{t,i}为在第t个周次运行时第i个控制点的磨损深度。

利用样条曲线模拟得到非均匀磨损后轴套轮廓的数学表达式为

\mathit{C}(\theta )={\displaystyle \sum _{i=\mathrm{1}}^N}{N}_{i,k}(\frac{\theta }{\mathrm{2}\mathrm{\pi }}){\mathit{P}}_i（11）

其中， C 为GCS中磨损轮廓边界函数；\theta为GCS中的圆周角；N_{i,k}(\theta /(\mathrm{2}\mathrm{\pi }))指第i个控制点{\mathit{P}}_i的k阶样条曲线基函数，其计算公式为

N_{i,k}(\theta /(\mathrm{2}\mathrm{\pi }))=\frac{\theta -{\theta }_i}{{\theta }_{i+k}-{\theta }_i}{N}_{i,k-\mathrm{1}}(\theta /(\mathrm{2}\mathrm{\pi }))+

\frac{{\theta }_{i+k+\mathrm{1}}-\theta }{{\theta }_{i+k+\mathrm{1}}-{\theta }_{i+\mathrm{1}}}{N}_{i+\mathrm{1},k-\mathrm{1}}(\theta /(\mathrm{2}\mathrm{\pi }))（12）

式中：{\theta }_i为全局坐标系中每个控制点的圆周角。

经过多周次的运动，铰链磨损更新完成后，可以得到轴套轮廓上每个控制点的半径计算公式：

R(\theta )=\sqrt[]{(\mathit{C}{(\theta ))}^{\mathrm{T}}\mathit{C}(\theta )}（13）

根据余弦公式，在图4b标注的\mathrm{\Delta }{O}_{\mathrm{1}}{O}_{\mathrm{2}}A和\mathrm{\Delta }{O}_{\mathrm{1}}{O}_{\mathrm{2}}B中可以得到A、B控制点与X轴的夹角分别为

{\theta }_A=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}(\frac{(\delta -R(\theta )+R_{\mathrm{1}}{)}^{\mathrm{2}}+(R({\theta }_A{))}^{\mathrm{2}}-{R_{\mathrm{1}}}^{\mathrm{2}}}{\mathrm{2}(\delta -R(\theta )+R_{\mathrm{1}})R({\theta }_A)})（14）

{\theta }_B=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}(\frac{(\delta -R(\theta )+R_{\mathrm{1}}{)}^{\mathrm{2}}+(R({\theta }_B{))}^{\mathrm{2}}-{R_{\mathrm{1}}}^{\mathrm{2}}}{\mathrm{2}(\delta -R(\theta )+R_{\mathrm{1}})R({\theta }_B)})（15）

由于轴套的磨损深度相对轴套的厚度非常小，\mathrm{\angle }A{O}_{\mathrm{2}}B也很小，因此\stackrel{̑}{AB}的长度可近似用|AB|来替代，\stackrel{̑}{AB}可以通过下式计算得到：

\stackrel{̑}{AB}=\sqrt[]{{\left|O_{\mathrm{1}}A\right|}^{\mathrm{2}}+{\left|O_{\mathrm{1}}B\right|}^{\mathrm{2}}-\mathrm{2}\left|O_{\mathrm{1}}A\right|\left|O_{\mathrm{1}}B\right|\mathrm{c}\mathrm{o}\mathrm{s}({\theta }_B-{\theta }_A)}（16）

则铰链的接触面积A_{j,i}的计算公式可写为

A_{j,i}=\stackrel{̑}{AB}·L（17）

式中：L为轴套的宽度。

在机构长期运行仿真中，单次工作循环引起的铰链磨损量级微小，每个运动周次更新磨损后轴套内轮廓的变化很小，如果要对机构进行几万周次的仿真，每一运动周次后更新磨损不仅会产生过高的计算负荷，还会形成大量过程性冗余数据。为优化计算流程，本文提出磨损外推系数T的建模方法：将单次循环的磨损量与系数T相乘，即可等效推算T次循环后的累积磨损量。这种处理方式在保证计算精度的前提下可显著降低计算资源消耗，有效解决长期仿真模拟效率低的问题。

为了更准确地模拟旋转铰链的接触碰撞过程，将建立的旋转铰链轴套内轮廓更换为样条曲线，根据磨损计算公式新建各样条点的磨损半径参数，并将新的磨损半径参数关联到对应的样条点，从而实现铰链轴套内轮廓的更新。使用样条线来模拟轴套内轮廓时，设置样条点个数为36，如图5所示。

2.4 扑翼驱动机构动力学仿真及结果

上文已建立含间隙铰链扑翼驱动机构的动力学模型，利用MATLAB编程实现了每运行完一个磨损外推周次后，输出36个样条点的磨损半径以及两摇臂角速度偏差值，然后对各铰链轴套内轮廓进行更新，继续下一磨损周次的仿真，如此循环反复直至达到规定的运行周次数。含间隙铰链的扑翼驱动机构的动力学仿真流程见图6。

扑翼驱动机构在运动过程中，由低速轴大齿轮带动曲柄轴以480 r/min的转速转动。本节将对转速为480 r/min、包含6个间隙铰链的扑翼驱动机构进行动力学仿真，研究机构运行的动力学响应与铰链磨损特性。

使用LMS Virtual. Lab软件进行动力学仿真，得到一个运动周次内的电机输出轴齿轮的驱动力矩、模拟的摇臂气动载荷如图7所示，满足机构的运行需求。

图8所示为第一个运动周次内机构的6个间隙铰链法向接触力随时间的变化趋势。从图中可以看到，所有铰链接触力的趋势都大致相同，且位于驱动机构左右对称位置的两个间隙铰链即铰链1和铰链2、铰链3和铰链4、铰链5和铰链6的接触力曲线更为接近，但铰链5和铰链6接触力曲线较其他铰链整体上凸度较大，在0.100 s附近铰链5和铰链6的接触力最大可达到9.27 N，其他铰链的接触力最大为5.94 N。这是由于气动载荷直接加载在摇臂两端，使得摇臂较其他构件受力更大、运动更加剧烈，铰链旋转产生的法向接触力也更大。虽然左右对称位置的两个间隙铰链接触力很接近，但还是存在一定差异，这是由于扑翼驱动机构采用的是单曲柄双摇杆结构，结构本身就存在一定的不对称性，经过结构优化后作为扑翼驱动机构，不对称误差在可接受范围内。

图9所示为机构在第一个运动周次内，各铰链轴与轴套相对转动的角度变化。可以看到，铰链1与铰链2的运动几乎同步，差异较小，铰链3与铰链4旋转角度有较大差异，这是由机构本身决定的，经过优化设计后得到的摇臂铰链5和铰链6的旋转角度差异减小，最大角度差不超过1°。

图10所示为机构在第一个运动周次内摇臂扑动角度变化，两摇臂扑动趋势高度同步，且扑动至最高点对应角度均为29.570°，扑动至最低点对应角度分别为-29.742°和-29.743°，这与扑翼驱动机构的实际工作情况相符合，验证了所建立的动力学模型的有效性。

3 扑翼驱动机构可靠性分析

3.1 随机因素分析

受加工、制造工艺的限制，在实际中，机构构件的不确定性无处不在，这些不确定性对机构功能有一定影响，因此，识别扑翼驱动机构的不确定性来源并对其进行适当量化是非常重要的。

3.1.1 制造公差

受机器精度的限制，零件尺寸参数的误差是不可避免的，本节将机构尺寸参数视为服从截尾正态分布的随机变量^［20］。

对于间隙铰链，轴与轴套配合产生的间隙如图11所示，用公式表达为

c=\frac{\mathrm{1}}{\mathrm{2}}(D-d)（18）

式中：D为轴套内孔直径；d为轴的直径。

间隙的均值和方差分别为

{\mu }_{\mathrm{c}}=\frac{\mathrm{1}}{\mathrm{2}}(T_D+T_d-\mathrm{2}{e}_{\mathrm{s}})（19）

{\sigma }_{\mathrm{c}}^{\mathrm{2}}=\frac{\mathrm{1}}{\mathrm{144}}(T_D^{\mathrm{2}}+T_d^{\mathrm{2}})（20）

其中，T_D为轴套直径公差，T_d为轴直径公差，e_{\mathrm{s}}为轴的上偏差，根据实际需要可参考ISO 286-2。

根据图12，构件的长度误差可表示为

T_{\mathrm{L}}=E_{\mathrm{S}}-E_{\mathrm{I}}（21）

式中：E_{\mathrm{S}}为构件长度上偏差；E_{\mathrm{I}}为构件长度下偏差。

假设构件长度公差在零线附近对称，根据\mathrm{6}\sigma原则^［21］，公差的均值和方差可表示为

{\mu }_{\mathrm{L}}=\frac{\mathrm{1}}{\mathrm{2}}\left|E_{\mathrm{S}}+E_{\mathrm{I}}\right|=\mathrm{0}（22）

{\sigma }_{\mathrm{L}}=\frac{\mathrm{1}}{\mathrm{6}}{T}_{\mathrm{L}}（23）

结合扑翼驱动机构的工作原理及构件的加工方式，最终确定尺寸随机变量为轴套宽度、铰链间隙、连杆长度、连杆倾斜加工角和摇臂长度，具体分布见表1。

3.1.2 材料参数误差

材料制造过程中的不确定性会导致材料的物理性质在一定范围内发生变化^［22］，因此铰链材料参数的随机分布设定基于多源实验数据与物理约束。由于材料参数误差波动符合制造批次统计规律，因此将其视为服从截尾正态分布的随机变量进行分析^［23］，但受物理边界限制，需要进行截断。此外参数均值均取自实测数据，标准差通过多组试样统计。在铰链磨损计算公式中，轴和轴套的材料属性参数主要为弹性模量、泊松比和磨损率这三个参数，具体分布如表2所示。

除了上述尺寸变量及材料特性变量外，在扑翼驱动机构的动力学仿真中，还考虑了模拟气动载荷弹簧力的弹性系数以及接触碰撞模型中的摩擦因数对机构的影响，同样将其视为服从截尾正态分布的随机变量，具体分布见表3。

3.2 可靠性功能函数

考虑不确定性因素的影响，扑翼驱动机构主要失效模式为铰链磨损失效和扑动运动不对称失效。磨损量一般以线性磨损量表示居多，通常指接触面上的法向磨损量。针对铰链磨损失效这一失效模式，当轴套的法向磨损深度超过最大允许磨损深度时，铰链视为失效，不能完成规定的任务，其失效函数为

G_{\mathrm{w}}={\delta }_{\mathrm{0}}(x_{\mathrm{1}},x_{\mathrm{2}},\dots ,x_n)-{\delta }_{\mathrm{m}\mathrm{a}\mathrm{x}}（24）

式中：{\delta }_{\mathrm{0}}(x_{\mathrm{1}},x_{\mathrm{2}},\dots ,x_n)为考虑不确定性因素的机构铰链磨损最大值；{\delta }_{\mathrm{m}\mathrm{a}\mathrm{x}}为最大允许磨损深度。

根据文献［24］，定义机构铰链最大允许磨损量{\delta }_{\mathrm{m}\mathrm{a}\mathrm{x}}=\mathrm{0.009}~\mathrm{0.018}D。在没有详细数据的情况下，建议将铰链的最大允许磨损量取为{\mu }_{{\delta }_{\mathrm{m}\mathrm{a}\mathrm{x}}}\approx \mathrm{0.01}D，{\sigma }_{{\delta }_{\mathrm{m}\mathrm{a}\mathrm{x}}}\approx \mathrm{0.0002}~\mathrm{0.0005}D。扑翼驱动机构各铰链轴套的标称直径均为3 mm，因此对应的轴套最大允许磨损量为0.03 mm，所以铰链磨损失效的失效函数可写为

G_{\mathrm{w}}={\delta }_{\mathrm{0}}(x_{\mathrm{1}},x_{\mathrm{2}}\dots ,x_n)-\mathrm{0.03}\mathrm{m}\mathrm{m}（25）

磨损失效的失效概率可表示为

（26）

针对扑翼运动不对称失效这一失效模式，由实际试验及外场数据得到当扑翼驱动机构两摇臂的扑动角度偏差约1°时，视为机构失效。在LMS Virtual. Lab软件中进行仿真模拟，得到当两侧摇臂扑动角度偏差超过1°时，两摇臂角速度的最大偏差约为25 r/min，如图13所示。

当扑翼驱动机构两摇臂的扑动角速度最大偏差超过25 r/min时，视为机构失效。则扑动运动不对称失效的失效函数为

G_{\mathrm{s}}={\gamma }_{\mathrm{m}\mathrm{a}\mathrm{x}}-\mathrm{25}\mathrm{r}/\mathrm{m}\mathrm{i}\mathrm{n}（27）

其中，{\gamma }_{\mathrm{m}\mathrm{a}\mathrm{x}}表示每个运动周次内两摇臂扑动角速度的最大偏差，扑动运动不对称的失效概率可表示为

（28）

扑翼驱动机构的铰链磨损失效属于零部件级失效，扑翼运动不对称属于系统级失效。在计算驱动机构可靠性时，先根据式（26）判断零部件失效模式即铰链是否磨损失效，若失效，则认为机构已经失效；若未失效，则继续根据式（28），即机构是否发生运动不对称失效，若系统级也失效，则认为机构失效，否则认为机构性能正常，仍能正常工作，机构的可靠性功能函数为

（29）

3.3 基于Improved-AK-MCS算法和动力学仿真分析的可靠性分析

3.3.1 Improved-AK-MCS改进机制

本文在传统AK-MCS方法^［16］基础上提出以下三重创新性改进机制。

1）双失效域协同采样机制。针对铰链磨损失效（G_{\mathrm{w}}）与运动不对称失效（G_{\mathrm{s}}）的竞争特性，创新性地构建耦合学习函数U_{\mathrm{c}}(x)：

U_{\mathrm{c}}(x)=\mathrm{m}\mathrm{i}\mathrm{n}(\frac{|{\mu }_{G_{\mathrm{w}}}(x)|}{{\sigma }_{G_{\mathrm{w}}}(x)},\frac{|{\mu }_{G_{\mathrm{s}}}(x)|}{{\sigma }_{G_{\mathrm{s}}}(x)})（30）

式中：x为扑翼驱动机构多源不确定性的随机变量；G_{\mathrm{w}}(x)为铰链磨损失效功能函数，当G_{\mathrm{w}}(x)\le \mathrm{0}时，机构因铰链磨损量超过阈值而失效；G_{\mathrm{s}}(x)为运动不对称失效功能函数，当G_{\mathrm{s}}(x)\le \mathrm{0}时，机构因运动轨迹不对称度超过允许范围而失效；{\mu }_{G_{\mathrm{w}}}(x)、{\mu }_{G_{\mathrm{s}}}(x)表示功能函数的预测均值；{\sigma }_{G_{\mathrm{w}}}(x)、{\sigma }_{G_{\mathrm{s}}}(x)表示功能函数的预测标准差。

该函数突破传统单失效域采样局限^［1］，通过同步评估两失效域的边界概率密度实现关键区域协同优化采样。

2）动态权重调整机制。为表征失效模式竞争关系，引入权重因子k：

k=\frac{\partial {\gamma }_{\mathrm{m}\mathrm{a}\mathrm{x}}}{\partial {\delta }_{\mathrm{0}}}\frac{{\sigma }_{G_{\mathrm{w}}}}{{\sigma }_{G_{\mathrm{s}}}}（31）

当|k|>\mathrm{0.5}时，表明以磨损失效为主导，系统自动提高磨损域采样密度。该机制会精准响应机构运行不同阶段的失效模式转换特征。

3）收敛准则强化机制。新增相对误差判据：

（32）

其中，R_k和R_{k-\mathrm{1000}}分别为当前和1000周次前的系统可靠度估计值。式（32）通过比较二者间的相对误差，构成一个严格的收敛判据，确保可靠性分析结果能反映系统的长期退化行为，避免过早停止计算，确保68 000周次后的可靠度预测误差小于3%。

3.3.2 Improved-AKMCS仿真分析流程

结合上述内容及扑翼驱动机构失效判据，可以得到磨损可靠度与不对称性可靠度为两种失效模式，其竞争的结果为总可靠度。基于考虑竞争失效的改进AK-MCS算法（improved-AKMCS）和LMS Virtual. Lab仿真分析的可靠性分析流程见图14，其主要步骤如下：

1）对3.1节中提到的30个随机变量进行蒙特卡罗抽样法抽取100 000个样本点，构建MC样本点集，再通过拉丁超立方抽样法在MC样本点集中选择200个样本点组成初始设计样本点集。

2）在LMS Virtual. Lab软件中根据DOE样本点集更改模型参数并进行仿真计算，得到DOE样本点集的铰链磨损预测值和摇臂角速度偏差，并构建初始铰链磨损Kriging模型和初始运动不对称Kriging模型，然后使用初始Kriging模型分别计算MC样本点集中所有样本点的铰链磨损预测值{\widehat{y}}_{\mathrm{w}}(x)和摇臂角速度偏差{\widehat{y}}_{\mathrm{s}}(x)及对应的方差{\sigma }_{{\widehat{y}}_{\mathrm{w}}}^{\mathrm{2}}和{\sigma }_{{\widehat{y}}_{{}_{\mathrm{s}}}}^{\mathrm{2}}。

3）在两个Kriging模型中，根据U学习函数分别选出MC样本点中最好的点，即U学习函数值最小的样本点。若其学习函数值小于2，则将其加入DOE样本点中重新构建Kriging模型，得到新的样本点预测值；若学习函数值大于等于2，则学习停止。

4）根据机构的失效判据和失效概率的变异系数判断机构的失效情况，并计算失效概率{\widehat{P}}_{\mathrm{f}}和变异系数{\delta }_{{\widehat{P}}_{\mathrm{f}}}。若得到的变异系数大于等于5%，则需要更新DOE样本点集，重新开始步骤1）；若变异系数小于5%，则满足条件，进行下一步。

5）输出精确的Kriging模型和失效概率。

3.3.3 Kriging模型参数设定

基于高维非线性特征，配置如下参数：

1）相关函数。采用高斯核函数适配磨损深度{\delta }_{\mathrm{0}}与角速度偏差{\gamma }_{\mathrm{m}\mathrm{a}\mathrm{x}}的非线性映射：

R(x^{(i_{\mathrm{1}})},x^{(i_{\mathrm{2}})};\theta )=\mathrm{e}\mathrm{x}\mathrm{p}(-{\displaystyle \sum _{m=\mathrm{1}}^{\mathrm{30}}}{\theta }_m|x_m^{(i_{\mathrm{1}})}-x_m^{(i_{\mathrm{2}})}{|}^{\mathrm{2}})（33）

式中：x^{(i_{\mathrm{1}})}、x^{(i_{\mathrm{2}})}为随机变量x的样本点，i_{\mathrm{1}}、i_{\mathrm{2}}代表从随机变量集合中抽取或生成的第i_{\mathrm{1}}个和第i_{\mathrm{2}}个样本点（i_{\mathrm{1}}\mathrm{、}{i}_{\mathrm{2}}=\mathrm{1,2},\dots ,N）；m为灵敏度分析参数编号。

2）超参数优化。通过最大似然估计（MLE）求解{\theta }_k^{\mathrm{*}}：

{\theta }_k^{\mathrm{*}}=\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{a}{\mathrm{x}}_{{\theta }_k}[-\frac{n}{\mathrm{2}}\mathrm{l}\mathrm{n}({\widehat{\sigma }}^{\mathrm{2}})-\frac{\mathrm{1}}{\mathrm{2}}\mathrm{l}\mathrm{n}\left|\mathit{R}\right|]（34）

式中：n为采用拉丁超立方抽样的初始样本，n=200； R 为Kriging模型的相关矩阵；\widehat{\sigma }为样本标准差。

3）回归模型。常数基函数\mathit{f}(x)=[\mathrm{1}]，方差估计为

{\widehat{\sigma }}^{\mathrm{2}}=\frac{\mathrm{1}}{n}{(\mathit{y}-\mathit{F}\widehat{\mathit{\beta }})}^{\mathrm{T}}{\mathit{R}}^{-\mathrm{1}}(\mathit{y}-\mathit{F}\widehat{\mathit{\beta }})（35）

式中： y 为仿真得到的输出值，包括铰链最大磨损值{\delta }_{\mathrm{0}}与摇臂角速度偏差{\gamma }_{\mathrm{m}\mathrm{a}\mathrm{x}}； F 为基函数矩阵；\widehat{\mathit{\beta }}为回归系数向量。

3.4 可靠性分析结果

对扑翼驱动机构进行72 000周次仿真，设置磨损外推系数为1000，采用3.3节所述方法，结合可靠性功能函数分别计算仅考虑磨损失效、仅考虑机构不对称失效以及同时考虑两种失效模式三种情况下的机构可靠性，即考虑竞争失效模式的可靠性分析。

3.4.1 考虑磨损的机构可靠性

图15和图16显示了机构运行64 000周次、66 000周次和68 000周次后铰链最大磨损值的概率分布情况，可以看到铰链的最大磨损值基本符合正态分布。可以观察到随着机构运行周次的增加，铰链最大磨损值概率分布峰值对应的横坐标逐渐向右移动，也就意味着铰链磨损的均值稳步增大，机构运行64 000周次后铰链最大磨损值分布均值约为0.0288 mm，运行66 000周次后最大磨损值分布均值约为0.0302 mm，运行68 000周次后最大磨损值分布均值约为0.0313 mm。

图17和图18分别为铰链最大磨损值的累积分布和考虑磨损的机构可靠度曲线，可以看到机构的磨损可靠度在运行63 000周次后开始从1逐渐下降，机构运行64 000周次的磨损可靠度为0.9652，运行66 000周次的磨损可靠度为0.4516，运行68 000周次的磨损可靠度为0.0357，直至运行70 000周次后机构磨损可靠度降为0。

3.4.2 考虑不对称失效的机构可靠性

图19所示为机构运行64 000周次、66 000周次和68 000周次后两摇臂角速度最大偏差值的概率分布情况，两摇臂角速度的最大偏差值也基本符合正态分布。机构运行66 000周次产生的摇臂角速度最大偏差的均值比运行64 000周次的最大偏差值略大，但远小于运行68 000周次产生的最大偏差值。机构运行64 000周次、66 000周次以及68 000周次后产生的最大偏差均值分别为21.75 r/min，24.14 r/min以及28.22 r/min，如图20所示。

图21和图22分别为角速度的最大偏差值的累积分布和考虑扑动不对称的机构可靠度曲线，机构的不对称可靠度在运行64 000周次后开始下降，机构运行64 000周次的不对称可靠度为0.9989，运行66 000周次的不对称可靠度为0.7690，运行68 000周次的不对称可靠度为0.0165，随后降为零。

3.4.3 考虑磨损和不对称的机构可靠性

根据第3.3节，同时考虑磨损失效和运动不对称失效的机构可靠性计算方法，计算得到机构可靠度结果如表4和图23所示。可以得到机构运行直到65 000周次时，同时考虑磨损和机构不对称的机构可靠度与磨损可靠度基本一致，说明在前65 000周次内磨损失效为机构的主要失效因素。随后受到机构不对称失效的影响，机构总可靠度与磨损可靠度之间的差距逐渐增大，在机构运行67 000周次时机构总可靠度曲线与不对称可靠度曲线之间的差距减小，反映了机构运行后期摇臂不对称失效对机构可靠度的影响更大。

3.5 可靠度灵敏度分析

在3.3节采用Improved-AKMCS方法预测扑翼驱动机构的铰链磨损值和摇臂角速度偏差值时，输入了30个随机变量，将其逐一进行编号，如表5所示。

3.5.1 铰链磨损失效的影响因素灵敏度分析

图24a所示为铰链磨损失效参数的均值灵敏度分析结果，从中可以得到铰链1（2）轴的弹性模量、铰链1（2）的磨损率、轴套1的初始半径、轴套2的初始半径、连杆总长度以及摇臂前端长度这6个随机变量均值对铰链磨损结果的影响较大，其中铰链1（2）的磨损率、轴套1和轴套2的初始半径均值、连杆总长度与铰链最大磨损值成正相关，即这4个参数值越大，铰链的最大磨损值越大，其余两个参数的均值与铰链最大磨损值成负相关，即这两个参数均值越小，铰链的最大磨损值越大。

图24b所示为铰链磨损失效参数的标准差灵敏度分析结果，得到标准差对铰链最大磨损值影响较大的参数与均值对磨损值影响较大的参数相同，与之不同的是这6个参数的标准差都与铰链最大磨损值成负相关，想要得到较小的铰链磨损值，应该将这6个参数设置较大的标准差。

3.5.2 摇臂角速度不对称失效的影响因素灵敏度分析

图25a所示为铰链摇臂角速度不对称失效参数的均值灵敏度分析结果，从图中可以得到轴套1、轴套2、轴套5和轴套6的初始半径，连杆总长度以及摇臂前端长度这6个参数的均值对摇臂角速度不对称性影响较大，其中连杆总长度与摇臂角速度不对称性成负相关，与其余参数成正相关，即连杆总长度的均值越小或其余5个参数均值越大，摇臂不对称性越大。图25b所示为铰链摇臂角速度不对称失效参数的标准差灵敏度分析结果，得到轴套5和轴套6的初始半径以及连杆总长度对摇臂角速度不对称性影响较大，且都成负相关。

4 实验验证

4.1 实验描述

为了验证理论的正确性，使用4架微型扑翼仿生飞行器在风洞中进行实验，获得了铰链间隙的实际磨损特性。整个实验在室温下进行，研究了微型扑翼飞行器在稳定来流下从正常运行到机构失效的间隙接头磨损情况。轴颈的材料为不锈钢，轴承的材料为钛合金。实际实验中扑翼机构的尺寸和质量参数均与模拟实验相同。

实验过程如下：①将扑翼驱动机构与其专用的机构一起固定在风洞实验平台上，调整迎角为10°，设置风洞来流速度为20 m/s；②等待风洞内的风速稳定后，再激活扑翼飞行器，采用直流电压为扑翼飞行器供电，并在实验中以500 r/min的速度稳定运行；③采用高清摄像机实时监控机构运行情况；④每隔一段时间后使用放大镜和显微镜观察铰链磨损情况，直至机构失效。

实验现场如图26a所示，当机构运行失效时，关节A处出现了一个可见间隙，如图26b所示。

4.2 实验结果

将内轴承的最左点作为磨损起点，标记为0°，然后将内轴承顺时针标记至360°。轴承的磨损深度在110°处达到第一个峰值，最大磨损深度为0.0249 mm，约为轴承半径的1.66%。此后，磨损深度保持在一个低值，直到在310°处达到第二个峰值0.0087 mm。该机构运行72 000次后的轴承轮廓如图27所示，它表明磨损深度集中在轴承的顶部和底部顶点。图28所示为机构运行10.4万次循环失效后的轴承，红线为轴承的初始轮廓。可以观察到，4个实验部件的轴承磨损轮廓几乎是均匀的。图27和图28的对比表明，仿真得到的磨损结果与实际实验轴承磨损趋势一致。

由上述实验结果可得铰链在运行64 000次、66 000次、68 000次后的磨损深度概率密度函数曲线（图29）。显然，随着运行次数的增加，磨损深度的最大概率值逐渐增大，从0.028 mm逐渐扩展至0.032 mm。

4.3 实验结果与仿真结果对比分析

将实验所得的磨损深度值与仿真所得的磨损深度值进行比较，可得仿真结果与实际实验结果有较高拟合程度，如图30所示。

由实验可得，可靠度随运行周次的变化如图31所示，与仿真所得的可靠度有较高拟合度。

通过Improved-AKMCS主动学习过程，铰链磨损功能函数（G_{\mathrm{w}}）初始样本点数为200，主动学习增加点数737，最终样本点数为937；运动不对称功能函数（G_{\mathrm{s}}）初始样本点数为200，主动学习增加点数692，最终样本点数为892。关键样本均集中在失效边界（）附近。

5 结论

本文针对扑翼飞行器核心部件——扑翼驱动机构的可靠性问题，基于LMS Virtual.Lab建立了考虑间隙铰链非均匀磨损的动力学仿真模型。核心创新在于：提出了融合铰链磨损与摇臂扑动角速度偏差的竞争失效可靠性建模方法，并据此改进主动学习Kriging模型，提出Improved-AKMCS算法以实现多失效模式竞争分析。研究结果表明，机构运行前65 000周次以磨损失效为主导；后续阶段（如67 000周次后）摇臂角速度不对称失效影响显著增强，成为可靠性下降主因。灵敏度分析揭示了关键影响参数：铰链最大磨损值主要受材料属性（弹性模量）、磨损特性（磨损率）及几何尺寸（轴套半径、连杆/摇臂长度）影响；角速度不对称性则对铰链初始间隙（轴套半径）及连杆长度更为敏感。最终，通过实验对比验证了理论与模型预测的准确性。

参考文献

[1] ZHANG Y， WANG Z， ZHENG H. Aerodynamic Analysis and Design Optimization of a Novel Flapping Wing Micro Air Vehicle in Hovering Flight［J］. Journal of Applied Fluid Mechanics， 2023， 16（3）： 505-517.

[2] 赵晓伟， 曾东鸿， 占英， 等. 仿鸟类扑翼飞行器研究进展［J］. 动力学与控制学报， 2024， 22（4）： 1-15.

[3] ZHAO Xiaowei， ZENG Donghong， ZHAN Ying， et al. Advances on System Design of Bio-inspired Flapping Wing MAV［J］. Journal of Dynamics and Control， 2024， 22（4）： 1-15.

[4] FLORES P. Modeling and Simulation of Wear in Revolute Clearance Joints in Multibody Systems［J］. Mechanism and Machine Theory， 2009， 44（6）： 1211-1222.

[5] ZHUANG Xinchen， SARAYGORD AFSHARI S， YU Tianxiang， et al. A Hybrid Model for Wear Prediction of a Single Revolute Joint Considering a Time-varying Lubrication Condition［J］. Wear， 2020， 442/443： 203124.

[6] ZHU Aibin， HE Shengli， ZOU Chao， et al. The Effect Analysis of Contact Stiffness on Wear of Clearance Joint［J］. Journal of Tribology， 2017， 139（3）： 031403.

[7] ERKAYA S， UZMAY I. Experimental Investigation of Joint Clearance Effects Onthedynamics of a Slider-crank Mechanism［J］. Multibody System Dynamics， 2010， 24（1）： 81-102.

[8] 喻天翔， 庄新臣， 宋笔锋， 等. 飞机连杆机构多铰链磨损寿命综合预测方法［J］. 航空学报， 2022， 43（8）： 127-136.

[9] YU Tianxiang， ZHUANG Xinchen， SONG Bifeng， et al. Integrated Wear Life Prediction Method of Multiple Joints in an Aircraft Linkage Mechanism［J］. Acta Aeronautica et Astronautica Sinica， 2022， 43（8）： 127-136.

[10] GENG Xinyu， LI Meng， LIU Yufei， et al. Non-probabilistic Kinematic Reliability Analysis of Planar Mechanisms with Non-uniform Revolute Clearance Joints［J］. Mechanism and Machine Theory， 2019， 140： 413-433.

[11] BAI Zhengfeng， SUN Yi. A Study on Dynamics of Planar Multibody Mechanical Systems with Multiple Revolute Clearance Joints［J］. European Journal of Mechanics-A， 2016， 60： 95-111.

[12] 庞欢， 喻天翔， 宋笔锋. 平面连杆机构运动精度可靠性及灵敏度分析［J］. 中国机械工程， 2014， 25（18）： 2415-2420.

[13] PANG Huan， YU Tianxiang， SONG Bifeng. Analyses of Kinematic Accuracy Reliability and Sensitivity for a Planar Linkage Mechanism［J］. China Mechanical Engineering， 2014， 25（18）： 2415-2420.

[14] ZHANG Dequan， ZHANG Ning， YE Nan， et al. Hybrid Learning Algorithm of Radial Basis Function Networks for Reliability Analysis［J］. IEEE Transactions on Reliability， 2021， 70（3）： 887-900.

[15] WU Jianing， YAN Shaoze， ZUO M J. Evaluating the Reliability of Multi-body Mechanisms： a Method Considering the Uncertainties of Dynamic Performance［J］. Reliability Engineering & System Safety， 2016， 149： 96-106.

[16] 黄鹏. 混合不确定性下工业机器人运动精度可靠性分析与优化设计［D］. 成都： 电子科技大学， 2021.

[17] HUANG Peng. Reliability Analysis and Optimization Design of Industrial Robot Kinematic Accuracy under Mixed Uncertainties［D］. Chengdu： University of Electronic Science and Technology of China， 2021.

[18] SUN Dongyang， CHEN Guoping. Kinematic Accuracy Analysis of Planar Mechanisms with Clearance Involving Random and Epistemic Uncertainty［J］. European Journal of Mechanics - A， 2016， 58： 256-261.

[19] 周长聪， 赵浩东， 常琦， 等. 飞机舱门泄压阀机构磨损可靠性与灵敏度分析［J］. 北京航空航天大学学报， 2021， 47（4）： 690-697.

[20] ZHOU Changcong， ZHAO Haodong， CHANG Qi， et al. Reliability and Sensitivity Analysis of Relief Valve Mechanism of Aircraft Door Considering Wear［J］. Journal of Beijing University of Aeronautics and Astronautics， 2021， 47（4）： 690-697.

[21] ECHARD B， GAYTON N， LEMAIRE M. AK-MCS： an Active Learning Reliability Method Combining Kriging and Monte Carlo Simulation［J］. Structural Safety， 2011， 33（2）： 145-154.

[22] ZHANG Dequan， ZHOU Pengfei， JIANG Chen， et al. A Stochastic Process Discretization Method Combing Active Learning Kriging Model for Efficient Time-variant Reliability Analysis［J］. Computer Methods in Applied Mechanics and Engineering， 2021， 384： 113990.

[23] LIU Yushan， LI Luyi， CHANG Zeming， et al. A Novel Optimization-based Adaptive Kriging Method for Structural Reliability Analysis［J］. Engineering with Computers， 2025， 41（5）： 2953-2967.

[24] ZHAN Hongyou， XIAO Ningcong， JI Yuxiang. An Adaptive Parallel Learning Dependent Kriging Model for Small Failure Probability Problems［J］. Reliability Engineering and System Safety， 2022， 222（C）：108403.

[25] WANG Jinge， ZHANG Junfu， DU Xiaoping. Hybrid Dimension Reduction for Mechanism Reliability Analysis with Random Joint Clearances［J］. Mechanism and Machine Theory， 2011， 46（10）： 1396-1410.

[26] DONG Zhenpeng， LIANG Xiao， HOU Yuanhang， et al. 6σ Uncertainty Optimisation Strategy for USV Navigation Performance Based on Reliability Analysis［J］. Ocean Engineering， 2022， 249： 110915.

[27] MICHAEL Ashby， KARA Johnson. Materials and Design： The Art and Science of Material Selection in Product Design［M］. London： Butterworth-Heinemann， 2014.

[28] BondAndrew， HarrisAndrew. Decoding Eurocode 7［M］. London： CRC Press， 2006.

[29] 庄表中. 结构随机振动［M］. 北京： 国防工业出版社， 1995.

[30] ZHUANG Biaozhong. Random Vibration of Structure［M］. Beijing： National Defense Industry Press， 1995.

基金资助

国家自然科学基金(52202507)

陕西省自然科学基金(2024JC-YBMS-371)

陕西省自然科学基础研究计划(2024JC-YBQN-0564)

中央高校基本科研业务费专项资金(300102224204)

中央高校基本科研业务费专项资金(300102225104)