计入冲击与受载齿轮副时变啮合刚度解析算法

Analytical Algorithm for Time Varying Mesh Stiffness of Gear Pairs Considering Impact and Loading

摘要

针对目前理论不能准确计算齿轮受载与冲击时刻刚度发生的变化，基于赫兹非线性接触理论，构建有限长线接触弹性理论解析模型，由三维接触斑大小与啮合过程中主曲率半径的变化，得到不同载荷作用下齿轮法向接触刚度；而后结合Weber能量法，提出考虑负载条件下时变啮合刚度修正算法，阐明受载变形与刚度硬化的映射关系；同时建立线外啮合冲击数学模型，根据冲击时刻齿轮副动态接触关系，探究啮入冲击引起的啮合刚度变化规律.研究结果表明，接触刚度曲线与主曲率半径变化规律一致，考虑啮入冲击时，齿轮副啮合刚度产生突变并呈现由大到小的变化规律，随着负载递增，法向接触与综合刚度非线性增大，冲击效应增强，且单齿啮合区域逐渐减小，冲击对刚度的影响越来越明显.

Abstract

Due to the inability of the current theory to accurately calculate the changes in the stiffness of the gear during loading and impact， the Hertz nonlinear contact theory is employed to construct an analytical model of length-limited line contact elasticity theory. The normal contact stiffness of the gear under different loads is obtained from the changes in the three-dimensional contact spot size and the principal curvature radius in the meshing process. Later， using the Weber energy method， a time-varying meshing stiffness correction algorithm considering the load condition is proposed to clarify the mapping relationship between load deformation and stiffness hardening. At the same time， the mathematical model of off-line meshing impact is established， and the change law of the meshing stiffness caused by the meshing impact is investigated according to the dynamic contact relationship of the gear pair at the impact time. The research findings demonstrate that the contact stiffness curve is consistent with the change law of the principal curvature radius. When the meshing impact is considered， the meshing stiffness of the gear pair changes abruptly from the large to the small. As the load increases， the contact and comprehensive meshing stiffness increase non-linearly， the single tooth meshing region gradually decreases， and the impact on the stiffness is more markedly.

Graphical abstract

关键词

刚度变化 / 赫兹接触 / Weber能量法 / 接触变形 / 冲击刚度

Key words

stiffness change / Hertz contact / Weber energy method / contact deformation / shock stiffness

齿轮传动系统具有传动精度高、啮合平稳等特点，被广泛应用于航空航天、医疗器械、工业机器等技术领域中.由于其使役工况复杂多变，运行负载变化范围明显，冲击碰撞问题尤为突出，这些会影响齿轮传动的安全与稳定.同时齿轮时变啮合刚度非线性特征变化显著，而能否明晰负载-冲击-刚度三者之间动态演变规律是突破现有齿轮研究的关键.因此，研究齿轮负载变化与线外啮合诱发的时变啮合刚度规律具有十分重要的意义.

齿轮系统啮合刚度计算受到国内外学者的广泛关注，目前主要研究方法分为有限元法^［1-3］与解析法.采用有限元法计算刚度时相当于模拟齿轮实际啮合情况，具有较高的计算精度，但因其求解结果基于建模水平、边界条件与工况模拟是否真实，且计算效率较低，故更多学者对高效率的解析法进行研究.刘文等^［4］基于势能法提出了一种考虑齿根圆与基圆不重合情形时齿轮时变啮合刚度修正算法.林腾蛟等^［5］运用切片法计算考虑温度情形下斜齿轮时变啮合刚度，同时采用不同工况参数对齿轮啮合刚度的影响进行分析.陈思宇等^［6］利用能量等效给出直齿锥齿轮啮合刚度求解的快速计算方法.Liu等^［7］采用解析几何法建立了内齿轮与外齿轮的齿尖碎裂失效模型，推导出故障齿轮的啮合刚度并进行了验证.Chen等^［8］依据直齿圆柱齿轮啮合刚度的计算，推导高阶定相齿轮副的啮合刚度并进行了优化.Sun等^［9］基于线性热膨胀理论，综合计算热膨胀下的综合啮合刚度.对于法向接触刚度的计算，一些研究者根据分形理论，建立法向刚度分型模型，考虑分形维数、表面粗糙度幅值与材料特性等情况下对法向接触刚度的影响^［10-12］.部分学者对啮合冲击的计算方法也进行了深入的研究^［13-14］.

综上所述，以往计算轮齿接触变形是将接触界面简化为沿啮合线方向随曲率半径而发生变化的光滑曲面，其载荷与接触变形为线性关系，并不能表征轮齿接触界面受载压平导致的接触变形非线性增大效应，同时未给出考虑负载与冲击的刚度解析计算公式，明晰三者之间的规律与变化.

因此，本文基于赫兹（Hertz）非线性接触理论，以点接触为基础构建线弹性接触解析计算模型，结合齿轮啮合过程中曲率半径与接触斑大小变化求解齿轮啮合过程中弹性变形量.提出了考虑负载条件下时变啮合刚度修正算法，计算不同载荷作用下齿轮法向接触与综合啮合刚度，并建立啮入冲击刚度耦合作用激励模型.推导出任意齿轮副不同工况下冲击诱发的刚度变化解析式，阐明了不同负载条件下冲击激励对刚度的作用影响.

1 赫兹非线性接触刚度解析算法

