齿轮传动广泛应用于船舶、航空航天及风力发电等领域，齿轮系统的动力学特性将直接影响机械结构的安全性和可靠性.最近几年来，国内外学者对齿轮系统动力学等内容开展了深入的研究^[1-2]. Kahraman等^[3-4]建立了单自由度直齿轮非线性动力学模型，考虑周期时变刚度的影响，用数值法得到了系统的动力学特性. Chen等^[5]推导了齿轮中心距误差与齿隙之间的关系，定义了动态齿隙，并对计及动态间隙、摩擦和时变刚度的多自由度非线性动态齿轮传动系统进行了研究. Wu等^[6]建立了复合行星齿轮非线性动力学模型，讨论了因磨损导致的齿侧间隙变化对系统的影响.

在以往研究中，因研究条件限制而忽略由制造误差和磨损润滑等造成的不确定性因素，但由于上述因素对系统响应有重要的影响，学者们开始重视对不确定性因素的研究^[7].王靖岳等^[8]通过正态分布规律来描述啮合刚度的不确定性，并运用4阶Runge-Kutta数值积分法对三自由度直齿轮副系统的稳定性问题进行了分析和研究.廖映华等^[9]则在啮合刚度的确定性表达式中增加随机分量，通过服从截尾正态分布规律的随机变量描述其随机不确定性.冯松等^[10]根据Timoshinko梁理论，分析和计算了点蚀和磨损情况下的齿轮啮合刚度以及其发展规律，但只是分析了刚度，并没有研究刚度对于整个齿轮-轴承系统的动力学特性的影响.陈奇等^[11]引入可用于接触强度分析的M-B分形接触新理论模型，建立两粗糙曲面的分形接触模型，为利用分形接触模型研究齿轮接触动力学提供了参考.

基于考虑齿侧间隙的齿轮-轴承系统模型^[12]，本文应用分形理论描述齿轮轮廓，计算和分析不同分形维数D下的时变啮合刚度，研究考虑齿面分形特性的齿轮啮合刚度对齿轮-轴承系统动力学特性的影响.旨在通过研究，探明齿侧间隙和时变刚度等因素对齿轮啮合的动力学特性的影响，为精确建模、齿轮设计和制造提供参考依据.

1 系统模型 1.1 齿轮-轴承系统

采用考虑齿侧间隙的齿轮-轴承系统模型^[12]如图 1所示，模型中齿轮为渐开线直齿轮，轴承为滑动轴承.考虑时变啮合刚度、综合传递误差、齿侧间隙以及滑动轴承非线性油膜力等因素，并考虑轴的径向振动以及两齿轮的扭转振动，但不考虑齿面摩擦力的影响，建立坐标系x_iO_iy_i (i=1或2，表示主、从动轮)，其中O₁、O₂为两轴承几何中心；O_p、O_g为主从动齿轮几何中心.

系统的无量纲动力学方程采用文献[12]中的式(6)；k_h(t)为时变啮合刚度，综合传递误差e(t)=ecos(ω_ht)；齿轮副间隙非线性函数

