齿轮的磨损会造成接触环境恶化，加大齿间动态载荷，而载荷的增大又会加重齿轮磨损，磨损会造成更严重的轮齿失效形式.揭示磨损与齿轮参数的关系可实现主动抑制磨损，降低齿轮的振动冲击，延长齿轮使用寿命.

目前，各类齿轮系统磨损机理的研究还不尽完善，加之齿轮的啮合及运动学参数复杂，很难建立起磨损与齿轮参数间的直接关系.以往的研究中给定齿轮参数后，都是通过确定啮合载荷和齿面滑动距离(速度)后，应用磨损机理来确定磨损量.在对齿轮磨损的数值计算中，国内外学者普遍采用Archard磨损模型^[1-5]，只是在确定齿轮的啮合载荷和齿面滑动距离时方法不尽相同，Anddersson^[1]通过齿轮几何学理论得到了啮合过程中滑动距离与啮合位置的关系. Flodin等^[2]通过单点观测法确定了滑动速度和啮合压力. Wu等^[3]在考虑了齿轮动力学和弹性流体动力润滑的影响下，推导了齿轮磨损计算公式.王晓笋等^[4]和Kuang等^[5]研究了齿面累积磨损量对轮齿啮合刚度的影响并进行了动力学分析.上述模型的齿间载荷只是进行了简单处理，并未考虑磨损带来的间隙对载荷的影响.

磨损是一个逐渐累积的渐变过程，达到极限值的时间历程很长，为了保证足够精度，数值仿真的方法给定的模拟步长又不能太大，这两个因素都会导致进行一次磨损量计算将耗费大量的时间，为了建立磨损与齿轮参数的关系，需要进行大样本模拟，若参数样本较多时，这种数值仿真算法将不再适用.

为了解决上述问题，可以使用近似模型来代替磨损数值模型.常用的近似逼近方法有响应面方法^[6]、自适应神经网络方法^[7]、Kriging方法^[8]、径向基方法^[9]等.响应面模型构造的函数形式固定，不适合解决复杂问题，自适应神经网络的稳健性低，同样精确度时的计算效率不如Kriging模型. Laurenceau等^[10]将Kriging模型用在飞机动力学设计中，并与响应面法进行了对比. Lucifredi等^[11]将动态Kriging模型应用于水电力系统的预测维护中，并对比了它与神经网络方法的计算精度与效率.本文采用Kriging方法对齿面磨损进行预测，用Kriging模型代替磨损仿真模型.

1 齿面磨损模型

齿轮在高接触压力的作用下，相互啮合的轮齿表面发生相对滚动和滑动.虽然大部分齿轮采用了油脂润滑，但是润滑状态还是处于边界润滑或混合润滑，啮合表面不能被润滑剂完全隔开，摩擦表面会产生金属间的直接摩擦.此时，相互啮合的轮齿表面在低速重载作用下，会形成材料的转移以及表层金属剥落，造成齿轮的磨损.

Flodin等^[2]基于广义的Archard磨损方程提出了直齿轮轻微磨损预测模型.假定在极短的时间内，齿廓上的点i的表面压力p_i和滑移速度v_i保持不变.因此，点i的磨损量可表示为

$ {{h}_{i, n}}={{h}_{i, n-1}}+\Delta tkN\sum\limits_{j=1}^{s}{{{p}_{i, j}}{{v}_{i, j}}}. $

式中：h_{i, n}为啮合点i处n次磨损后的磨损深度，h_{i, n-1}为啮合点i处n-1次磨损后的磨损深度，n为当前磨损次数，Δt为时间步长，k为磨损系数，N为每次磨损运行间隔的转数，在此期间齿轮轮廓将不进行更新，s为赫兹接触半径内的点. p_{i, j}为在j点接触时i点的表面压力，v_{i, j}为点i的滑移速度.

1.1 表面压力

基于Winkler弹性模型可得到在啮合点j接触时处于接触半径内齿廓点i的表面压力p_{i, j}:

$ {{p}_{i, j}}({{x}_{i, j}})=\frac{(0.59{{E}^{*}})}{Ra}({{a}^{2}}-x_{i, j}^{2}). $

(1)

式中：j为啮合点，a为接触半径，x_{i, j}为接触半径内的一点i到啮合点j的距离，E^*为等效弹性模量，R为两个圆柱的综合曲率半径.式(1)中的间接变量计算公式为