1.1 赫兹点接触与线接触理论

在实际工程运用中，赫兹接触理论通常运用于计算两构件产生的变形、最大接触应力与应力变化的分布规律.应用赫兹接触分析问题时需要满足一些条件：

1） 假设两构件属于完全弹性变形体，材料均匀分布；

2） 接触体保持各向同性，接触体表面光滑连续，并且接触区域远小于接触体尺寸；

3） 接触时压力分布与接触区域保持一致，同时在初始接触点附近两弹性体可视为弹性半空间.

根据这些条件处理一般光滑弹性体接触问题时，由于接触形状的不同，赫兹接触问题分为点接触与线接触两种，图1与图2表示点接触时接触区域与接触应力分布图.

根据赫兹理论，点接触时在法向力F_n的作用下，三维弹性体接触区域形状为椭圆，椭圆长短轴半径a、b分别为：

a=g{\left[\frac{\mathrm{3}\mathrm{\pi }{F}_{\mathrm{n}}\left({\gamma }_{\mathrm{1}}+{\gamma }_{\mathrm{2}}\right)}{\mathrm{4}(X+Y)}\right]}^{\mathrm{1}/\mathrm{3}}（1）

b=h{\left[\frac{\mathrm{3}\mathrm{\pi }{F}_{\mathrm{n}}\left({\gamma }_{\mathrm{1}}+{\gamma }_{\mathrm{2}}\right)}{\mathrm{4}(X+Y)}\right]}^{\mathrm{1}/\mathrm{3}}（2）

\mathrm{2}(Y+X)=\frac{\mathrm{1}}{R_{\mathrm{1}}^{\text{'}}}+\frac{\mathrm{1}}{R_{\mathrm{1}}^{″}}+\frac{\mathrm{1}}{R_{\mathrm{2}}^{\text{'}}}+\frac{\mathrm{1}}{R_{\mathrm{2}}^{″}}（3）

\mathrm{2}(Y-X)=\left[{\left(\frac{\mathrm{1}}{R_{\mathrm{1}}^{\text{'}}}-\frac{\mathrm{1}}{R_{\mathrm{1}}^{″}}\right)}^{\mathrm{2}}+{\left(\frac{\mathrm{1}}{R_{\mathrm{2}}^{\text{'}}}-\frac{\mathrm{1}}{R_{\mathrm{2}}^{″}}\right)}^{\mathrm{2}}+\right.

{\left.\mathrm{2}\left(\frac{\mathrm{1}}{R_{\mathrm{1}}^{\text{'}}}-\frac{\mathrm{1}}{R_{\mathrm{1}}^{″}}\right)\left(\frac{\mathrm{1}}{R_{\mathrm{2}}^{\text{'}}}-\frac{\mathrm{1}}{R_{\mathrm{2}}^{″}}\right)\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{2}\theta \right]}^{\frac{\mathrm{1}}{\mathrm{2}}}（4）

式中：R_{\mathrm{1}}^{\text{'}}与R_{\mathrm{1}}^{″}、R_{\mathrm{2}}^{\text{'}}与R_{\mathrm{2}}^{″}分别为接触体1与接触体2在接触点的纵向、横向主曲率半径；θ为两接触体纵向曲率半径R_{\mathrm{1}}^{\text{'}}与R_{\mathrm{2}}^{\text{'}}法向平面之间的夹角；X、Y为弹性赫兹接触理论系数.

椭圆偏心率e与系数X、Y之间的关系式为：

\frac{X}{Y}=\left(\mathrm{1}-e^{\mathrm{2}}\right)\frac{F(e)-H(e)}{H(e)-\left(\mathrm{1}-e^{\mathrm{2}}\right)F(e)}（5）

g={\left[\frac{\mathrm{2}H(e)}{\mathrm{\pi }\left(\mathrm{1}-e^{\mathrm{2}}\right)}\right]}^{\frac{\mathrm{1}}{\mathrm{3}}}（6）

h={\left[\frac{\mathrm{2}H(e)}{\mathrm{\pi }}\right]}^{\frac{\mathrm{1}}{\mathrm{3}}}{\left(\mathrm{1}-e^{\mathrm{2}}\right)}^{\frac{\mathrm{1}}{\mathrm{6}}}（7）

{\gamma }_{\mathrm{1}}=\frac{\left(\mathrm{1}-{v_{\mathrm{1}}}^{\mathrm{2}}\right)}{\mathrm{\pi }{E}_{\mathrm{1}}};{\gamma }_{\mathrm{2}}=\frac{\left(\mathrm{1}-{v_{\mathrm{2}}}^{\mathrm{2}}\right)}{\mathrm{\pi }{E}_{\mathrm{2}}}（8）

F(e)={\int }_{\mathrm{0}}^{\frac{\mathrm{\pi }}{\mathrm{2}}}\frac{\mathrm{d}\phi }{\sqrt[]{\mathrm{1}-e^{\mathrm{2}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\phi }}（9）

H(e)={\int }_{\mathrm{0}}^{\frac{\mathrm{\pi }}{\mathrm{2}}}\sqrt[]{\mathrm{1}-e^{\mathrm{2}}\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}\phi }\mathrm{d}\phi（10）

式中E₁、E₂、v₁、v₂分别为两接触体弹性模量与泊松比，F（e）、H（e）是以φ为积分变量的第一、第二类完整椭圆积分.

