CN105205276A - 关节轴承磨损失效物理建模与分析方法 - Google Patents

Abstract

本发明属于机械工程技术领域，公开了一种关节轴承磨损失效物理建模与分析方法，包括：建立关节轴承失效物理模型、获取关节轴承结构和动态磨损数据、动态磨损过程拐点辨识、磨损阶段磨损常数辨识以及非线性磨损过程失效物理描述。本发明首次提出了融合关节轴承的磨损特性参数、结构参数、材料参数和载荷工况参数的磨损模型，参数之间的物理关系明确，能够十分准确地描述自润滑关节轴承起始于微小间隙的磨损退化过程，并用于关节轴承的结构优化设计和磨损分析中。

Description

关节轴承磨损失效物理建模与分析方法

技术领域

本发明涉及基于失效物理的关节轴承非线性磨损建模与分析方法，属于机械工程技术领域。

背景技术

关节轴承是一种特殊结构的滑动轴承，其结构简单，主要由一个有外球面的内圈和一个有内球面的外圈组成，自润滑关节轴承有衬垫。自润滑关节轴承运动时，衬垫层与内圈或外圈相对滑动而产生摩擦，其摩擦系数较小，工作时无需补充润滑剂，因而成为高端关节轴承的首选结构和润滑方案，广泛应用于航空、核电、舰船等装备的操纵与传动系统中。关节轴承的主要失效模式为磨损，其磨损分析方法对自润滑关节轴承结构设计和寿命预测具有至关重要的作用。

磨损是一个非线性过程，如图1所示，可分为三个阶段：磨合磨损阶段(I)，稳定磨损阶段(II)，急剧磨损阶段(III)。在磨合磨损阶段，磨损率随时间增加而逐渐降低；在稳定磨损阶段，摩擦表面经磨合以后达到稳定运行状态，磨损率基本不变，这是摩擦副的正常工作时期，也是决定摩擦副磨损寿命的最关键的时期；急剧磨损阶段，磨损率随时间迅速增加，使摩擦副工作条件急剧恶化，直至零件完全失效。关节轴承聚四氟乙烯(PTFE)织物衬垫的磨损过程也存在三个磨损阶段。

失效物理方法已经成功应用在机械、民用产品和航空航天结构中。

发明内容

针对关节轴承结构设计和磨损分析的需要，本发明提出一种关节轴承磨损失效物理建模与分析方法，是一种关节轴承非线性磨损过程的失效物理建模与关节轴承非线性磨损过程的分析方法，其中，磨损的分析方法融入建模的过程之中，关节轴承的失效物理模型即为磨损分析的核心算法。

本发明的技术方案是：

一种关节轴承磨损失效物理建模与分析方法，其特征在于，包括以下步骤：

步骤1：确定材料磨损模型

在自润滑关节轴承中，内环为钢材料，减磨的自润滑衬垫层为高分子复合材料，磨损主要发生在衬垫层上。关节轴承的寿命取决于自润滑衬垫层的磨损寿命。

磨粒磨损的Archard磨损公式为

$V=k_s\frac{F_N}{H}x---\left(1\right)$

黏着磨损的Archard磨损公式为

$V=k_s\frac{F_N}{3{\mathrm{\σ}}_s}x---\left(2\right)$

其中，V为磨损量，F_N为法向载荷，x为滑动距离，H为摩擦副中较软材料的硬度，σ_s为接触摩擦副中较软材料的受压屈服极限(对于自润滑关节轴承，较软材料为自润滑衬垫层)，k_s称为磨损常数。磨损常数取决于粗糙表面接触条件，如材料和润滑，分为磨粒磨损常数和黏着磨损常数。Archard滑动磨损计算方法认为磨损量与法向载荷和滑动距离成正比，与较软材料硬度(或受压屈服极限)成反比。在自润滑关节轴承中，较软材料为自润滑衬垫层。

由式(1)、(2)可知，磨粒磨损与黏着磨损的Archard公式具有相同的表达形式。自润滑关节轴承衬垫的磨损同时存在磨粒磨损和黏着磨损过程，在磨损计算时难以将二者分开。

在步骤一中，本发明将强度作为衡量自润滑衬垫层耐磨损性能的量，所以自润滑关节轴承磨损公式为

$V=k_s\frac{F_N}{{\mathrm{\σ}}_s}x---\left(3\right)$