$ a={{\left( \frac{4WR}{\mathsf{ π} {{E}^{*}}}~ \right)}^{1/2}}. $

式中：W为W₁或W₂，W₁为第一对轮齿啮合副分担的法向力，W₂为第二对轮齿啮合副分担的法向力，由式(2)和式(3)确定，当第1对齿轮先接触时，有

$ \left\{ \begin{align} &{{W}_{1}}=\frac{{{k}_{1}}\left( t \right)\left( F+\left| \mathit{\Delta } \right|\cdot {{k}_{2}}\left( t \right) \right)\text{ }}{{{k}_{1}}\left( t \right)+{{k}_{2}}\left( t \right)}, \\ &{{W}_{2}}=\frac{{{k}_{2}}\left( t \right)\left( F-\left| \mathit{\Delta } \right|\cdot {{k}_{1}}\left( t \right) \right)}{{{k}_{1}}\left( t \right)+{{k}_{2}}\left( t \right)}. \\ \end{align} \right. $

(2)

当第2对齿轮先接触时，

$ ~\left\{ \begin{align} &{{W}_{1}}=\frac{{{k}_{1}}\left( t \right)\left( F-\left| \mathit{\Delta } \right|\cdot {{k}_{2}}\left( t \right) \right)}{{{k}_{1}}\left( t \right)+{{k}_{2}}\left( t \right)}\text{, } \\ &{{W}_{2}}=\frac{{{k}_{2}}\left( t \right)\left( F+\left| \mathit{\Delta } \right|\cdot {{k}_{1}}\left( t \right) \right)}{{{k}_{1}}\left( t \right)+{{k}_{2}}\left( t \right)}. \\ \end{align} \right. $

(3)

式中：F为单位宽度的法向总载荷，k₁(t)和k₂(t)为齿对1和齿对2单位齿宽的啮合刚度. Δ为齿轮的几何侧隙.

当啮合点在B₁C段时，

$ \left| \mathit{\Delta } \right|=\left| {{h}_{1}}\left( y \right)+{{h}_{2}}\left( y \right)-\left( {{h}_{1}}\left( y+{{p}_{\text{b}}} \right)+{{h}_{2}}\left( y+pb \right) \right) \right|; $

当啮合点在DB₂段时，

$ \left| \mathit{\Delta } \right|=\left| {{h}_{1}}\left( y \right)+{{h}_{2}}\left( y \right)-\left( {{h}_{1}}\left( y-{{p}_{\text{b}}} \right)+{{h}_{2}}\left( y-{{p}_{\text{b}}} \right) \right) \right|. $

如图 1所示，y为接触点i和节点i’沿着啮合线的距离，h(y)为磨损深度.

1.2 滑动速度

应用齿轮啮合理论可确定点i的滑动速度：

$ \begin{align} &\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {{v}_{i}}=\left( {{\omega }_{1}}+{{\omega }_{2}} \right)\times {{y}_{i}}, \\ &{{y}_{i}}=i{{i}^{\prime }}={{r}_{\text{b1}}}\left( \text{tan}\ {{\beta }_{1i}}-\text{tan}\ \phi \right)={{r}_{\text{b2}}}(\text{tan}\ \phi -\text{tan}\ {{\beta }_{2i}}). \\ \end{align} $

式中：r_b1r_b2为主从动轮的基圆半径，β_1iβ_2i分别为啮合点i在齿轮1和齿轮2上对应的压力角，ϕ为啮合角.

2 Kriging预测模型

上述磨损模型的工作参数、尺寸参数均是确定的，属于确定性模型.当N=1 000转，磨损次数为450时，计算时间大约为1 h，若给定另一组参数，程序需要重新计算并消耗同样的时间.在大量样本面前，直接计算的方法将不再适用，采用Kriging模型可以建立响应与模型参数的关系，当给定样本后可以快速精确地预测出磨损量(响应值).

同时，磨损是一个逐渐累积的过程，当磨损量达到极限时，齿轮失效，由齿面磨损的计算公式可以看出，在计算过程中假定了在N转内齿廓不进行更新，如果N的取值太大会造成计算精度低，因此在允许的范围内使N的取值尽可能小，但会造成计算时间长.在Kriging模型中将时间参数考虑进去，可以快速精确地预测磨损的变化趋势.

Kriging模型将模型参数看成随机变量，将磨损量看成随机响应，并将随机响应模型描述为两部分，参数化的线性回归模型和非参数化的随机过程模型.为了满足模拟过程的无偏性，可得随机响应的估计值为

$ \mathit{\boldsymbol{\hat h}}\left( \mathit{\boldsymbol{x}} \right) = \mathit{\boldsymbol{f}}{\left( \mathit{\boldsymbol{x}} \right)^{\rm{T}}}{\mathit{\boldsymbol{\beta }}^ * } + \mathit{\boldsymbol{r}}{\left( \mathit{\boldsymbol{x}} \right)^{\rm{T}}}{\mathit{\boldsymbol{R}}^{ - 1}}(\mathit{\boldsymbol{H}} - \mathit{\boldsymbol{F}}{\mathit{\boldsymbol{\beta }} ^ * }). $

(4)

式中：x={x₁, x₂, …, x_n+1}为待测点，n+1为变量数；f(x)^T={f₁(x), f₂(x), …, f_p(x)}为向量x的多项式，p为多项式数；β^*={β₁^*, β₂^*, …, β_p^*}^T为模拟的广义最小二乘估计；假设有m个样本，分别为x_s1, x_s2, …, x_sm；r(x)={R(θ_k, x, x_s1), R(θ_k, x, x_s2), …, R(θ_k, x, x_sm)}^T为待测点x和样本x_s1, x_s2, …, x_sm之间的相关向量；R_m×m为样本间的相关矩阵；H为响应样本值；F为样本的多项式；β^*为模拟的最小二乘估计，

$ {\mathit{\boldsymbol{\beta }}^ * } = {\left( {{\mathit{\boldsymbol{F}}^{\rm{T}}}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{F}}} \right)^{ - 1}}{\mathit{\boldsymbol{F}}^{\rm{T}}}{\mathit{\boldsymbol{R}}^{ - 1}}\mathit{\boldsymbol{Y}}. $

式中：R(θ_k, x_i, x_j)是样本点中任意两个样本点x_i和x_j的空间相关方程，i, j∈s₁, …, s_m由用户自定义，其中计算效果最好，被广泛采用的是高斯相关方程.

$ R({\theta _k}, {\mathit{\boldsymbol{x}}_i}, {\mathit{\boldsymbol{x}}_j}) = \prod\limits_k^{n + 1} {{\rm{exp}}( - {\theta _k}\left| {x_k^i - x_k^j} \right|{^2})} . $

由R(θ_k, x_i, x_j)可组集为相关矩阵

$ \mathit{\boldsymbol{R = }}{\left[{\begin{array}{*{20}{c}} {\sigma _1^2}&{R\left( {{\mathit{\boldsymbol{x}}_1}, {\mathit{\boldsymbol{x}}_2}} \right)}& \cdots &{R\left( {{\mathit{\boldsymbol{x}}_1}, {\mathit{\boldsymbol{x}}_m}} \right)}\\ {R\left( {{\mathit{\boldsymbol{x}}_2}, {\mathit{\boldsymbol{x}}_1}} \right)}&{\sigma _2^2}&{}&{}\\ \vdots &{}&{}&{}\\ {R\left( {{\mathit{\boldsymbol{x}}_m}, {\mathit{\boldsymbol{x}}_1}} \right)}&{}&{}&{\sigma _m^2} \end{array}} \right]_{\left( {m \times m} \right).}} $

应用最大似然估计法可得到相关方程中θ^T={θ₁, …, θ_k, …, θ_n+1}^T的最优解θ^*为

${\mathit{\boldsymbol{\theta }}^ * } = \mathop {{\rm{min}}}\limits_\theta \left\{ {{{\left| \mathit{\boldsymbol{R}} \right|}^{\frac{1}{m}}}{\sigma ^2}} \right\}. $

式中, m为样本数，|R|为相关矩阵R的行列式.

当样本给定后，可确定出β^*、R、H、F的具体表达式，当给出待测点x以及f(x)的回归模型，以及r(x)的相关方程，即可通过式(4)得到$\mathit{\boldsymbol{\hat h}}\left( \mathit{\boldsymbol{x}} \right)$.

3 算例

表 1为计算一对齿轮磨损需要的参数，材料为45号钢.齿轮尺寸参数、工作参数和磨损参数均参考自文献[2].虽然文献[2]只进行了理论分析，但是被大量学者认可，并在文献[12]中进行了实验验证.

定义磨损次数为150次，从第n-1次磨损到第n次磨损之间，齿轮转过了N=1 000转，并且1, 2, …, N转的任意一转中，齿轮的轮廓是恒定的，啮合点的压力和速度都相等.

由图 2和图 3可知，节点处的磨损最小，齿根处的磨损要大于齿顶的磨损，主动轮(小齿轮)的磨损大于从动轮(大齿轮)的磨损，这与参考文献[12-14]得到的结论一致.最大磨损出现在主动轮齿廓点233处和从动轮263处.如图 1所示，一对齿的啮合从B₁开始至B₂结束，最大磨损出现在啮入位置和啮出位置处.由图 2和图 3可知，由于考虑间隙的影响，在单双齿啮合转换处，磨损量出现先增大再减小的现象，不同于文献[2]中的一直减小的趋势.

若齿形角α_n、齿宽b、外加力矩T₁、小齿轮速度n₁以及磨损系数k为变量，时间t为1 000 min，采用拉丁超立方试验设计方法，分别得到变量和真实磨损组成的样本50、100、500组，建立3种不同样本量下的Kriging近似模型(用Kriging-50，Kriging-100，Kriging-500表示).回归模型中多项式p=0，相关模型为高斯方程. 3种Kriging模型的具体参数如表 2所示.

图 4、5为齿轮最大磨损量变化曲线，由图 4、5可知，磨损量随时间逐渐累积增加.在样本量为50、100和500时，建立Kriging模型的时间只需几秒钟，并且模型对真实磨损量的逼近精度很高，同时，真实磨损量只提供了t₀时间内的样本，Kriging模型可预测t＞t₀时的磨损量，并且在样本量为100组以上时，误差非常小.

为说明Kriging模型的优越性，本文同时采用人工神经网络模型(ANN)^[15]，建立了6-13-1形式的3层BP神经网络，隐含层的传递函数选用S型对数函数，输出层的传递函数选用purelin线性传递函数，训练函数选择trainlm函数，网络的训练误差达到10^-6时结果收敛. 图 6、7为不同样本量时ANN模型的模拟精度.

采用50组样本训练得到的ANN模型的精确程度不高，采用100、500组的模型精确程度都很高，并且对未来趋势的预测也比较准确.同样采用50、100和500组样本建立的Kriging近似模型，3个Kriging模型的精确程度都很高，对未来趋势的预测也比较准确.由此可看出，与ANN方法相比，Kriging方法可使用更少的样本得到较高的精度.

为了检验模型的通用性，另外取50组样本分别带入上述3种Kriging模型进行检验，将样本的真实值与Kriging模型的预测值进行对比，结果表明Kriging-50模型的拟合优度为0.764 5，Kriging-100模型的拟合优度为0.998 6，Kriging-500模型的拟合优度为0.999 9，可见当采用Kriging-100模型进行检验时，拟合优度可近似达到1，因此，只需要100组样本模拟得到的Kriging模型可近似代替磨损预测模型. 图 8为采用Kriging-100模型的磨损量预测值和真实值.

样本量为10⁵时，采用Kriging模型进行磨损量预测所需要的时间只有几秒钟，而用磨损仿真模型需要3.8年.可以看出，采用Kriging模型比直接采用磨损预测模型大大节约了时间和成本.

4 结论

本文考虑了磨损对齿廓表面几何形状的影响，推导了几何侧隙下的齿间载荷分配公式，并以此为基础，引入Archard磨损模型和Winkler弹性力学模型，建立了低速重载齿轮的磨损量计算模型.基于Kriging方法得到了近似逼近模型，建立了齿轮参数与磨损的直接关系.为减少齿轮振动冲击提供了途径.得到的结论如下：

1) 节点处的磨损最小，齿根处的磨损要大于齿顶的磨损，主动轮(小齿轮)的磨损大于从动轮(大齿轮)的磨损.

2) 节点到齿顶的磨损分布较均匀，而齿根到节点的磨损变化较剧烈.重点应该关注齿根附近的磨损.

3) 使用很少的样本建立的Kriging模型可以完全代替磨损数值模型，计算效率和拟合优度都非常高.可用于对未知磨损的预测.