线接触是两弹性圆柱接触体在受外力作用下的接触情况，当接触线受到挤压时成为面接触，形成长为L，宽为2b₁的矩形接触面积.图3与图4为两接触体在线接触时的接触区域与接触应力分布图.其中ρ₁、ρ₂为两圆柱体接触时的等效曲率半径.

利用赫兹接触理论，计算接触区宽b₁与最大接触应力σ_max.

b_{\mathrm{1}}=\sqrt[]{\frac{\mathrm{4}{F}_{\mathrm{n}}}{\mathrm{\pi }L}\frac{\frac{\mathrm{1}-v_{\mathrm{1}}^{\mathrm{2}}}{E_{\mathrm{1}}}+\frac{\mathrm{1}-v_{\mathrm{2}}^{\mathrm{2}}}{E_{\mathrm{2}}}}{\frac{\mathrm{1}}{{\rho }_{\mathrm{1}}}\pm \frac{\mathrm{1}}{{\rho }_{\mathrm{2}}}}}（11）

{\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}=\sqrt[]{\frac{F_{\mathrm{n}}}{\mathrm{\pi }L}\frac{\frac{\mathrm{1}}{{\rho }_{\mathrm{1}}}+\frac{\mathrm{1}}{{\rho }_{\mathrm{2}}}}{\frac{\mathrm{1}-v_{\mathrm{1}}^{\mathrm{2}}}{E_{\mathrm{1}}}+\frac{\mathrm{1}-v_{\mathrm{2}}^{\mathrm{2}}}{E_{\mathrm{2}}}}}（12）

1.2 点-线弹性接触理论解析计算

在处理齿轮啮合弹性体接触问题时，由于齿轮在啮合过程中横向主曲率半径R_{\mathrm{1}}^{\text{'}}、R_{\mathrm{2}}^{\text{'}}近似为无穷大，因此应将线接触的弹性接触变形问题建立在点接触理论的基础上进行求解，从而根据点线接触转换即可计算齿轮啮合时齿宽为有限长度的轮齿接触界面接触斑大小与长短轴半径^［15］.

将文献［15］中的解析式进行简化得到：

\sum d=\mathrm{2}(X+Y)（13）

\sum d_{\mathrm{1}}=\frac{\mathrm{1}}{R_{\mathrm{1}}^{\text{'}}}+\frac{\mathrm{1}}{R_{\mathrm{2}}^{\text{'}}}（14）

点接触区压力呈半椭球分布，最大接触压力p₀=3F_n/（2πab）.

当接触体横向曲率半径R″→∞时，接触区将变为线接触，采用线接触模型，接触区半宽为：

b_{\mathrm{1}}=\sqrt[]{\frac{\mathrm{4}\left({\gamma }_{\mathrm{1}}+{\gamma }_{\mathrm{2}}\right)q}{\sum d_{\mathrm{1}}}}（15）

计算所得b₁与短半轴解析式，求出线接触时单位长度法向力关系式q：

q=\frac{h^{\mathrm{2}}\left(\mathrm{\Sigma }{d}_{\mathrm{1}}\right){\left\{\frac{\mathrm{3}\mathrm{\pi }\left({\gamma }_{\mathrm{1}}+{\gamma }_{\mathrm{2}}\right)F_{\mathrm{n}}}{\mathrm{2}\sum d}\right\}}^{\frac{\mathrm{2}}{\mathrm{3}}}}{\mathrm{4}\left({\gamma }_{\mathrm{1}}+{\gamma }_{\mathrm{2}}\right)}（16）

根据线接触时接触区压力分布，得出最大接触压力p₁=2q/（πb₁），将q代入该式中与点接触时最大接触压力相比得到：

\frac{P_{\mathrm{1}}}{P_{\mathrm{0}}}=\frac{\mathrm{\pi }}{\mathrm{2}}\cdot \left(\frac{\sum d_{\mathrm{1}}}{\sum d}\right)\cdot g{h}^{\mathrm{2}}（17）

由于齿轮接触时，一对齿轮横向主曲率半径为无穷大即：