实际使用中将关节轴承内外球面间的最大允许间隙作为其磨损失效阈值。关节轴承的结构间隙s定义为内外球面间的间隙；初始间隙u₀为产品出厂时的内外球面间隙。结构间隙为初始间隙与磨损深度u的叠加，s＝u₀+u。接触分为非协调接触与协调接触，当接触面的尺寸相对于接触物体的尺寸不可忽略时定义为协调接触，关节轴承是典型的球面协调接触情形，其关节轴承接触压力采用公式(8)和(9)计算。在关节轴承接触压力计算公式(8)和(9)中，R₁为内球面半径，R₂为外球面半径。在关节轴承磨损分析中，R₁和R₂为相对值，在磨损的过程中认为R₂＝R＝d_k/2为恒定值，而R₁随着磨损深度而变化；其中，d_k为关节轴承协调接触表面的直径。即有

R₁＝R-(u₀+u)/2(4)

ΔR＝s/2＝(u₀+u)/2(5)

由于磨损量与载荷成正比，所以最大磨损量取决于接触载荷最大的点。关节轴承中的接触应力最大点p₀处于接触区域的中心点O，接触压力随着磨损深度的变化而变化。在接触区域的中心，取一个微元区，其面积为A_a。设这个微元区为平面，且微元区上的接触压力均匀，都是p₀，所以，

F_N＝p₀A_a,V＝uA_a(6)

由式(3)可得用磨损深度表示的Archard磨损模型为

$u=k_s\frac{p_0\left(u\right)}{{\mathrm{\σ}}_s}x---\left(7\right)$

其中，p₀(u)为当磨损深度为u时的最大关节轴承接触压力，关节轴承接触压力的计算方法如步骤2所示。

步骤2：确定关节轴承接触压力计算方法

关节轴承接触压力计算方法如下，

(1)当0＜a＜h时，关节轴承中的接触压力分布采用完整球面协调接触的模型计算，即