$ {\bar f_{\rm{b}}} = \left\{ {\begin{array}{*{20}{c}} {{{\bar x}_{\rm{h}}} - {{\bar b}_{\rm{h}}}, }&{{{\bar x}_{\rm{h}}} > {{\bar b}_{\rm{h}}};}\\ {0, }&{\left| {{{\bar x}_{\rm{h}}}} \right| \le {{\bar b}_{\rm{h}}};}\\ {{{\bar x}_{\rm{h}}} + {{\bar b}_{\rm{h}}}, }&{{{\bar x}_{\rm{h}}} < - {{\bar b}_{\rm{h}}}.} \end{array}} \right. $

模型中，应用Capone^[13]修正模型计算滑动轴承的非线性油膜力，滑动轴承的油膜力分布一般由雷诺方程确定，无量纲的油膜力分量如文献[12]中的式(3)所示.

1.2 考虑齿面分形特性的齿轮啮合刚度

微观上，齿轮的表面因呈现出随机性、多尺度性和自仿射性，具有分形的基本特征，因此可采用分形理论来描述齿轮表面^{[12, 14-15]}.

以渐开线直齿圆柱齿轮为具体对象，建立齿轮齿廓方程^[10]：

$ \left\{ {\begin{array}{*{20}{l}} {{x_i} = {r_i}\cos {\gamma _i}, }\\ {{y_i} = {r_i}\sin {\gamma _i}.} \end{array}} \right. $

(1)

式中：γ_i=s/(2r)－(inv α_i－inv α₀)；s为分度圆齿厚；r为分度圆半径；r_i为渐开线齿廓上任一点的回转半径；α₀为分度圆压力角，按常用习惯取α₀=20°；α_i为渐开线齿廓上任一点对应的压力角.使用W-M模型^[16]描述齿轮表面的形貌：

$z(x) = L{\left( {\frac{G}{L}} \right)^{D - 1}}\sum\limits_{n = 0}^\infty {\frac{{\cos \left[ {\left( {2\pi {\gamma ^n}x} \right)/L} \right]}}{{{\gamma ^{(2 - D)n}}}}} . $

(2)

式中:z(x)为齿轮形貌高度, x为轮廓的位置坐标, D为分形维数, G为特征尺度系数, γⁿ齿轮轮廓的空间频率, L为取样长度.

综上，根据式(1)和式(2)可得具有分形特性的齿轮齿廓形貌高度(图 2)和考虑齿面分形特性的齿轮齿廓的渐开线方程，如式(3)所示.从图 2中可以看出，随着分形维数的增加，齿轮形貌高度的幅值变小，并且表现出齿轮粗糙齿面中更多的细节.同时由于齿面形貌的变化，在齿轮啮合的过程中使得轮齿变形量发生变化，进而导致接触刚度发生变化.

$ \left\{ {\begin{array}{*{20}{l}} {{x_i} = {r_i}\cos {\gamma _i}, }\\ {{y_i} = {r_i}\sin {\gamma _i} + z\left( {{x_i}} \right).} \end{array}} \right. $

(3)

根据式(3)，可得考虑齿面分形特性的粗糙齿面轮廓线，为了更加明显地表现考虑齿面分形特性的齿廓，对其形貌高度进行放大处理，得到如图 3所示的与光滑齿面轮廓线对比示意图.

采用Weber-Banaschek方法^[17]，考虑轮齿的弯曲、剪切变形δ_Z，轮齿基础部分的变形δ_R和齿轮的接触变形δ_pW，计算齿轮的啮合刚度.

$ \begin{array}{*{20}{c}} {{\delta _{\rm{Z}}} = \frac{{{F_{\rm{N}}}}}{{Eb}}{{\cos }^2}{\alpha _x}\left\{ {10.92\int_0^{{h_x}} {\frac{{{{\left( {{h_x} - y} \right)}^2}}}{{{{(2x)}^3}}}} {\rm{d}}y + } \right.}\\ {3.1\left( {1 + 0.294{{\tan }^2}{\alpha _x}} \right)\int_0^{{h_x}} {\frac{{{\rm{d}}y}}{{2x}}} \} , }\\ {{\delta _{\rm{R}}} = \frac{{{F_{\rm{N}}}}}{{Eb}}{{\cos }^2}{\alpha _x}\left\{ {5.2\frac{{h_x^2}}{{S_{\rm{F}}^2}} + \frac{{{h_x}}}{{{S_{\rm{F}}}}} + 1.4\left( {1 + 0.294{{\tan }^2}{\alpha _x}} \right)} \right\}.} \end{array} $

式中:x、y分别为轮齿啮合法向力作用点处的坐标，α_x为作用点的啮合角，F_N为接触点压力，E为综合弹性模量，b为齿宽.

接触产生的变形量如式(4)所示，a为赫兹接触半径，其表达式如式(5)所示，进而得到齿轮的啮合刚度k，如式(6)所示.

$ {\delta _{{\rm{pW}}}} = 0.58\frac{{{F_{\rm{N}}}}}{{Eb}}\left( {{\mathop{\rm In}\nolimits} \frac{{2{h_1}}}{a} + {\mathop{\rm In}\nolimits} \frac{{2{h_2}}}{a} - 0.429} \right), $

(4)

$ a = 1.52\sqrt {\frac{{{F_{\rm{N}}}}}{{Eb}}} \sqrt {\frac{{{\rho _1}{\rho _2}}}{{{\rho _1} + {\rho _2}}}} , $

(5)

$ \begin{array}{*{20}{c}} {{\delta _\Sigma } = {\delta _{{\rm{Z}}1}} + {\delta _{{\rm{Z}}2}} + {\delta _{{\rm{R}}1}} + {\delta _{{\rm{R}}2}} + {\delta _{{\rm{pW}}}}, }\\ {k = {F_{\rm{N}}}/{\delta _\Sigma }.} \end{array} $

(6)

式中，ρ₁和ρ₂分别为两个齿轮轮齿的齿面在接触点处的曲率半径，其他具体参数如图 4所示.

由式(6)可知，刚度是轮廓位置坐标的函数，将随其变化而发生变化，又由于齿轮转速将轮廓位置坐标和时间相关联，故，刚度将随时间变化而变化，具有时变特性.

综上，得到不考虑齿面分形特性和不同分形维数D下的齿轮啮合刚度，如图 5所示.从图 5中可以看出，随着分形维数的增加，齿轮的时变啮合刚度的整体趋势相同，但是在考虑分形理论之后，齿轮轮廓的形貌高度随轮廓的位置发生变化，啮合刚度随之产生波动，并且随着分形维数D的增加，齿轮表面越来越精细，形貌高度变小，刚度波动减小，趋向图 5(a)中所示正常刚度.另外，不同分形维数D下的刚度的平均值如表 1所示，随着D的减小，齿轮齿面粗糙度增大，啮合刚度降低，与文献[18]中结论相似.

2 数值仿真

将时变啮合刚度模型代入考虑齿面分形特性的齿轮系统，采用4-5阶Runge-Kutta方法进行仿真分析，可得到系统的分岔图、三维频谱图、时域图以及均方根图等.为了验证其仿真结果的准确性，由于目前国内外公开发表的文献中未能找到与其相关的可参考信息，后续将采用利用齿轮的振动位移等效啮合点位移的方法进行实验分析验证.

2.1 分形维数D对系统响应的影响

为了研究分形维数D对系统响应的影响，取系统参数k_i=1.2，b_h=0.2 mm，e_p=e_g=0.01，e=0.1，ζ_h=0.01，ζ₁=ζ₂=ζ_p=ζ_g=0.01，F_m=0.105，F_a=0，D=1.5，关于齿轮-轴承系统的其他结构参数如表 2所示，仿真结果如图 6所示.

由图 6(b)中可以看出，当ω < 1.5时，系统处于单周期或双周期运动，系统的运动状态比较稳定；当1.5 < ω < 2.1时，系统进入运动区域A₁，单周期分岔为三周期，之后进入混沌，持续一段时间后又变为三周期运动，最后回归稳定的单周期运动，随之进入运动区域A₂；当ω > 2.4时，先后经历两周期、拟周期、四周期、两周期、拟周期和混沌等运动，最终回归单周期运动.从图 6(b)的三维频谱中，可以看出，在相对低速区域A₁，系统的频率成分主要集中在1/3f_m，2/3f_m，f_m和2f_m次谐响应；在相对高速区域A₂，系统的频率成分主要集中在1/2f_m和f_m次谐响应.

从图 6可以看出，当D=1.1时，由于刚度的波动加剧，系统振动的相对位移稍变大，运动区域A₁和A₂中的周期运动变为拟周期运动，并且进入混沌区域的时间点提前，混沌和拟周期运动区域增大，系统的不稳定性增大.

为了研究单一转速下，刚度对系统响应的影响情况，令ω=1.7，k_i=1.6，其余参数保持不变.为了凸显分形维数D的变化对系统的影响，更明显地进行对比分析，只列举D分别取1.1、1.5和1.9时，系统的相图和Poincaré截面图，如图 7所示；D分别取1.1、1.3、1.5和1.9时，系统的时域图及幅值均方根图如图 8和图 9所示.

从图 7中可以看出，随着D的减小，齿轮系统从标准的三周期运动向拟周期运动过渡，Poincaré截面图中由最开始的3点到后来的3簇，逐渐扩散，系统趋向于不稳定的拟周期运动.

从图 8中可以看出，随着分形维数D的增大，系统振幅变化的整体趋势基本相同，其中D=1.1和D=1.3的幅值波动情况更为明显，结合表 3所示的振动位移均方根值可以发现，虽然D=1.5和D=1.9时的时域图变化不大，但D=1.5时的均方根值稍大，也就是说，随着D的减小，齿面更加精细，系统振动位移减小，系统运动相对更加稳定.

从图 9(a)中可以看出，当D＝1.1和D=1.3时的均方根值变化幅度较大，D=1.5和D=1.9时变化幅度变小，这和图 5刚度的变化情况相对应；图 9(b)为低速时均方根值的放大图，从图中可以看出，当D=1.1时，在低速区出现多个峰值，系统相对不稳定，随着D的增大，系统均方根曲线变得平稳，在ω=0.3附近出现了峰值1、2，在ω=0.5附近也出现了峰值，与文献[19]中情况类似，符合齿轮的运动规律.

2.2 含随机扰动的时变刚度对系统响应的影响

为了验证本文提出的考虑齿面分形特性的时变刚度模型对系统响应分析的敏感性和优点，与含服从正态分布的随机扰动分量的时变啮合刚度模型进行对比.含随机扰动分量的刚度模型表达式为^[8]

$ {k_{\rm{h}}} = 1 - \left( {{\varepsilon _{\rm{r}}} + {\omega _\Delta }} \right)\cos \left( {{\mathit{\Omega} _{\rm{h}}}\tau } \right). $

(7)

式中：ε_r为量纲一齿轮啮合刚度、ω_△为服从高斯分布的量纲一啮合刚度的随机扰动量，Ω_h=ω_h/ω_n.

将由式(6)得到的不同分形维数下的刚度值与正常刚度值的差值作为刚度波动数组，该数组的均值和方差如表 4所示，代入式(7)，得到含随机扰动的时变啮合刚度，进而进行系统响应的分析，相关结果如图 10和图 11所示.

从图 10中的时域图可以看出，系统在不同随机扰动下的位移变化差别不大，但从表 4和图 11中看出，随着均值和方差的变化，系统振动的位移减小；与上节进行对比可知，随着分形维数D的减小，刚度波动降低，计算得到的方差变小，在服从正态分布的随机扰动下的系统振动位移减小，系统趋于稳定，但系统对于考虑齿面分形特性的时变啮合刚度的变化更加敏感，体现出由于齿面变化导致系统响应变化中的更多细节.

2.3 阻尼比对系统响应的影响

为了研究阻尼比对齿轮-轴承系统动力学特性的影响，取参数D=1.5，k_i=1.2，阻尼比ζ分别为0.01、0.05和0.10，结果如图 12所示.

从图 12可以看出，在ζ=0.01时，对于ω < 1.2区域，系统一直保持着单周期运动，之后存在两个不稳定区域；随着ζ的增大，当ζ=0.05时，两个不稳定区域均趋于稳定，低速区域由周期运动与混沌运动的混合区域转换为较小区域的拟周期及混沌运动，高速区域整体收缩；当ζ=0.1时，系统混沌区域消失，变为单周期运动.

3 结论

1) 在考虑齿面分形特性的齿轮-轴承系统中，随着分形维数D的增大，齿轮轮齿表面更加精细，时变啮合刚度的波动降低，系统趋于相对稳定的周期运动.

2) 齿轮-轴承系统对于考虑齿面分形特性的时变啮合刚度的变化更加敏感，由分形维数D的变化导致的齿轮齿廓的变化对系统响应有较大的影响.

3) 随着阻尼比的增大，系统趋向稳定，当阻尼比增大到一定程度时，系统会变为单周期稳定运动.