\underset{R^{″}\to \mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\left(\frac{\sum d_{\mathrm{1}}}{\sum d}\right)=\mathrm{1}（18）

当R_{\mathrm{2}}^{″}→∞时，b/a→0，通过g与h的极限关系式得到：

\underset{\frac{b}{a}\to \mathrm{0}}{\mathrm{l}\mathrm{i}\mathrm{m}}\left(g{h}^{\mathrm{2}}\right)=\frac{\mathrm{2}}{\mathrm{\pi }}（19）

根据上述关系式得到p₁/p₀=1，求得极限情况下法向力F_n的关系式为：

\frac{F_{\mathrm{n}}}{\mathrm{2}a}=\frac{\mathrm{2}q}{\mathrm{3}}（20）

齿轮接触时为有限长度，代入解析关系式得到有限长线接触长短半轴等效计算参数：

I=\frac{b}{a}=\frac{\mathrm{4}b}{\mathrm{3}L}（21）

式中：当L→∞时，I→0，变为前述接触体无限长线接触情形；当L为定值时，I变为较小数值，可表示有限长线接触.根据所得计算公式，可以计算有限长度接触椭圆偏心率e：

e=\sqrt[]{\mathrm{1}-\frac{b^{\mathrm{2}}}{a^{\mathrm{2}}}}=\sqrt[]{\mathrm{1}-\frac{\mathrm{16}{b}^{\mathrm{2}}}{\mathrm{9}{L}^{\mathrm{2}}}}（22）

1.3 赫兹非线性法向接触刚度

计算接触变形δ与法向力F_n之间的解析关系式，同时根据刚度定义给出符合赫兹非线性接触条件通用法向接触刚度计算公式^［16］：

F_{\mathrm{n}}=K_{\mathrm{n}\mathrm{r}}{\delta }^{\frac{\mathrm{3}}{\mathrm{2}}}（23）

K_{\mathrm{n}\mathrm{r}}=\frac{\mathrm{4}}{\mathrm{3}}\frac{q_{\mathrm{k}}}{\left({\gamma }_{\mathrm{1}}+{\gamma }_{\mathrm{2}}\right)\sqrt[]{X+Y}}（24）

K_{\mathrm{h}\mathrm{n}}=\frac{\mathrm{d}{F}_{\mathrm{n}}}{\mathrm{d}\delta }=\frac{\mathrm{3}}{\mathrm{2}}{K}_{\mathrm{n}\mathrm{r}}{\delta }^{\frac{\mathrm{1}}{\mathrm{2}}}=\frac{\mathrm{3}}{\mathrm{2}}{K_{\mathrm{n}\mathrm{r}}}^{\frac{\mathrm{2}}{\mathrm{3}}}{F_{\mathrm{n}}}^{\frac{\mathrm{1}}{\mathrm{3}}}（25）

式中：q_k=\sqrt[]{\mathrm{\pi }}［g/2F（e）］；K_nr为赫兹接触常数.从解析关系式看出法向接触刚度不仅与接触体的材料、椭圆接触斑大小以及主曲率半径有关，同时也与法向载荷相关，当载荷增大时法向接触刚度也随之增大.

2 齿轮非线性时变啮合刚度计算

本文采用赫兹非线性接触理论得到接触刚度，考虑不同负载情况下刚度非线性变化趋势，计算齿轮的综合时变啮合刚度.

2.1 时变啮合刚度修正算法

齿轮在啮合力作用下产生变形，采用材料力学中悬臂梁模型，基于Weber能量法计算齿轮弯曲变形、剪切变形、轴向压缩变形，进而计算齿轮时变啮合刚度.齿轮刚度计算模型如图5所示.

给出齿轮弯曲、剪切、轴向压缩变形计算解析式：

\begin{array}{l}{\delta }_{\mathrm{b}}=\frac{F_j}{E_{\mathrm{e}}}\left[\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}{\beta }_j\left(\frac{L_i^{\mathrm{3}}+\mathrm{3}{L}_i^{\mathrm{2}}{S}_{ij}+\mathrm{3}{L}_i{S}_{ij}^{\mathrm{2}}}{\mathrm{3}{\overline{I}}_i}\right)\right.-\\ \left.\mathrm{c}\mathrm{o}\mathrm{s}{\beta }_j\mathrm{s}\mathrm{i}\mathrm{n}{\beta }_j\left(\frac{L_i^{\mathrm{2}}{Y}_j+\mathrm{2}{L}_i{Y}_j{S}_{ij}}{\mathrm{2}{\overline{I}}_i}\right)\right]\end{array}（26）

{\delta }_{\mathrm{s}}=\frac{F_j}{E_{\mathrm{e}}}\left[\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}{\beta }_j\left(\frac{\mathrm{12}(\mathrm{1}+\nu )L_i}{\mathrm{5}{\overline{A}}_i}\right)\right]（27）

{\delta }_{\mathrm{c}}=\frac{F_j}{E_{\mathrm{e}}}\left[\mathrm{s}\mathrm{i}{\mathrm{n}}^{\mathrm{2}}{\beta }_j\left(\frac{L_i}{{\overline{A}}_i}\right)\right]（28）

式中：F_j 为任意啮合点的载荷；Y_j 与S_ij 分别表示载荷作用点在y轴与微元在x轴方向上投影的距离；L_i 为微元厚度；β_j 是F_j 与y轴之间的载荷角；\nu为泊松比.

{\overline{I}}_i=\frac{\mathrm{2}{I}_i{I}_{i+\mathrm{1}}}{I_i+I_{i+\mathrm{1}}}（29）

I_i=\frac{\mathrm{1}}{\mathrm{12}}B{\left(\mathrm{2}{h}_i\right)}^{\mathrm{3}}（30）

{\overline{A}}_i=\frac{\mathrm{2}{A}_i{A}_{i+\mathrm{1}}}{A_i+A_{i+\mathrm{1}}}（31）

A_i=\mathrm{2}B{h}_i（32）

（33）

式中：E为弹性模量；B为齿宽；H_p为节圆齿厚；h_i 表示微元到x轴之间的距离；A_i 与I_i 则分别表示为当量截面面积与当量惯性矩.当B与H_p之比小于5时为窄齿，其余则为宽齿.

计算齿轮啮合过程基体变形^［17］：

对于窄齿，有

{\delta }_{\mathrm{f}}=\frac{F_{\mathrm{n}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}{\beta }_j}{BE}\left[\mathrm{1.534}\left(\mathrm{1}+\frac{\mathrm{0.416}\mathrm{7}\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}{\beta }_j}{\mathrm{1}+\nu }\right)\right.+