$\left\{\begin{array}{cc}{p}_0=(n+1)\frac{F}{{\mathrm{\πa}}^2}& \left(a\right)\\ a={\[\frac{4{\mathrm{BR}}_1{R}_{2}F}{{\mathrm{\π}}^2{E}^{*}(g\mathrm{\Δ}R+c)}(n+1/2)(n+1)\]}^{1/3}& \left(b\right)\end{array}\right.---\left(8\right)$

(2)当h＜a≤R₂时，其求解方法为

$\left\{\begin{array}{cc}Q={\∫}_h^{a}r\·{(1-r^2/a^2)}^n\·arccos\frac{h}{r}\·dr& \left(c\right)\\ F_0=\frac{4(n+1)\·Q\·F}{{\mathrm{\πa}}^2-4(n+1)\·Q},F_t=F+F_0& \left(d\right)\\ p_0=(n+1)\frac{F_t}{{\mathrm{\πa}}^2}& \left(e\right)\\ a={\[\frac{4{\mathrm{BR}}_1{R}_{2t}F}{{\mathrm{\π}}^2{E}^{*}(g\mathrm{\Δ}R+c)}(n+1/2)(n+1)\]}^{1/3}& \left(f\right)\end{array}\right.---\left(9\right)$

对于以上两式，都有

$\left\{\begin{array}{cc}n=0.5-0.24\exp \[-15.08(1-a/R_2)\]& \left(g\right)\\ B=\frac{\sqrt{\mathrm{\π}}\mathrm{\Γ}(n+1)}{2\mathrm{\Γ}(3/2+n)},c=\frac{3.8304{\mathrm{BF}}_t}{{\mathrm{\π}}^2{E}^{*}{R}_2}& \left(h\right)\\ g\left(a\right)=2/\mathrm{\π}+{(a/R_2)}^2& \left(j\right)\end{array}\right.---\left(10\right)$

其中：r为球面上的点到接触区域中心的投影距离；a为接触区域的投影半径；h为关节轴承半宽；p₀为最大关节轴承接触压力；F为关节轴承承受的外部施加载荷；n为压力分布指数；F₀为等效附加载荷；Γ为gamma函数；B，c为中间参数，用公式(10h)定义。E^*为等效弹性模量，设其中，E₁,E₂分别为关节轴承内圈与外圈材料的弹性模量，μ₁,μ₂为相应的泊松比。

步骤3：确定动态磨损量表达式和数值求解过程

采用关节轴承接触压力计算公式(8)和(9)可计算出每个磨损深度u对应的p₀(u)值。如果最大接触压力p₀恒定，那么磨损深度与接触压力成正比关系。但是在关节轴承的寿命周期中，磨损量逐渐增大会导致最大接触压力逐渐增大。但在一个微小的相对滑动距离dx内，其最大接触压力可认为是恒定的，那么

$du=k_s\frac{p_0}{{\mathrm{\σ}}_s}\·dx---\left(11\right)$

du为滑动距离dx内磨损量的增量。

对于关节轴承，其工作过程中的运动方式一般为摆动，设摆动角度为±α，摆动频率为f_s。所以，在时间dt内，关节轴承内外环之间的相对滑动距离为

dx＝2R·α(t)·f_s(t)·dt(12)

在实际使用中，摆动角度α(t)，摆动频率f_s(t)可能随着任务的不同而产生变化，故而它们是工作时间的函数。故

$du=\frac{2{\mathrm{Rk}}_s}{{\mathrm{\σ}}_s}{p}_0\left(u\right)\·\mathrm{\α}\left(t\right)\·f_s\left(t\right)\·dt---\left(13\right)$

时间磨损率w为单位时间内的磨损深度，即

$w=\frac{du}{dt}=\frac{2{\mathrm{Rk}}_s}{{\mathrm{\σ}}_s}{p}_0\left(u\right)\·\mathrm{\α}\left(t\right)\·f_s\left(t\right)---\left(14\right)$

从磨损率表达式(14)可知：在其他因素不变的情况下，磨损率与摆动角度和摆动频率成正比；如果磨损常数、摆动频率和摆动角度都为定值，且摩擦副表面状态不变，磨损率并不会保持常数，而是随磨损深度的增加而增加。公式(14)也反映出磨损率与磨损常数理论上成正比，在相同的工况下，如果磨损过程中磨损常数发生变化，则磨损率也成比例发生变化。

由公式(14)得到的时间磨损率反映的并不是自润滑衬垫材料的磨损率，而是自润滑关节轴承的磨损率，其与关节轴承的结构参数R有关。而磨损常数只与材料和摩擦副表面接触特性有关，所以关节轴承的磨损常数k_s是比磨损率w更基本的磨损特征量。同种衬垫的自润滑关节轴承，磨损常数遵循相同的物理规律。

由于关节轴承的一个摆动周期时间短，磨损增量非常微小，在理论建模过程中认为单个摆动周期内的最大接触压力保持恒定。由于最大关节轴承接触压力是磨损量的函数，而磨损量为时间的函数，故而最大关节轴承接触压力也是时间的函数，那么从时间t₀到t_T内的累计磨损量为

$u=\frac{2{\mathrm{Rk}}_s}{{\mathrm{\σ}}_s}{\∫}_{t_0}^{t_T}{p}_0\left(u\right(t\left)\right)\·\mathrm{\α}\left(t\right)\·f_s\left(t\right)\·dt---\left(15\right)$

如果关节轴承的摆动角度和运动频率都是恒定值，那么在时间t₀到t_T内的累计磨损量为

$u=2{\mathrm{R\αf}}_s\frac{k_s}{{\mathrm{\σ}}_s}{\∫}_{t_0}^{t_T}{p}_0\left(u\right(t\left)\right)\·dt---\left(16\right)$

式(16)成立的前提条件是磨损常数k_s在时间段t₀到t_T内保持不变。如果摩擦副的磨损常数不是恒定不变的量，则时间t₀到t_T内的累计磨损量为

$u=2{\mathrm{R\αf}}_s\frac{1}{{\mathrm{\σ}}_s}{\∫}_{t_0}^{t_T}{k}_s\left(t\right)\·p_0\left(u\right(t\left)\right)\·dt---\left(17\right)$

磨损常数是摩擦副的本质属性，与材料特性和摩擦副的接触状况有关。如果磨损常数不恒定，就其本质而言，磨损常数将随着磨损深度而变化。即磨损常数是磨损量的函数。所以，上式变为

$u=2{\mathrm{R\αf}}_s\frac{1}{{\mathrm{\σ}}_s}{\∫}_{t_0}^{t_T}{k}_s\left(u\right(t\left)\right)\·p_0\left(u\right(t\left)\right)\·dt---\left(18\right)$

当u＝u_m时，认为寿命终结，对应的工作时间即为关节轴承的寿命T。

从公式(8)和(9)可知，由于关节轴承协调接触模型确定的关节轴承最大接触压力并不是间隙s的显示表达式，故式(18)的积分求解需要利用数值方法迭代求解。迭代的方程为

$u_{i+1}-u_i=\frac{2R}{{\mathrm{\σ}}_s}{k}_s\left(u_i\right)\·p_0\left(u_i\right)\·{\mathrm{\α}}_i\·f_{s,i}\·(t_{i+1}-t_i),i=0,1,2...\left(19\right)$

当摆动角度和运动频率都是恒定值时，迭代求解的方程为

$u_{i+1}-u_i=2{\mathrm{R\αf}}_s\frac{1}{{\mathrm{\σ}}_s}{k}_s\left(u_i\right)\·p_0\left(u_i\right)(t_{i+1}-t_i),i=0,1,2...\left(20\right)$

当磨损常数在区间u∈(u_x,u_y]内保持不变，那么迭代方程为

$u_{i+1}-u_i=2{\mathrm{R\αf}}_s\frac{k_s}{{\mathrm{\σ}}_s}\·p_0\left(u_i\right)(t_{i+1}-t_i),u_i\∈(u_x,u_y\],i=0,1,2...$

ΔR_i＝(u₀+u_i)/2，t₀＝0(21)

u₀为关节轴承初始间隙，通过测量得到。

步骤4：选取关节轴承样本，在关节轴承试验台架上开展关节轴承磨损寿命试验，获取关节轴承结构和动态磨损数据，得到间隙随磨损时间的关系并绘制动态磨损曲线。

该步骤实际就是数据采样的过程，其步骤方法是工程实践中众所周知的方法。具体的操作方法是：选取关节轴承样本，在关节轴承试验台架上开展关节轴承磨损寿命试验，试验在恒温条件下开展，温度波动范围为5％，安装温度传感器进行温度跟踪。试验过程中应用精密位移传感器采集关节轴承内外圈间隙变化的动态磨损数据，采用等时间间隔采样，间隔为1-3分钟。台架的参数根据相关厂家的产品确定。试验可根据特定需求采用恒定载荷试验和变载荷试验。当采用变载荷试验时，需要明确施加变化载荷的时间点，应用步骤2中的方法计算不同载荷下的接触压力。全寿命试验的截止条件为最大累计磨损深度小于关节轴承允许的极限磨损深度。记录间隙变化随磨损时间的关系并画图，即得到动态磨损曲线。

步骤5：动态磨损过程拐点辨识。

在实际测试中，由于系统误差和随机误差的存在，摩擦副的动态磨损曲线不可能光滑，且可能存在较大波动，此时准确辨识离散曲线的临界“拐点”成为磨损分析的难点之一。本发明应用多项式逼近法寻找动态磨损曲线的拐点，利用n次高阶多项式函数y＝f(t)＝p₁tⁿ+p₂t^n-1+…p_nt+p_n+1逼近动态磨损曲线，其中p₁…p_n+1多项式系数。若动态磨损曲线波动过大，可采用分段拟合的方法。磨损曲线拐点的判定准则是：动态磨损曲线上某点的曲率K_ρ在转换区域中取得极大值。

$K_{\mathrm{\ρ}}=\frac{\left|f^{\′\′}\left(t\right)\right|}{{\[1+f^{\′}{\left(t\right)}^2\]}^{3/2}}---\left(22\right)$

其中，f′(t)为拟合多项式的一阶导数，也是该磨损曲线的动态磨损率；f″(t)为多项式的二阶导数值。

步骤6：磨损阶段磨损常数辨识。

当其他条件不变时，自润滑关节轴承寿命周期内磨损率出现阶段性波动是由磨损常数变化和接触压力增大耦合作用的结果。由于决定磨损常数的本质原因是材料特性和摩擦副接触特性，所以可认为每个磨损阶段内的磨损常数保持恒定。

在磨损失效物理方法中，关节轴承自润滑衬垫的磨损常数k_s需要通过磨损试验确定其变化规律。通过关节轴承的磨损曲线可计算出动态磨损过程的磨合磨损常数、稳定磨损常数和急剧磨损常数。由于试验测量的磨损过程存在随机误差，所以磨损曲线存在较大波动。基于拐点辨识的结果，磨损常数的求解方法为分阶段逼近，求得使试验动态磨损曲线与失效物理模型得到理论曲线之间的误差平方和SSE取最小的k_s。

$SSE=\underset{i=1}{\overset{N}{\Σ}}{(u_{PoF}\[i\]-u_{Test}\[i\])}^2---\left(23\right)$

其中N为动态磨损数字采样点数，u_PoF[i]为失效物理模型计算结果，u_Test[i]为试验测试结果。

这里，本发明磨损常数辨识采用的最小误差平方和方法也是数学上的一种通用方法。

步骤7：非线性磨损过程失效物理描述。

在步骤6中认为每个磨损阶段内的磨损常数保持恒定，在此基础上，分别定义磨合磨损常数、稳定磨损常数和急剧磨损常数为k_s,I，k_s,II,k_s,III。则动态磨损过程用失效物理方程描述为

$\begin{array}{c}u=2{\mathrm{R\αf}}_s\frac{1}{{\mathrm{\σ}}_s}\{k_{s,I}\·{\∫}_0^{t_I}{p}_0\left(u\right(t\left)\right)\·dt+k_{s,II}.{\∫}_{t_I}^{t_{II}}{p}_0\left(u\right(t\left)\right)\·dt+\\ k_{s,III}\·{\∫}_{t_{II}}^T{p}_0\left(u\left(t\right)\right)\·dt\}\end{array}---\left(24\right)$

其中，t_I,为磨合磨损阶段转化为稳定磨损阶段的时间点，等效为磨损量达到稳定磨损阶段的临界工作时间；t_II为稳定磨损阶段进入急剧磨损阶段的时间点，等效为磨损量达到急剧磨损阶段的临界工作时间。T为关节轴承寿命。

由于关节轴承产品的一致性不同，每个产品的磨损过程都呈现一定的差异，对于某些关节轴承，急剧磨损阶段不明显或不会出现，此时公式(24)中无第三项。

与现有数据驱动的关节轴承磨损过程分析方法相比，本方法具有以下显著优点：

(1)本发明提出了融合关节轴承的磨损特性参数、结构参数、材料参数和载荷工况参数的失效物理模型，参数之间的物理关系明确，能够用于关节轴承的结构优化设计和磨损分析中。相关研究填补了相关研究领域的空白。

(2)本发明建立的失效物理模型能十分准确地描述自润滑关节轴承起始于微小间隙的磨损退化过程，得到磨损的连续动态变化，提高了磨损分析的准确性。

(3)本发明提出的磨损阶段拐点辨识算法考虑了磨损过程的特征，符合工程实际，使拐点的辨识大大减少了主管因素的影响，更加客观准确。

本发明提出关节轴承滑动磨损的失效物理模型，应用基于失效物理模型的动态磨损过程分段分析法描述自润滑关节轴承的非线性磨损退化过程。本方法融合了关节轴承磨损特性、结构参数、材料参数和承载工况参数，能够非常准确地描述其磨损过程，并为关节轴承的设计分析提供支撑。

附图说明

图1典型磨损过程示意图

图2具体实施例的流程图

图3自润滑关节轴承动态磨损曲线的拐点辨识

图3a和3b分别为两个试验样本(M11和M14)的动态磨损曲线的高阶多项式逼近和曲率曲线，其中波动的折线为试验测得的动态磨损曲线，平滑实线为高阶多项式逼近的动态磨损曲线，虚线为高阶多项式的曲率。由图可知，应用步骤5中拐点辨识方法能准确辨识出磨损曲线存在的两个拐点，这两个拐点将磨损过程分为三个阶段。

图4自润滑关节轴承磨损过程的失效物理描述

图4a和4b给出了利用失效物理模型公式(24)，采用步骤1至步骤6的结果，分析载荷24kN的作用下五个样本的动态磨损过程，其中，图4a为样本M11、M12、M13，图4b为样本M14、M15。

具体实施方式

下面结合自润滑关节轴承恒定载荷磨损试验数据，按照图2所示流程说明本发明一种关节轴承磨损失效物理建模与分析方法的具体实施方式。

步骤1：建立自润滑关节轴承失效物理模型。

关节轴承磨损失效物理模型的建立过程如发明内容中步骤1-步骤7所示，最终得到的失效物理方程表达式为公式(24)。

步骤2：获取关节轴承结构和动态磨损数据。

本发明专利的实施方式中选择的关节轴承案例为两面带密封圈的PTFE衬垫单缝外环自润滑向心关节轴承，型号为GE20ET-2RS。此型号关节轴承的相关参数如表1所示。

表1试验轴承GE20ET-2RS参数表

试验采用的恒定应力载荷为24kN，磨损失效阈值为0.3mm，试验参数与试验结果如表2所示。5个样本的试验编号分别为M11-M15，试验样本摆动频率为32.7次/min，摆动振幅为在平衡位置±18.3°

表2关节轴承加速试验参数与试验结果

在关节轴承加速磨损试验中，通过位移传感器实时采样测量并记录关节轴承的间隙变化。由于在载荷的作用下试验机械系统会产生变形，且不同载荷下的变形量不同，所以要对动态磨损曲线进行归零处理，消除由载荷产生的变形带来的系统误差。消除此项误差的方法是将每个样本磨损试验初始10min内所测得的位移量的平均值作为磨损量的起始点。

本发明开展的关节轴承磨损试验的衬垫材料为PTFE芳纶织物，通过粘接工艺成型后的衬垫力学性能得到较大强化。本发明采用425MPa作为PTFE芳纶织物衬垫的压缩屈服强度估计值进行磨损计算。

步骤3：动态磨损过程拐点辨识。

采用拐点辨识算法辨识磨损曲线的拐点。图3为样本M11和M14的动态磨损曲线的高阶多项式逼近和曲率曲线，其中波动的折线为试验测得的动态磨损曲线，平滑实线为高阶多项式逼近的动态磨损曲线，虚线为高阶多项式的曲率。由图可知，应用公式(22)所确定的方法能准确辨识出磨损曲线存在的两个拐点，这两个拐点将磨损过程分为三个阶段：磨合磨损、稳定磨损和急剧磨损。需要指明的是，拐点的判定还需结合动态磨损曲线的形状而确定，在某些时刻，虽然曲线的曲率在某时刻局部取得极大值，但该点并不是拐点。只有当曲率在某时刻取极大值，且动态磨损曲线上同时出现相对明显的弯曲变化时，此点才是磨损阶段转换的临界拐点。

应用多项式拟合的方法可从图中准确地判定曲率极大值点所在的时刻，如图中椭圆圈所示，该时刻也是对应的动态磨损曲线的拐点所在位置，如箭头所示。相对其它拐点辨识方法，多项式拟合的方法不受磨损采样分辨率的限制，判定的拐点客观准确。

采用拐点辨识得到的各样本的拐点坐标如表3所示。

步骤4：磨损阶段磨损常数辨识。

采用磨损常数辨识算法，得到动态磨损曲线的磨合磨损常数、稳定磨损常数和急剧磨损常数如表3所示。

表3样本动态磨损过程拐点辨识与磨损常数计算结果

表中，“×”代表该样本未出现拐点2，未出现拐点2并不是由于试验时间过短所致的，而是样本的特殊性导致急剧磨损阶段与稳定磨损阶段重合，此时，认为急剧常数与稳定磨损常数相等。

步骤5：非线性磨损过程失效物理描述

利用失效物理模型公式(24)，采用前面步骤所提供的数据，分析24kN下关节轴承GE20ET-2RS的动态磨损过程如图4所示。

Claims (1)

1.一种关节轴承磨损失效物理建模与分析方法，其特征在于，包括以下几个步骤：

步骤1、确定材料磨损模型

将强度作为衡量自润滑衬垫层耐磨损性能的量，所以关节轴承滑动磨损的磨损体积公式为

$V=k_s\frac{F_N}{{\mathrm{\σ}}_s}x---\left(1\right)$

其中，V为磨损量，F_N为法向载荷，x为滑动距离，σ_s为接触摩擦副中较软材料的受压屈服极限，k_s称为磨损常数；

实际使用中将关节轴承内外球面间的最大允许间隙作为其磨损失效阈值；关节轴承的结构间隙s定义为内外球面间的间隙；初始间隙u₀为产品出厂时的内外球面间隙；结构间隙为初始间隙与磨损深度u的叠加，s＝u₀+u；在关节轴承接触压力计算公式(8)和(9)中，R₁为内球面半径，R₂为外球面半径；在关节轴承磨损分析中，R₁和R₂为相对值，在磨损的过程中认为R₂＝R＝d_k/2为恒定值，而R₁随着磨损深度而变化；其中，d_k为关节轴承协调接触表面的直径；即有

R₁＝R-(u₀+u)/2(2)

ΔR＝s/2＝(u₀+u)/2(3)

由于磨损量与载荷成正比，所以最大磨损量取决于接触载荷最大的点；关节轴承中的接触应力最大点p₀处于接触区域的中心点O，接触压力随着磨损深度变化而变化；在接触区域的中心，取一个微元区，其面积为A_a；设微元区为平面，且微元区上的接触压力均匀，都是p₀，所以，

F_N＝p₀A_a,V＝uA_a(4)

由式(1)可得用磨损深度表示的Archard磨损模型为

$u=k_s\frac{p_0\left(u\right)}{{\mathrm{\σ}}_s}x---\left(5\right)$

其中，p₀(u)为当磨损深度为u时的最大关节轴承接触压力；

步骤2、确定关节轴承接触压力计算方法

关节轴承接触压力计算方法如下，

(1)当0＜a＜h时，关节轴承中的接触压力分布采用完整球面协调接触的模型计算，即

$\left\{\begin{array}{cc}{p}_0=(n+1)\frac{F}{{\mathrm{\πa}}^2}& \left(a\right)\\ a={\[\frac{4{\mathrm{BR}}_1{R}_{2}F}{{\mathrm{\π}}^2{E}^{*}(g\mathrm{\Δ}R+c)}(n+1/2)(n+1)\]}^{1/3}& \left(b\right)\end{array}\right.---\left(6\right)$

(2)当h＜a≤R₂时，其求解方法为

$\left\{\begin{array}{cc}Q={\∫}_h^{a}r\·{(1-r^2/a^2)}^n\·arccos\frac{h}{r}\·dr& \left(c\right)\\ F_0=\frac{4(n+1)\·Q\·F}{{\mathrm{\πa}}^2-4(n+1)\·Q},F_t=F+F_0& \left(d\right)\\ p_0=(n+1)\frac{F_t}{{\mathrm{\πa}}^2}& \left(e\right)\\ a={\[\frac{4{\mathrm{BR}}_1{R}_2{F}_t}{{\mathrm{\π}}^2{E}^{*}(g\mathrm{\Δ}R+c)}(n+1/2)(n+1)\]}^{1/3}& \left(f\right)\end{array}\right.---\left(7\right)$

对于以上两式，都有

$\left\{\begin{array}{cc}n=0.5-0.24\exp \[-15.08(1-a/R_2)\]& \left(g\right)\\ B=\frac{\sqrt{\mathrm{\π}}\mathrm{\Γ}(n+1)}{2\mathrm{\Γ}(3/2+n)},c=\frac{3.8304B_{t}F}{{\mathrm{\π}}^2{E}^{*}{R}_2}& \left(h\right)\\ g\left(a\right)=2/\mathrm{\π}+{(a/R_2)}^2& \left(j\right)\end{array}\right.---\left(8\right)$

其中：r为球面上的点到接触区域中心的投影距离；a为接触区域的投影半径；h为关节轴承半宽；p₀为最大关节轴承接触压力；F为关节轴承承受的外部施加载荷；n为压力分布指数；F₀为等效附加载荷；Γ为gamma函数；B，c为中间参数，用公式(8h)定义；E^*为等效弹性模量，设其中，E₁,E₂分别为关节轴承内圈与外圈材料的弹性模量，μ₁,μ₂为相应的泊松比；

步骤3、推导动态磨损量表达式与数值求解方法

一个微小的相对滑动距离dx内，其最大关节轴承接触压力可认为是恒定的，那么

$du=k_s\frac{p_0}{{\mathrm{\σ}}_s}\·dx---\left(9\right)$

其中du为滑动距离dx内磨损量的增量；

对于关节轴承，其工作过程中的运动方式一般为摆动，设摆动角度为±α，摆动频率为f_s；所以，在时间dt内，关节轴承内外环之间的相对滑动距离为

dx＝2R·α(t)·f_s(t)·dt(10)

在实际使用中，摆动角度α(t)，摆动频率f_s(t)可能随着任务的不同而产生变化，故而其是工作时间的函数；故

$du=\frac{2{\mathrm{Rk}}_s}{{\mathrm{\σ}}_s}{p}_0\left(u\right)\·\mathrm{\α}\left(t\right)\·f_s\left(t\right)\·dt---\left(11\right)$

时间磨损率w为单位时间内的磨损深度，即

$w=\frac{du}{dt}=\frac{2{\mathrm{Rk}}_s}{{\mathrm{\σ}}_s}{p}_0\left(u\right)\·\mathrm{\α}\left(t\right)\·f_s\left(t\right)---\left(12\right)$

由于关节轴承的一个摆动周期时间短，磨损增量非常微小，在理论建模过程中认为单个摆动周期内的最大接触压力保持恒定；由于最大关节轴承接触压力是磨损量的函数，而磨损量为时间的函数，故而最大关节轴承接触压力也是时间的函数，那么从时间t₀到t_T内的累计磨损量为

$u=\frac{2{\mathrm{Rk}}_s}{{\mathrm{\σ}}_s}{\∫}_{t_0}^{t_T}{p}_0\left(u\right(t\left)\right)\·\mathrm{\α}\left(t\right)\·f_s\left(t\right)\·dt---\left(13\right)$

如果关节轴承的摆动角度和运动频率都是恒定值，那么在时间t₀到t_T内的累计磨损量为

$u=2{\mathrm{R\αf}}_s\frac{k_s}{{\mathrm{\σ}}_s}{\∫}_{t_0}^{t_T}{p}_0\left(u\right(t\left)\right)\·dt---\left(14\right)$

式(14)成立的前提条件是磨损常数k_s在时间段t₀到t_T内保持不变；如果摩擦副的磨损常数不是恒定不变的量，则时间t₀到t_T内的累计磨损量为

$u=2{\mathrm{R\αf}}_s\frac{1}{{\mathrm{\σ}}_s}{\∫}_{t_0}^{t_T}{k}_s\left(u\right(t\left)\right)\·p_0\left(u\right(t\left)\right)\·dt---\left(15\right)$

当u＝u_m时，认为寿命终结，对应的工作时间即为关节轴承的寿命T；

从公式(6)和(7)可知，由于关节轴承协调接触模型确定的关节轴承最大接触压力并不是间隙s的显示表达式，故式(15)的积分求解需要利用数值方法迭代求解；迭代的方程为

$u_{i+1}-u_i=\frac{2R}{{\mathrm{\σ}}_s}{k}_s\left(u_i\right)\·p_0\left(u_i\right)\·{\mathrm{\α}}_i\·f_{s,i}\·(t_{i+1}-t_i),i=0,1,2...---\left(16\right)$

当摆动角度和运动频率都是恒定值时，迭代求解的方程为

$u_{i+1}-u_i=2{\mathrm{R\αf}}_s\frac{1}{{\mathrm{\σ}}_s}{k}_s\left(u_i\right)\·p_0\left(u_i\right)(t_{i+1}-t_i),i=0,1,2...---\left(17\right)$

当磨损常数在区间u∈(u_x,u_y]内保持不变，那么迭代方程为

$u_{i+1}-u_i=2{\mathrm{R\αf}}_s\frac{k_s}{{\mathrm{\σ}}_s}\·p_0\left(u_i\right)(t_{i+1}-t_i),u_i\∈(u_x,u_y\],i=0,1,2...$

ΔR_i＝(u₀+u_i)/2,t₀＝0(18)

u₀为关节轴承初始间隙；

步骤4、选取关节轴承样本，在关节轴承试验台架上开展关节轴承磨损寿命试验，获取关节轴承结构和动态磨损数据，得到间隙随磨损时间的关系并绘制动态磨损曲线；

步骤5，非线性磨损过程的拐点辨识；

应用多项式逼近法寻找动态磨损曲线的拐点，利用n次高阶多项式函数逼近动态磨损曲线；若动态磨损曲线波动过大，可采用分段拟合的方法；磨损曲线拐点的判定准则是：动态磨损曲线上某点的曲率K_ρ在转换区域中取得极大值；

$K_{\mathrm{\ρ}}=\frac{\left|f^{\′\′}\left(t\right)\right|}{{\[1+f^{\′}{\left(t\right)}^2\]}^{3/2}}---\left(19\right)$

其中，f′(t)为拟合多项式的一阶导数，也是该磨损曲线的动态磨损率；f″(t)为多项式的二阶导数值；

步骤6，磨损阶段磨损常数辨识；

基于非线性磨损过程的拐点辨识的结果，磨损常数的求解方法为分阶段逼近，求得使动态磨损曲线与失效物理模型得到理论曲线之间的误差平方和SSE取最小的k_s；

$SSE=\underset{i=1}{\overset{N}{\Σ}}{(u_{PoF}\[i\]-u_{Test}\[i\])}^2---\left(20\right)$

其中N为动态磨损数字采样点数，u_PoF[i]为失效物理模型计算结果，u_Test[i]为试验测试结果；

步骤7，动态磨损过程关节轴承磨损失效物理模型；

认为每个磨损阶段内的磨损常数保持恒定，分别定义磨合磨损常数、稳定磨损常数和急剧磨损常数为k_s,I，k_s,II,k_s,III；则动态磨损过程用失效物理模型为

$u=2{\mathrm{R\αf}}_s\frac{1}{{\mathrm{\σ}}_s}\{k_{s,I}\·{\∫}_0^{t_I}{p}_0\left(u\right(t\left)\right)\·dt+k_{s,II}\·{\∫}_{t_I}^{t_{II}}{p}_0\left(u\right(t\left)\right)\·dt+k_{s,III}\·{\∫}_{t_{II}}^T{p}_0\left(u\right(t\left)\right)\·dt\}---\left(21\right)$

其中，t_I,为磨合磨损阶段转化为稳定磨损阶段的时间点，等效为磨损量达到稳定磨损阶段的临界工作时间；t_II为稳定磨损阶段进入急剧磨损阶段的时间点，等效为磨损量达到急剧磨损阶段的临界工作时间；p₀(u(t))采用公式(6)和(7)计算；T为关节轴承寿命。