\left.\mathrm{5.306}{\left(\frac{L_{\mathrm{f}}}{H_{\mathrm{f}}}\right)}^{\mathrm{2}}+\mathrm{2}(\mathrm{1}-\nu )\left(\frac{L_{\mathrm{f}}}{H_{\mathrm{f}}}\right)\right]（34）

对于宽齿，有

{\delta }_{\mathrm{f}}=\frac{F_{\mathrm{n}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{2}}{\beta }_j}{BE}\left(\mathrm{1}-{\nu }^{\mathrm{2}}\right)\left[\mathrm{2}\left(\frac{\mathrm{1}-\nu -\mathrm{2}{\nu }^{\mathrm{2}}}{\mathrm{1}-{\nu }^{\mathrm{2}}}\right)\frac{L_{\mathrm{f}}}{H_{\mathrm{f}}}\right.+

\left.\mathrm{5.306}{\left(\frac{L_{\mathrm{f}}}{H_{\mathrm{f}}}\right)}^{\mathrm{2}} +\mathrm{1.534}\left(\mathrm{1}+\frac{\mathrm{0.416}\mathrm{7}\mathrm{t}\mathrm{a}{\mathrm{n}}^{\mathrm{2}}{\beta }_j}{\mathrm{1}+\nu }\right)\right]（35）

由有效长度L_e，计算求出L_f和H_f：

\begin{array}{l}{L}_{\mathrm{f}}=X_j-X_M-Y_j\mathrm{t}\mathrm{a}\mathrm{n}{\beta }_j\\ H_{\mathrm{f}}=\mathrm{2}{Y}_M\end{array}（36）

式中：X_j 表示载荷作用点在x轴方向上的距离；X_M 与Y_M 为基点M位于横纵坐标轴上的大小.

给出轮齿接触变形传统计算公式^［18］：

{\delta }_{\mathrm{p}}=\frac{\mathrm{2}{F}_{\mathrm{n}}(\mathrm{1}-{\nu }^{\mathrm{2}})}{\mathrm{\pi }BE}\left(-\mathrm{0.429}+\mathrm{l}\mathrm{n}\frac{\mathrm{4}{\rho }_{\mathrm{1}}{\rho }_{\mathrm{2}}}{b^{\mathrm{2}}}\right)（37）

由式（37）可知，目前接触变形与载荷为线性关系，未能表征不同载荷作用下接触刚度的非线性变化规律.因此，本文综合考虑齿轮啮合过程中接触点曲率半径、载荷、接触斑大小、接触宽度及齿面材料等因素，计算非线性赫兹接触变形δ_p.

综上，求得单对轮齿啮合刚度K：

K=F_j/\left[{\left({\delta }_{\mathrm{b}}+{\delta }_{\mathrm{s}}+{\delta }_{\mathrm{c}}+{\delta }_{\mathrm{f}}\right)}_{\mathrm{1}}+{\delta }_p+{\left({\delta }_{\mathrm{b}}+{\delta }_{\mathrm{s}}+{\delta }_{\mathrm{c}}+{\delta }_{\mathrm{f}}\right)}_{\mathrm{2}}\right]（38）

2.2 非线性法向接触刚度

采用建立的考虑非线性赫兹接触条件下轮齿啮合接触刚度算法，分析不同负载条件下啮合刚度变化特性.引用文献中所用齿轮副参数与材料属性如表1所示^［19-20］.

图6表示啮入至啮出过程中齿轮副曲率半径随啮合角变化的规律曲线图，可以看出：由于主、从动轮齿数关系及几何参数不同，两齿轮副啮合曲率半径变化趋势也不同.齿轮副1主动轮齿数小于从动轮，呈上升趋势，曲率半径随啮合角的增大而增大；而齿轮副2主动轮齿数大于从动轮，呈下降趋势，曲率半径在啮合转角为0.32 rad达到最大后逐渐减小.

将表1中的参数代入齿轮法向接触刚度计算，结果如图7所示.

负载条件下不同齿轮副的法向接触刚度与传统赫兹接触刚度的对比如图7所示，与传统赫兹接触刚度为定值相比，法向接触刚度会随着啮合角的变化而变化，接触刚度曲线与综合曲率半径变化规律一致.随转矩递增，法向接触刚度非线性增大.当转矩从100 N·m增大到2 000 N·m时，齿轮副1与齿轮副2从啮入至啮出过程中法向接触刚度分别增加了（1.39~2.58）×10⁹ N/m与（3.10~2.48）×10⁹ N/m.

2.3 负载条件下时变啮合刚度

计算齿轮副1与齿轮副2不同输入转矩下齿轮时变啮合刚度，结果如图8所示.

图8为输入不同转矩后时变啮合刚度变化规律曲线，与文献［19-20］中所采用的刚度算法的结果一致，验证了本文计算方法的准确性.可以看出：随着输入转矩增大，齿轮副1与齿轮副2啮合刚度非线性增大.当转矩从100 N·m增大到2 000 N·m时，齿轮副1单、双齿啮合刚度分别增加1.8×10⁷ N/m和2.9×10⁷ N/m，齿轮副2单、双齿啮合刚度分别增加3.8×10⁷ N/m和2.7×10⁷ N/m.随输入转矩递增，齿轮时变啮合刚度增幅逐渐减小，这是由于负载增大使轮齿受载压平，轮齿接触界面所受接触面积的增大能力逐渐减弱，齿轮变形受弹性变形约束的影响变化幅度逐渐减小.当载荷增大到一定值时，啮合刚度基本保持不变.

3 计入冲击的齿轮副时变啮合刚度解析算法

本文基于冲击理论给出冲击过程中任意点接触位置与载荷角的变化解析关系式，并根据冲击时刻动态接触关系，探究冲击对刚度影响的确切变化规律.

3.1 线外冲击位置的确定与离散化

确定线外啮入点冲击位置后，基于圆柱碰撞理论可以计算出冲击力大小，线外啮合冲击位置如图9所示^［21］.

确定线外冲击点位置后，计算冲击过程中任意啮合点坐标，本节给出任意点位置计算方法.

如图10所示，L₁、L₂为齿轮副理论啮合点；B₁、B₂为齿轮副实际啮合点；N、S分别为主、从动轮在啮合线单齿啮合区间所对应的点.

理想情况下，前一齿位于N点时，后一齿则在S点开始进入啮合.而考虑线外啮合后，当前一齿在N^{\text{'}}点时，后一齿已经在线外啮入点C产生冲击作用.因此将C点等效为啮合线上的B_{\mathrm{2}}^{\text{'}}点，根据冲击理论求出主动轮从C点到B₂点所转过的角度θ，设B_{\mathrm{2}}^{\text{'}}B₂与N^{\text{'}}N距离相等，则\theta×r_b1=B_{\mathrm{2}}^{\text{'}}B₂=N^{\text{'}}N，将B_{\mathrm{2}}^{\text{'}}B₂离散后不同的等效啮合点来确定与之对应啮合冲击点位置，该方法即可用于后续任意线外啮合点冲击位置到圆心的半径与载荷角的计算，等效后实际啮合线将变为B_{\mathrm{2}}^{\text{'}}B₁.

3.2 任意冲击位置半径计算

主、从动轮线外啮合冲击数学模型如图11所示.其中r_a1、r_a2为主、从动轮齿顶圆半径，r_b1、r_b2为主、从动轮基圆半径，r_{\mathrm{b}\mathrm{2}}^{\text{'}}为线外啮合时刻从动轮假想基圆半径.齿轮副在冲击时刻从动轮啮合点到圆心的距离始终为齿顶圆半径r_a2.

\mathrm{\angle }R{O}_{\mathrm{1}}P=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}\frac{O_{\mathrm{1}}{C}^{\mathrm{2}}+{\left(r_{\mathrm{1}}+r_{\mathrm{2}}\right)}^{\mathrm{2}}-{r_{\mathrm{a}\mathrm{2}}}^{\mathrm{2}}}{\mathrm{2}\times \left(r_{\mathrm{1}}+r_{\mathrm{2}}\right)\times O_{\mathrm{1}}C}-\frac{\Delta }{r_{\mathrm{b}\mathrm{1}}}（39）

\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{\angle }R{O}_{\mathrm{1}}P=\frac{O_{\mathrm{1}}{R}^{\mathrm{2}}+{\left(r_{\mathrm{1}}+r_{\mathrm{2}}\right)}^{\mathrm{2}}-{r_{\mathrm{a}\mathrm{2}}}^{\mathrm{2}}}{\mathrm{2}\times \left(r_{\mathrm{1}}+r_{\mathrm{2}}\right)\times O_{\mathrm{1}}R}（40）

式中：O₁C为初始啮合点半径；r₁与r₂为主、从动轮分度圆半径；Δ为冲击时转动角度在啮合线上的等效距离.结合上式可得到任意冲击时刻主动轮啮合半径O₁R.

3.3 线外啮合冲击载荷角计算

载荷角区与压力角表示不同，载荷角是指轮齿对称线垂线与载荷方向之间的夹角，本节给出主、从动轮冲击时载荷角计算方法.同时参考图11与 图12，其中α₁、α₂分别为线外啮合时主、从动轮啮合点压力角，θ₁、θ₂分别为主、从动轮啮合半径与中心线之间的夹角，β₁为主动轮线外啮合载荷角.

\left\{\begin{array}{l}{\alpha }_{\mathrm{1}}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}\frac{r_{\mathrm{b}\mathrm{1}}}{O_{\mathrm{1}}R}\\ {\alpha }_{\mathrm{2}}={\alpha }_{\mathrm{1}}+{\theta }_{\mathrm{1}}+{\theta }_{\mathrm{2}}\end{array}\right.（41）

{\theta }_{\mathrm{1}}=\mathrm{\angle }R{O}_{\mathrm{1}}P=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}\frac{{\left(r_{\mathrm{1}}+r_{\mathrm{2}}\right)}^{\mathrm{2}}+O_{\mathrm{1}}{R}^{\mathrm{2}}-{r_{\mathrm{a}\mathrm{2}}}^{\mathrm{2}}}{\mathrm{2}\times \left(r_{\mathrm{1}}+r_{\mathrm{2}}\right)\times O_{\mathrm{1}}R}（42）

{\theta }_{\mathrm{2}}=\mathrm{\angle }R{O}_{\mathrm{2}}P=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}\frac{{\left(r_{\mathrm{1}}+r_{\mathrm{2}}\right)}^{\mathrm{2}}+{r_{\mathrm{a}\mathrm{2}}}^{\mathrm{2}}-O_{\mathrm{1}}{R}^{\mathrm{2}}}{\mathrm{2}\times \left(r_{\mathrm{1}}+r_{\mathrm{2}}\right)\times r_{\mathrm{a}\mathrm{2}}}（43）

\left\{\begin{array}{l}{\beta }_{\mathrm{1}}={\alpha }_{\mathrm{1}}-{\phi }_{\mathrm{2}}={\alpha }_{\mathrm{1}}-\left({\phi }_{\mathrm{1}}+{\alpha }_{\mathrm{f}}\right)\\ {\alpha }_{\mathrm{f}}=\mathrm{i}\mathrm{n}\mathrm{v}\alpha -\mathrm{i}\mathrm{n}\mathrm{v}{\alpha }_{\mathrm{1}}\end{array}\right.（44）

式中：α_f为线外啮合点与分度圆的夹角；分度圆与对称线夹角φ₁=π/（2z₁）；φ₂为线外啮合点与对称线的夹角.根据上述关系式即可计算冲击时刻主动轮载荷角β₁，相同方法求得从动轮载荷角β₂.

3.4 计入冲击与受载的齿轮副时变啮合刚度结果分析

考虑线外啮合冲击时刻接触位置、载荷角及接触曲率等变化，综合分析计入冲击与受载对不同齿轮副时变啮合刚度的影响，解析计算如图13所示.

图13为计入冲击时不同输入转矩的时变啮合刚度曲线图.可见：随着输入转矩增大，齿轮副1与齿轮副2的啮合刚度非线性增大，且由轻载至重载变化过程中因加剧了冲击效应，刚度产生硬化，使得刚度突变值更大.考虑冲击时初始啮合点刚度高于理论值，啮合刚度曲线呈现由大到小的变化规律，随着齿轮啮合平稳后趋于正常.当转矩从100 N·m增大到2 000 N·m时，齿轮副1与齿轮副2因实际重合度增大使单齿啮合区间减小，单齿啮合时间缩短为正常啮合时的33%与59%.

4 结 论

本文给出满足赫兹非线性接触理论的弹性体接触斑大小和刚度计算方法.建立了极限情况下齿轮副点线接触转换模型，并根据轮齿接触变形、材料属性以及接触几何及载荷之间的相互关系，提出了考虑负载条件下时变啮合刚度修正算法；而后给出了一种考虑线外啮合冲击任意时刻冲击半径与载荷角的数学模型，根据冲击时刻动态接触关系，探究不同负载下冲击对刚度的影响，得出以下结论.

1）与传统赫兹接触刚度为定值相比，接触刚度规律变化不同，接触刚度曲线与齿轮综合主曲率半径变化规律一致，当不断增大输入转矩后，接触刚度递增且增大幅度逐渐减小.

2）将所得齿轮法向接触刚度代入时变啮合刚度中计算，负载增大齿轮啮合刚度非线性增大，当载荷增大到一定值时，齿轮啮合刚度基本保持不变.

3）综合分析计入线外啮合冲击时刚度产生突变，刚度会随着时间变化呈现由高到低的变化规律.随着负载递增，刚度产生硬化，冲击效应增强，齿轮单齿啮合时间因实际重合度增大而减小.刚度的变化会进一步引起齿轮产生振动与噪声，上述研究结果为齿轮系统动力学响应精准分析提供了基础.

参考文献

[1] ANDERSSON A， VEDMAR L ．A dynamic model to determine vibrations in involute helical gears［J］．Journal of Sound and Vibration，2003，260（2）：195-212．

[2] WANG J， HOWARD I. The torsional stiffness of involute spur gears［J］. Proceedings of the Institution of Mechanical Engineers，Part C： Journal of Mechanical Engineering Science， 2004， 218（1）： 131-142．

[3] 唐进元，蒲太平. 基于有限元法的螺旋锥齿轮啮合刚度计 算［J］. 机械工程学报， 2011， 47（11）： 23-29．

[4] TANG J Y， PU T P ．Spiral bevel gear meshing stiffness calculations based on the finite element method［J］．Journal of Mechanical Engineering， 2011，4 7（11）： 23-29．（in Chinese）

[5] 刘文，李锐，张晋红， 等 ．斜齿轮时变啮合刚度算法修正及影响因素研究［J］．湖南大学学报（自然科学版）， 2018，45（2）：1-10．

[6] LIU W， LI R， ZHANG J H，et al ．Study on correction algorithm of time-varying mesh stiffness of helical gears and its influencing factors［J］．Journal of Hunan University （Natural Sciences），2018，45（2）：1-10．（in Chinese）

[7] 林腾蛟， 赵子瑞， 江飞洋， 等 ．考虑温度效应的斜齿轮时变啮合刚度解析算法［J］．湖南大学学报（自然科学版）， 2020， 47（2）：6-13．

[8] LIN T J， ZHAO Z R， JIANG F Y，et al ．An analytic algorithm of time-varying mesh stiffness of helical gears considering temperature effect［J］．Journal of Hunan University （Natural Sciences）， 2020， 47（2）： 6-13．（in Chinese）

[9] 陈思宇，谭儒龙，郭晓东，等 ．基于能量等效的直齿锥齿轮时变啮合刚度计算［J］．航空动力学报，2022，37（2）：320-329．

[10] CHEN S Y， TAN R L， GUO X D， et al ．Time-varying meshing stiffness calculation of straight bevel gears based on energy equivalence［J］．Journal of Aerospace Power， 2022， 37（2）： 320-329．（in Chinese）

[11] LIU Y H， SHI Z Q， SHEN G J，et al ．Evaluation model of mesh stiffness for spur gear with tooth tip chipping fault［J］．Mechanism and Machine Theory， 2021， 158： 104238．

[12] CHEN X Q， WEI J， ZHANG J X，et al ．A novel method to reduce the fluctuation of mesh stiffness by high-order phasing gear sets：theoretical analysis and experiment［J］．Journal of Sound and Vibration，2022，524：116752．

[13] SUN Z， CHEN S Y， HU Z H，et al ．Analytical models for thermal deformation and mesh stiffness of spur gears under steady temperature field［J］．Engineering Failure Analysis，2022，133：105972．

[14] 张伟，张学良， 温淑花， 等 ．考虑微凸体基体变形和相互作用的结合面法向接触刚度模型［J］．西安交通大学学报， 2020，54（6）： 115-121．

[15] ZHANG W， ZHANG X L， WEN S H，et al ．A normal contact stiffness model of joint surfaces considering interaction of deformations of substrate and asperity［J］．Journal of Xi’an Jiaotong University， 2020， 54（6）： 115-121．（in Chinese）

[16] 兰国生，孙万，谭文兵，等 ．基于圆锥微凸体的结合面法向刚度分形模型研究［J］．振动与冲击，2021，40（15）：207-215．

[17] LAN G S， SUN W， TAN W B，et al ．Cone elastoplastic fractal model of two contact rough surfaces［J］．Journal of Vibration and Shock，2021，40（15）： 207-215．（in Chinese）

[18] 兰国生，冀成龙，李祥，等 ．混合润滑结合面法向接触刚度模型研究［J］．振动与冲击， 2023， 42（17）： 220-227．

[19] LAN G S， JI C L， LI X，et al ．Normal contact stiffness model of mixed lubrication joint surface［J］. Journal of Vibration and Shock， 2023， 42（17）： 220-227．（in Chinese）

[20] 武宝林， 杨素君， 姚俊红. 齿轮传动中啮合冲击的理论分析［J］. 机械科学与技术， 2003， 22（1）： 55-57．

[21] WU B L， YANG S J， YAO J H ．Theoretical analysis on meshing impact of involute gears［J］. Mechanical Science and Technology，2003，22（1）： 55-57．（in Chinese）

[22] HE Z Y， ZHANG T， LIN T J ．Novel mathematical modelling method for meshing impact of helical gear［J］. Mechanism and Machine Theory， 2020， 152： 103949．

[23] 丁长安， 张雷， 周福章， 等. 线接触弹性接触变形的解析算 法［J］.摩擦学学报， 2001， 21（2）： 135-138.

[24] DING C A， ZHANG L， ZHOU F Z， et al. Theoretical formula for calculation of line-contact elastic contact deformation［J］. Tribology， 2001，21（2）： 135-138. （in Chinese）

[25] 关庆华，赵鑫，温泽峰，等 ．基于Hertz接触理论的法向接触刚度计算方法［J］．西南交通大学学报，2021，56（4）：883-890．

[26] GUAN Q H， ZHAO X， WEN Z F，et al ．Calculation method of Hertz normal contact stiffness［J］．Journal of Southwest Jiaotong University，2021，56（4）：883-890．（in Chinese）

[27] CORNELL R W. Compliance and stress sensitivity of spur gear teeth［J］. Journal of Mechanical Design， 1981， 103（2）：447-459．

[28] 张敬彩，刘红旗，赵向飞，等 ．直齿轮的轮齿变形分析及其快捷计算方法［J］．机械传动，2015，39（10）：48-55．

[29] ZHANG J C， LIU H Q， ZHAO X F，et al ．Theoretical analysis and quick calculation of spur gear tooth deformation［J］．Journal of Mechanical Transmission，2015，39（10）：48-55．（in Chinese）

[30] 张建宇，刘鑫博，张飞斌 ．裂纹齿轮啮合刚度的劣化特性［J］．振动·测试与诊断， 2014， 34（6）： 1124-1128．

[31] ZHANG J Y， LIU X B， ZHANG F B ．Degradation characteristics of meshing stiffness of spur gears with root crack［J］．Journal of Vibration，Measurement & Diagnosis，2014，34（6）：1124-1128．（in Chinese）

[32] 杨勇，王家序，周青华，等 ．考虑摩擦的磨损和修形齿轮啮合刚度计算［J］．工程科学与技术，2018，50（2）：212-219．

[33] YANG Y， WANG J X， ZHOU Q H，et al ．Stiffness calculation considering friction for a spur gear pair with tooth wear and profile modification［J］．Advanced Engineering Sciences，2018，50（2）：212-219．（in Chinese）

[34] 王承登， 何泽银，杨震， 等 ．受载齿轮副啮入冲击激励计算方法［J］．航空动力学报，2024，39（5）：239-248．

[35] WANG C D， HE Z Y， YANG Z，et al ．Calculation method of meshing impact excitation of loaded gear pair［J］．Journal of Aerospace Power，2024， 39（5）： 239-248．（in Chinese）

基金资助

国家重点研发计划资助项目(2023YFB3406301)

National Key Research and Development Program(2023YFB3406301)

重庆市自然科学基金资助项目(cstc2021jcyj-msxmX0562)

Chongqing Natural Science