齿轮齿条机构因其结构简单紧凑、传动精度高、支撑刚度高、运动换向能力强被广泛应用于各类军事和民用装备中. 充分了解齿轮副动态特性，对减少传动过程中的振动和噪声、延长轮齿系统的寿命至关重要，因此齿轮齿条传动过程中的非线性动态特性一直是国内外的研究热点^[1]. Kahraman等^[2]利用单自由度集中参数模型描述轮齿间的时变动态特性，并试验验证了其模型的准确性^[3]. 为了更准确地表征复杂的齿轮特性，多自由度的集中参数模型被研究和应用于不同工况的齿轮系统中^[4-5]. 集中参数模型虽然可以提供极高的计算效率，但是该模型基于确定的定义参数来获得系统的动态响应，使得计算精度在很大程度上受参数选择的影响. 时变啮合刚度作为关键参数是影响模型预测准确性的重要因素，准确、高效、方便地确定不同结构与工况的轮齿啮合刚度也是获得齿轮系统动态特性的关键^[6].

在实际工作环境中，啮合齿副间存在润滑作用，润滑状况和齿轮系统动力传动间存在强耦合关系^[7]. 为了研究齿轮的润滑机理及其对齿轮动态性能的影响，许多学者开展了基于润滑摩擦原理的建模理论研究. 黄立等^[8]将渐开线斜齿圆柱齿轮啮合等效为长椭圆接触，建立了齿轮弹流润滑计算的数学模型；针对平稳运行和启动工况，研究了油膜受瞬态挤压效应的影响. De la Cruz等^[9]提出新颖的多啮合传动摩擦动力学模型，用以研究减速器的噪声和振动响应. 该模型建立时考虑了瞬态条件下的热弹性流体动力润滑、轮系动力学和粗糙体的界面作用等因素的耦合作用. 菅光霄等^[10]基于六自由度集中参数齿轮摩擦动力模型，探究齿轮啮合过程中接触冲击现象对齿轮弹流润滑特性的影响. 该研究表明，接触冲击振动使油膜厚度和油膜压力产生较大的突变. Tao等^[11]建立三维线接触混合弹流动力润滑模型，表征在接触面积和形貌影响下大模数齿轮齿条的润滑行为，并基于此模型对齿面的形态进行优化设计，以期在稀油供应条件下获得更好的传动性能. 上述模型大多是在固定负载和转速情况下考虑啮合刚度变化引起的啮合齿接触负载变化，这种先验地分配啮合齿负载的方法存在武断性^[7]. 同时，大多数齿轮齿条机构须频繁地启停转向，这会引起齿轮齿条啮合处速度和负载大范围变化，对轮齿间油脂润滑性能有极大的影响^[12].

本研究基于渐开线标准直齿轮传动的啮合原理，结合齿轮齿条增程机构的传动特性，考虑齿轮啮合时的结构时变啮合刚度和瞬态热弹流润滑刚度的耦合效应影响，建立脂润滑齿轮齿条增速机构的动力学模型. 根据测定剪切稳定状态下的润滑脂流变参数和热力学参数，研究传动过程中啮合处脂润滑热成膜与摩擦特性以及脂润滑对齿轮齿条传动性能的影响.

如图1所示为齿轮齿条增程机构示意图. 图中， $ {x_{{\rm{p}}} } $、 $ {\dot x_{{\rm{p}}} } $、 $ {\ddot x_{{\rm{p}}} } $分别为齿轮平动位移、速度和加速度； $ {\theta _{{\rm{p}}} } $、 $ {\dot \theta _{{\rm{p}}} } $、 $ {\ddot \theta _{{\rm{p}}} } $分别为齿轮转动角位移、角速度、角加速度； $ {r_{{\rm{b}}} } $为齿轮分度圆半径； $ {m_{{\rm{r}}} } $、 $x_{{\rm{r}}} $、 ${\dot x_{\rm{r}}} $、 $ {\ddot x_{{\rm{r}}} } $分别为上侧齿条质量（包含输送负载）、平动位移、平动速度和平动加速度； $ {F_{{{\rm{g}}} 1,j}} $、 $ {F_{{{\rm{f}}} 1,j}} $分别为下侧啮合处的法向啮合力和切向摩擦力； $ {F_{{{\rm{g}}} 2,j}} = {F'_{{{\rm{g}}} 2,j}} $、 $ {F_{{{\rm{f}}} 2,j}} = {F'_{{{\rm{f}}} 2,j}} $分别表示上侧啮合处的法向啮合力和切向摩擦力； $ {o_1} $、 $ {o_2} $分别为下侧和上侧啮合线与齿轮基圆的切点， $ {\xi _{2,j}} $、 $ {\xi _{1,j}} $分别为上侧和下侧啮合点到对应切点的距离. 齿轮在执行装置提供的力 $ {F_{{\rm{T}}} } $作用下开始运动，并通过与下侧固定齿条、上侧运动齿条的啮合，实现齿轮的转动与上齿条的平动. 齿轮齿条动态啮合过程中存在3种运动状态：齿面啮合、齿背啮合和齿轮脱啮，为了简洁反映啮合时的受力情况，图1仅标注齿面啮合状态的受力. 当啮合状态为齿背啮合时，可以通过更改啮合作用线方向（ $ {o_1}{x_1} \to {o'_1}{x'_1} $， $ {o_2}{x_2} \to {o'_2}{x'_2} $）来更改受力情况；当为齿轮脱啮状态时，啮合处无作用力. 定义 $ \bar D $为齿侧间隙的一半； $ \delta $为动态传动误差，计算上侧齿副动态传动误差时， $\delta = ( x_{{\rm{p}}} - x_{{\rm{r}}} )\cos \alpha +{r_{{\rm{b}}} }{\theta _{{\rm{p}}} }$，计算下侧齿副动态传动误差时， $ \delta = {x_{{\rm{p}}} }\cos \alpha - {r_{{\rm{b}}} }{\theta _{{\rm{p}}} } $. 不同啮合状态的判断条件如下：1）当 $ \delta > \bar D $时，为齿面啮合状态；2）当 $ \delta < - \bar D $时，为齿背啮合状态；3）当 $ - \bar D \leqslant \delta \leqslant \bar D $时，为齿轮脱啮状态. 本研究关注齿轮齿条机构内部传动的动态性能，忽略上侧运动齿条运行过程中滑移轨道处的摩擦作用，得到齿轮齿条传动机构的运动方程为

(1) $ \left. \begin{array}{l} {m_{{\rm{r}}} }{{\ddot x}_{{\rm{r}}} } = \displaystyle \sum\limits_{j = 1}^{{N_{{\rm{U}}} }} {\left( {{{F}^{'}_{{{\rm{g}}} 2,j}}\cos \alpha - {{F}^{'}_{{f} 2,j}}\sin \alpha } \right)} , \\ {m_{{\rm{p}}} }{{\ddot x}_{{\rm{p}}} } = {F_{{\rm{T}}} } + \displaystyle \sum\limits_{j = 1}^{{N_{{\rm{L}}} }} {\left( { - {F_{{{\rm{g}}} 1,j}}\cos \alpha - {F_{{{\rm{f}}} 1,j}}\sin \alpha } \right)} + \\ \quad \quad \quad \displaystyle \sum\limits_{j = 1}^{{N_{{\rm{U}}} }} {\left( { - {F_{{{\rm{g}}} 2,j}}\cos \alpha +{F_{{{\rm{f}}} 2,j}}\sin \alpha } \right)} , \\ {I_{{\rm{p}}} }{{\ddot \theta }_{{\rm{p}}} } = \displaystyle \sum\limits_{j = 1}^{{N_{{\rm{U}}} }} {\left( {{F_{{{\rm{g}}} 2,j}}{r_{{\rm{b}}} } - {F_{{{\rm{f}}} 2,j}}{\xi _{2,j}}} \right)} - \\ \quad \quad \quad \displaystyle \sum\limits_{j = 1}^{{N_{{\rm{L}}} }} {\left( {{F_{{{\rm{g}}} 1,j}}{r_{{\rm{b}}} }+{F_{{{\rm{f}}} 1,j}}{\xi _{1,j}}} \right)} . \\ \end{array} \right\} $

式中： $ {m_{{\rm{p}}} } $为齿轮质量（包含推动杆）； $ {I_{{\rm{p}}} } $为齿轮转动惯量；α为两侧齿条的压力角， $ \alpha = {20^ \circ } $；j为第j对啮合的齿对数目；N_U、N_L分别为上侧和下侧啮合齿对的数目，N_U=1、2，N_L=1、2.

大多数已有研究基于预设的载荷分配因子在确定负载的情况下进行齿轮的弹流润滑分析，忽略了非稳态运动下脂膜刚度对总时变啮合刚度的影响. 本研究考虑齿轮啮合时的结构时变啮合刚度和热弹流润滑（thermal elastohydrodynamic lubrication，TEHL）的耦合效应，提出结构−脂膜耦合啮合刚度模型.

利用势能理论建立啮合结构刚度模型，由于重合度大于1，将多齿同时啮合时引起的邻近齿耦合柔性引入啮合刚度模型中，表达式为

(2) $ \left. \begin{array}{l} \begin{split} \frac{1}{{K_{{\rm{S}}} ^i}} =& \frac{1}{{K_{{\rm{B}}} ^i}}+\frac{1}{{K_{{\rm{C}}} ^i}} \text{；} \\& \frac{1}{{K_{{\rm{B}}} ^i}} = \sum\limits_{i = 1}^l {\left( {\frac{1}{{k_{{{\rm{t}}} ,i}^{{\rm{p}}} }}+\frac{1}{{k_{{{\rm{f}}} ,i}^{{\rm{p}}} }}+\frac{1}{{k_{{{\rm{t}}} ,i}^{{\rm{r}}} }}+\frac{1}{{k_{{{\rm{f}}} ,i}^{{\rm{r}}} }}} \right)} \text{，} \\& \frac{1}{{K_{{\rm{C}}} ^i}} = \sum\limits_{i = 1}^l {\sum\limits_{j = 1}^l {\left( {\frac{1}{{k_{{{\rm{f}}} ,ij}^{{\rm{p}}} }}+\frac{1}{{k_{{{\rm{f}}} ,ij}^{{\rm{r}}} }}} \right)} } ,{\text{ }}i \ne j \text{；} \\& \qquad \qquad \quad \quad \frac{1}{{{k_{{\rm{t}}} }}} = \frac{1}{{{k_{{\rm{b}}} }}}+\frac{1}{{{k_{{\rm{s}}} }}}+\frac{1}{{{k_{{\rm{a}}} }}} . \end{split} \end{array} \right\} $

式中： $ i $为啮合的齿对数目， $ i=1、2 $； $ {K_{{\rm{S}}} } $为啮合结构刚度； $ {K_{{\rm{B}}} } $为未考虑耦合效应的结构刚度； $ {K_{{\rm{C}}} } $为耦合刚度； $ l $为瞬时啮合的齿对数量， $ l=1、2 $； $ {\text{p}} $、 $ {\text{r}} $分别表示齿轮和齿条； $ {k_{{\rm{b}}} } $、 $ {k_{{\rm{s}}} } $和 $ {k_{{\rm{a}}} } $分别为弯曲、剪切和轴向压缩刚度，参考文献[6]计算得到； $ {k_{{{\rm{f}}} ,i}} $、 $ {k_{{{\rm{f}}} ,ij}} $分别为基体变形刚度和基体耦合变形刚度，计算式分别为

(3) $ \frac{1}{{k_{{{\rm{f}}} ,i}^{{\rm{p}}} }} = \frac{{{{\cos }^2}\,{\alpha _1}}}{{{{Eb}}}}\left( {L_1^{{\rm{p}}} {{\left( {\frac{{{u_1}}}{{{s_{{\rm{f}}} }}}} \right)}^2}+M_1^{{\rm{p}}} \left( {\frac{{{u_1}}}{{{s_{{\rm{f}}} }}}} \right)+{\rm{P}}_1^{{\rm{p}}} \left( {1+Q_1^{{\rm{p}}} {{\tan }^2}\,{\alpha _1}} \right)} \right) \text{，} $

(4) $ \begin{split} \dfrac{1}{{k_{{\rm{f}},21}^{\rm{p}}}} =& \dfrac{{\cos {\alpha _1}\cos {\alpha _2}}}{{Eb}}\Bigg( {L_2^{\rm{p}}\left( {\dfrac{{{u_1}{u_2}}}{{s_{\rm{f}}^2}}} \right) + \left( {M_2^{\rm{p}}\tan {\alpha _2} + P_2^{\rm{p}}} \right)\left( {\dfrac{{{u_1}}}{{{s_{\rm{f}}}}}} \right) + } \\& \left( {Q_2^{\rm{p}}\tan {\alpha _1} + R_2^{\rm{p}}} \right)\left( {\dfrac{{{u_2}}}{{{s_{\rm{f}}}}}} \right) + \left( {S_2^{\rm{p}}\tan {\alpha _1} + } \right.\left. {T_2^{\rm{p}}} \right)\tan {\alpha _2} + \\& U_2^{\rm{p}}\tan {\alpha _1} + V_2^{\rm{p}} \Bigg) , \\[-15pt] \end{split}$

(5) $ \begin{split} \frac{1}{{k_{{\rm{f}},12}^{\rm{p}}}} =& \frac{{\cos {\alpha _1}\cos {\alpha _2}}}{{Eb}} \Bigg( {L_3^{\rm{p}}\left( {\frac{{{u_1}{u_2}}}{{s_{\rm{f}}^2}}} \right) + \left( {M_3^{\rm{p}}\tan {\alpha _1} + P_3^{\rm{p}}} \right)\left( {\frac{{{u_2}}}{{{s_{\rm{f}}}}}} \right) + } \\& \left( {Q_3^{\rm{p}}\tan {\alpha _2} + R_3^{\rm{p}}} \right)\left( {\frac{{{u_1}}}{{{s_{\rm{f}}}}}} \right) + \left( {S_3^{\rm{p}}\tan {\alpha _2} + \left. {T_3^{\rm{p}}} \right)\tan {\alpha _2} + } \right.\\& U_3^{\rm{p}}\tan {\alpha _1} + V_3^{\rm{p}} \Bigg) \text{，} \\[-15pt] \end{split}$

(6) $ \begin{split} \frac{1}{{k_{{{\rm{f}}} ,i}^{{\rm{r}}} }} =& \frac{{{{\cos }^2}{\alpha _i}}}{{Eb}}\Bigg( {L_1^{{\rm{r}}} {{\left( {\frac{{{u_i}}}{{\rm{a}}}} \right)}^2}+M_1^{{\rm{r}}} \left( {\frac{{{u_i}}}{{\rm{a}}}} \right)+LH_1^{{\rm{r}}} \left( {\frac{{u_i^2}}{{a{H_0}}}} \right)+ } \\ &{ MH_1^{{\rm{r}}} \left( {\frac{{{u_i}}}{{{H_0}}}} \right)+P_1^{{\rm{r}}} {{\tan }^2}{\alpha _i}+Q_1^{{\rm{r}}} } \Bigg) \text{，} \end{split}$

(7) $ \begin{split} \frac{1}{{k_{{\rm{f}},21}^{\rm{r}}}} = &\frac{{\cos {\alpha _1}\cos {\alpha _2}}}{{Eb}}\left( {L_2^{\rm{r}}\left( {\frac{{{u_1}{u_2}}}{{{a^2}}}} \right) + M_2^{\rm{r}}\tan {\alpha _2}\left( {\frac{{{u_1}}}{a}} \right) + } \right.\\& Q_2^{\rm{r}}\tan {\alpha _1}\left( {\frac{{{u_2}}}{a}} \right) + MH_2^{\rm{r}}\left( {\frac{{{u_2}}}{{{H_0}}}} \right) + QH_2^{\rm{r}}\tan {\alpha _1}\left( {\frac{{{u_2}}}{{{H_0}}}} \right) + \\& \left( {S_2^{\rm{r}}\tan {\alpha _1} + } \right. \left. {T_2^{\rm{r}}} \right)\tan {\alpha _2} + U_2^{\rm{r}}\tan {\alpha _1} + V_2^{\rm{r}} \Bigg) ， \\[-15pt] \end{split} $

(8) $ \begin{split} \frac{1}{{k_{{\rm{f}},12}^{\rm{r}}}} =& \frac{{\cos {\alpha _1}\cos {\alpha _2}}}{{Eb}}\Bigg( {L_3^{\rm{r}}\left( {\frac{{{u_1}{u_2}}}{{{a^2}}}} \right) + M_3^{\rm{r}}\tan {\alpha _1}\left( {\frac{{{u_2}}}{{\rm{a}}}} \right) + }\\& Q_3^{\rm{r}}\tan {\alpha _2}\left( {\frac{{{u_1}}}{a}} \right) + MH_3^{\rm{r}}\left( {\frac{{{u_1}}}{{{H_0}}}} \right) + QH_3^{\rm{r}}\tan {\alpha _1}\left( {\frac{{{u_1}}}{{{H_0}}}} \right) +\\& \left( {S_3^{\rm{r}}\tan {\alpha _1} + } \right. \left. {T_3^{\rm{r}}} \right)\tan {\alpha _2} + U_3^{\rm{r}}\tan {\alpha _1} + V_3^{\rm{r}} \Bigg) . \\[-15pt] \end{split} $

式中： $ E $为弹性模量，a为齿根宽度的一半，b为齿宽；各结构刚度项对总刚度的贡献情况，如图2所示；常系数与齿轮齿条结构有关，在本研究中其值如表1所示.

如图3所示，双齿啮合时，根据变形协调条件^[6,13-14]，即同时接触的2个齿对的基节应相等，总啮合力 $ {F_{{{\rm{tm}}} }} = {F_1}+{F_2} $，其中F₁、F₂分别表示啮合齿对1和2上的法向啮合力. 由于啮合齿对的变形导致实际啮合的起始与终止位置不同于理论值，即出现线外啮合现象^[14-16]，引入仅由结构参数决定^[16]的无负载分离距离 $ {S_{\rm{a}}} $（接近过程）、 $ {S_{\rm{r}}} $（离去过程）来反映线外啮合现象. 考虑结构−脂膜耦合的动态传动误差（沿啮合线方向）表示为1）接近过程：

(9) $ \left. {\begin{array}{*{20}{l}} {\varDelta _{{{\rm{TE}}} }^1 = {{E_{{{\rm{e}}} 1}} - {h_{{{\rm{o}}} 1}}} +{{{F_1}}}/{{K_{\rm{B}}^1}}+{{{F_2}}}/{{K_{\rm{C}}^1}}{\text{ }}}, \\ {\varDelta _{{{\rm{TE}}} }^2 = {{E_{{{\rm{e}}} 2}} - {h_{{{\rm{o}}} 2}}} +{{{F_1}}}/{{K_{\rm{C}}^2}}+{{{F_2}}}/{{K_{\rm{B}}^2}}+{S_{\rm{a}}}} . \end{array}} \right\} $

2）离去过程：

(10) $ \left. {\begin{array}{*{20}{l}} {\varDelta _{{\rm{TE}}}^1 = {{E_1} - {h_{{{\rm{o}}} 1}}} +{{{F_1}}}/{{K_{\rm{B}}^1}}+{{{F_2}}}/{{K_{\rm{C}}^1}}+{S_{\rm{r}}}} ,\\ {\varDelta _{{\rm{TE}}}^2 = {{E_2} - {h_{{{\rm{o}}} 2}}} +{{{F_1}}}/{{K_{\rm{C}}^2}}+{{{F_2}}}/{{K_{\rm{B}}^2}}{\text{ }}} . \end{array}} \right\} $

式中： $ \varDelta _{{{\rm{TE}}} }^1 = \varDelta _{{{\rm{TE}}} }^2 = \delta $为动态传动误差； $ {h_{{{\rm{o}}} 1}} $、 $ {h_{{{\rm{o}}} 2}} $分别为啮合齿对1和2处的刚体中心脂膜厚度； $ {E_{{{\rm{e}}} 1}} = {e_{{{\rm{g}}} 1}}+{e_{{{\rm{r}}} 1}} $、 $ {E_{{{\rm{e}}} 2}} = {e_{{{\rm{g}}} 2}}+{e_{{{\rm{r}}} 2}} $，其中 $ {e_{{{\rm{g}}} 1}} $、 $ {e_{{{\rm{g}}} 2}} $分别为啮合齿对1和2处的齿轮齿形误差， $ {e_{{{\rm{r}}} 1}} $、 $ {e_{{{\rm{r}}} 2}} $分别为啮合齿对1和2处的齿条齿形误差，本研究只考虑齿侧误差对齿形误差的影响，故有 $ {E_{{{\rm{e}}} 1}} = {E_{{{\rm{e}}} 2}} $. 当啮合齿对处于理论双齿啮合区域时， $ {S_{{\rm{a}}} } $、 $ {S_{{\rm{r}}} } $=0. 当 $ {S_{{\rm{r}}} } \geqslant \delta - \left( {{E_{{{\rm{e}}} 1}} - {h_{{{\rm{o}}} 1}}} \right) $、 $ {S_{{\rm{a}}} } \geqslant \delta - \left( {{E_{{{\rm{e}}} 2}} - {h_{{{\rm{o}}} 2}}} \right) $同时满足时，相应接触位置实际上处于单齿啮合状态，动态误差表示为

(11) $ \varDelta _{{{\rm{TE}}} }^1 = \left( {{E_1} - {h_{{{\rm{o}}} 1}}} \right)+{{{F_1}}}/{{K_{{\rm{B}}} ^1}} . $

根据渐开线直齿轮的啮合特点，将接触齿间的流体润滑视为线接触进行建模，如图4所示，本研究考虑在剪切稳定后的润滑情况，因此忽略两侧的微量端泄效应. 结合润滑脂Power-law流变模型，推导瞬态啮合齿对的线接触Relynolds方程为

(12) $\left. \begin{array}{l} \begin{split} & \frac{{\partial \rho {G_{{\rm{R}}} }}}{{\partial x}} - {v_{{\rm{m}}} }\frac{\partial }{{\partial x}}\left( {\rho h} \right) - \frac{\partial }{{\partial t}}\left( {\rho h} \right) = 0 \text{，} \\& {G}_{\mathrm{R}}={0.5}^{\tfrac{n+1}{n}} \frac{n}{2n+1} {h}^{\tfrac{2n+1}{n}} {\left(\frac{1}{\phi }\frac{\partial p}{\partial x}\right)}^{\tfrac{1}{n}} . \end{split} \end{array} \right\} $

式中： $\; \rho $为润滑脂密度； $ p $为脂膜压应力； $ \varphi $为塑性黏度； $ n $为塑性指数； $ {v_{{\rm{m}}} } $为润滑脂卷吸速度，计算式为1）上侧齿：

(13) $ {v_{{\rm{m}}} } = \frac{{{u_{{\rm{g}}} }+{u_{{\rm{r}}} }}}{2} = \frac{{{{\dot \theta }_{{\rm{p}}} }{\xi _{1,i}}+{{\dot x}_{{\rm{p}}} }\sin \alpha }}{2} . $

2）下侧齿：

(14) $ {\text{ }}{v_{{\rm{m}}} } = \frac{{{{\dot \theta }_{{\rm{p}}} }{\xi _{2i}}+({{\dot x}_{{\rm{r}}} } - {{\dot x}_{{\rm{p}}} })\sin \alpha }}{2} . $

$ h $为局部脂膜厚度(或间隙)，是与时间相关的函数，表达式为

(15) $ h(x,t) = {h_{{\rm{o}}} }\left( t \right)+\frac{{{x^2}}}{{2{R_x}\left( t \right)}}+v(x,t) . $

式中： $ {R_x}\left( t \right) $为当量曲率半径； $ E' $为当量弹性模型； $ v $为表面弹性变形，计算式为

(16) $ v(x,t) = - \frac{2}{{{{{\text{π}}}} E'}}\int_{{x_{{{\rm{in}}} }}}^{{x_{{{\rm{out}}} }}} {p(s,t)} \ln \; {(s - x)^2}{\rm{d}}s . $

润滑脂密度 $ \;\rho $、塑性黏度 $ \varphi $对压力和温度敏感，在此假设 $ \;\rho $、 $ \varphi $与温度和压力的关系式^[17-19]. 严格来说，塑性指数 $ n $也是压力和温度的函数，但是通常将其简化成常数处理^[12,18].

(17) $ \rho = {\rho _0}\left( {1+\frac{{0.6 \times {{10}^{ - 9}}p}}{{1+1.7 \times {{10}^{ - 9}}p}} - {\rho _T}\left( {T - {T_0}} \right)} \right) \text{，} $

(18) $ \begin{split} \varphi = &{\varphi _0}\exp \left( {\left( {\ln \; {\varphi _0}+9.67} \right)} \right. \times \\ & \left. {\left( {{{\left( {1+5.1 \times {{10}^{ - 9}}p} \right)}^{{z_0}}}{{\left( {\frac{{T - 138}}{{{T_0} - 138}}} \right)}^{{s_0}}} - 1} \right)} \right) . \end{split} $

式中： $ \;{\rho _0} $为环境温度和压力下的密度， $ \;{\rho _T} $为热膨胀系数， $ {T_0} $为环境温度， $ T $为润滑脂温度， $ {\varphi _0} $为环境温度和压力下的塑性黏度，取 $ {z_0} $=0.68、 $ {s_0} $=−1.1.

不计热辐射和体积力的影响，忽略 $ z $方向的热对流和 $ x $方向的热传导，非稳态弹流接触中润滑膜的能量方程为^[12]

(19) $ \begin{gathered} \rho {c_{{\rm{R}}} }\left( {\frac{{\partial T}}{{\partial t}}+u\frac{{\partial T}}{{\partial x}}+w\frac{{\partial T}}{{\partial x}}} \right) = {k_{R} }\frac{{{\partial ^2}T}}{{\partial {z^2}}} - \\ \qquad {\text{ }}\frac{T}{\rho }\frac{{\partial \rho }}{{\partial T}}\left( {\frac{{\partial p}}{{\partial t}}+u\frac{{\partial p}}{{\partial x}}} \right)+\varphi {\left( {\frac{{\partial u}}{{\partial z}}} \right)^{1+1/n}} . \\ \end{gathered} $

式中： $ {c_{{\rm{R}}} } $、 $ {k_{{\rm{R}}} } $为润滑脂等压比热容和导热系数； $ u $、 $ w $分别为润滑脂在 $ x $和 $ z $方向的流速，其中 $ u $的计算式为

(20) $ \begin{gathered} u = \left( {{{ \displaystyle \int_0^z {\left( {\dfrac{z}{\varphi }} \right)} }^{1/n}}{\rm{d}}z - \dfrac{{{{ \displaystyle \int_0^h {\left( {\dfrac{z}{\varphi }} \right)} }^{1/n}}{\rm{d}}z}}{{{{ \displaystyle \int_0^h {\left( {\dfrac{1}{\varphi }} \right)} }^{1/n}}{\rm{d}}z}}{{ \displaystyle \displaystyle \int_0^z {\left( {\dfrac{1}{\varphi }} \right)} }^{1/n}}{\rm{d}}z} \right) \times \\ {\left( {\dfrac{{{\rm{d}}p}}{{{\rm{d}}x}}} \right)^{1/n}}+\dfrac{{{u_{{\rm{g}}} } - {u_{{\rm{r}}} }}}{{{{ \displaystyle \int_0^h {\left( {\dfrac{1}{\varphi }} \right)} }^{1/n}}{\rm{d}}z}}{ \displaystyle \int_0^z {\left( {\dfrac{1}{\varphi }} \right)} ^{1/n}}{\rm{d}}z+{u_{{\rm{r}}} } . \\ \end{gathered} $

$ w $由 $ u $代入流量连续性条件 $ {{\partial \left( {\rho u} \right)}}/{{\partial x}}+{{\partial \left( {\rho w} \right)}}/{{\partial z}} = 0 $求得.

能量方程的边界条件，即上下界面温度的计算式^[12]为

(21) $ \left. \begin{array}{l} {{T_{{\rm{g}}} } = \dfrac{{{k_{R} }}}{{\sqrt {{\text{π}} {\rho _{{\rm{g}}} }{c_{{\rm{g}}} }{k_{{\rm{g}}} }{u_{{\rm{g}}} }} }} \displaystyle \int_{ - \infty }^x {{{\left. {\dfrac{{\partial T}}{{\partial z}}} \right|}_{z = h}}} \dfrac{{{\rm{d}}s}}{{\sqrt {x - s} }}+{T_0}} ,\\ {{T_{{\rm{r}}} } = \dfrac{{{k_{R} }}}{{\sqrt {{\text{π}} {\rho _{{\rm{r}}} }{c_{{\rm{r}}} }{k_{{\rm{r}}} }{u_{{\rm{r}}} }} }} \displaystyle \int_{ - \infty }^x {{{\left. {\dfrac{{\partial T}}{{\partial z}}} \right|}_{z = 0}}} \dfrac{{{\rm{d}}s}}{{\sqrt {x - s} }}+{T_0}{\text{ }}} . \end{array} \right\} $

式中： $ {\;\rho _{{\rm{g}}} } $， $ {c_{{\rm{g}}} } $和 $ {k_{{\rm{g}}} } $分别为齿轮的密度、等压比热容和导热系数， $ {\;\rho _{{\rm{r}}} } $， $ {c_{{\rm{r}}} } $和 $ {k_{{\rm{r}}} } $分别为齿条的密度、等压比热容和导热系数.

在工程实际中，齿轮齿条的表面不可能是完全平滑的，且在重载与频繁启停的工况下，接触表面间形成的润滑油膜很薄，不能完全分离2个粗糙的表面，此时不可避免地存在微观粗糙峰间的接触. 因此将啮合处的摩擦力 $ {F_{{\rm{f}}} } $定为由2个部分组成^[1,20]：1）粗糙表面直接接触产生的边界摩擦 $ {f_{{{\rm{asp}}} }} $，2）脂膜黏性力 $ {f_{{\rm{v}}} } $. 这2个部分均与传动过程中的脂膜厚度有关，计算式分别为

(22) $ \left. \begin{array}{l} \begin{split} {F_{{\rm{f}}} } = &{f_{{{\rm{asp}}} }}+{f_{{\rm{v}}} } ; \\& {f_{{\rm{v}}} } = {\int_{{x_{{{\rm{in}}} }}}^{{x_{{{\rm{out}}} }}} {\varphi \left. {{{\left( {\frac{{\partial u}}{{\partial z}}} \right)}^n}} \right|} _{z = 0/h}}{\rm{d}}x , \\& {f_{{{\rm{asp}}} }} = {k_{{{\rm{asp}}} }}\iint\limits_{{A_{{\rm{c}}} }} {{p_{{{\rm{asp}}} }}}{\rm{d}}{A_{{\rm{c}}} } . \end{split} \end{array} \right\} $

式中： $ {k_{{{\rm{asp}}} }} $为边界摩擦的摩擦系数， $ {k_{{{\rm{asp}}} }} $∈[0.08, 0.13]^[21]，取 $ {k_{{{\rm{asp}}} }} $=0.12； $ {A_{{\rm{c}}} } $为赫兹接触面积随着啮合位置和啮合力变化而变化的量； $ {p_{{{\rm{asp}}} }} $为微观粗糙峰接触压力，计算式^[22-23]为

(23) $ {p_{{{\rm{asp}}} }} = \frac{{8\sqrt 2 {\text{π}} }}{{15}}{\left( {\eta \gamma \sigma } \right)^2}E'\sqrt {\frac{\sigma }{\gamma }} . $

式中： $\sigma $为复合粗糙度， $ \sigma = {\left( {\sigma _{\rm{r}}^2 + \sigma _{\rm{g}}^2} \right)^{0.5}} $，其中 $ {\sigma _{{\rm{g}}} } $、 $ {\sigma _{{\rm{r}}} } $分别为齿轮和齿条的表面粗糙度，根据设计值，取 $ {\sigma _{\rm{g}}} = {\sigma _r} = 0.4{\text{ }} $ μm； $ \eta \gamma \sigma $、 $ {\sigma }/{\gamma } $均为与几何参数和表面粗糙度有关的参数，参考文献[11]，取 $ \eta \gamma \sigma = 0.04 $， $ {\sigma }/{\gamma } = 0.001 $. 当假设微观粗糙峰高度为高斯分布时，确定统计参数 ${F_{2.5}}\left( {{H_\text{σ}}} \right)$的计算式^[20]为

(24) $ {F_{2.5}}\left( {{H_\text{σ} }} \right) = \left\{ {\begin{array}{*{20}{l}} {A{{\left( {4 - {H_\text{σ} }} \right)}^Z},{\text{ }}{H_\text{σ} } \leqslant 4} ;\\ {0{\text{ , }} \qquad \qquad {H_\text{σ} } > 4} . \end{array}} \right. $

式中： ${H_\text{σ} }$为平均油膜厚度 $ \bar h $与复合粗糙度 $ \sigma $的比值；A、Z均为统计回归系数， $ A = 4.406\;8 \times {10^{ - 5}} $、 $ Z = 6.804 $^[20].

齿轮齿条的结构参数如表2所示，润滑脂的流变参数和齿轮齿条传动系统的热力学参数如表3所示. 润滑脂为多功能通用润滑脂，其主要流变参数由流变仪测定.

如图5所示为脂润滑齿轮齿条动力学模型求解的流程，其步骤如下.1）设置基本参数（包括齿轮结构参数与润滑脂参数）以及齿轮齿条初始状态. 2）将根据当前时刻的齿轮齿条状态得到的脂卷吸速度 $ {v_{{\rm{m}}} } $、当量曲率半径 $ {R_x}\left( t \right) $和动态误差 $ \delta $代入脂润滑求解程序中. 3）依次迭代求解式（12）和式（19），直至压力、温度和误差平衡项均达到收敛条件，输出脂润滑相关状态量. 4）根据脂润滑状态量计算啮合力和摩擦力. 5）获得当前时刻的式（1），利用4~5阶的Runge-Kutta-Fehlberg法求解动力学方程，得到下个时刻的齿轮齿条状态. 6）更新时间参量 $ t = t+\Delta t $，如果大于设定的仿真时长，则结束程序；否则返回步骤2）.

在脂润滑求解过程中，压力的求解采用5层多重网格法，最密层的节点数为961. 在每层压力修正过程中，Gauss-Seidel、Jacobi双极子迭代法被分别应用在低压和高压区域的压力迭代中，表面弹性变形的计算采用DC-FFT方法. 温度场仅在最密层进行求解，温度场 $ oz $方向的离散节点数为30. 在进行润滑计算时，载荷是待求项，因此引入误差平衡项 ${\varepsilon _\delta } = {{\left| {\varDelta _{{{\rm{TE}}} }^2 - \varDelta _{{{\rm{TE}}} }^1} \right|}} / ({{\varDelta _{{{\rm{TE}}} }^2+\varDelta _{{{\rm{TE}}} }^1}})$来替代载荷平衡条件. 压力、温度和误差平衡项的迭代误差分别小于 $ 1 \times {10^{ - 4}} $、 $ 1 \times {10^{ - 3}} $和 $ 1 \times {10^{ - 4}} $.

齿轮齿条增程机构的前向（图1中X轴的正方向为前向）输送过程为主工况，与反向收回过程相比，前向输送的运行条件更为苛刻，因此将向前输送过程中的动态特性作为主要研究对象. 在该工况下，执行装置提供的力规律为

(25) $ {F_{{\rm{T}}} } = \left\{ {\begin{array}{*{20}{l}} \quad {1\; 220{\text{ }}, \;\; {\text{ 0}} \leqslant t \leqslant {{27} \mathord{\left/ {\vphantom {{27} {340}}} \right. } {340}}} ;\\ { - 1\; 220{\text{ , }} \; \; {{27} \mathord{\left/ {\vphantom {{27} {340}}} \right. } {340}}{\text{ < }}t \leqslant {{54} \mathord{\left/ {\vphantom {{54} {340}}} \right. } {340}}} ;\\ \qquad \;\; {0{\text{ }}, \; \, {{{\text{ }}54} \mathord{\left/ {\vphantom {{{\text{ }}54} {340}}} \right. } {340}} < t \leqslant 0.16} . \end{array}} \right. $

在整个工作周期内，齿轮与齿条的运动状态不仅受到外负载的影响，还与啮合点处的结构刚度和脂膜刚度有关，齿轮与齿条的瞬时运动状态影响润滑脂的成膜特性，这种齿轮齿条动力特性与热弹流体动压润滑特性的耦合作用，使啮合处的受力情况更为复杂. 上、下齿条的运动状态不同，导致同时刻上、下齿条啮合处的成膜特性不同，但在每个齿的啮合周期内的成膜规律是相同的，且上齿条作为输出端，其动态特性更为研究者所关注，因此仅提取并分析上齿条与齿轮的啮合动态特性.

如图6所示为上齿条的位移S、速度v随时间变化的关系曲线. 可以看出，上齿条实际位移与理论位移存在偏差，说明在传动过程中有能量被耗散. 速度曲线整体上满足先均匀加速后均匀减速的规律，但是在较小的时间域里存在明显波动. 波动的存在表明，传动系统是高频的强非线性系统. 传动过程中啮合状态的转换、总啮合刚度随运动状态的变化都是导致速度波动的原因. 由图还可以发现，减速段的速度波动明显强于加速段. 说明在减速阶段，啮合刚度的变化更强烈，还可能存在更频繁的冲击效应.

选取2个啮合齿对，分别为第2、8啮合齿，讨论它们在各自啮合周期内的动态传动误差、动态啮合力、总啮合刚度与脂膜特性的变化. 第2啮合齿是初始加速阶段对应的啮合齿，第8啮合齿是传动从加速阶段切换到减速阶段过程对应的啮合齿.

如图7所示为第2啮合齿在啮合周期内的状态特性. 由图7（a）~（c）可知，啮合周期内仅包含2种啮合状态：齿面啮合、齿轮脱啮。随着啮合齿从啮入运行到啮出（接触点到节点的转角偏差 $ \Delta \theta $从正到负），动态误差速率和总啮合力的波动幅值下降，速率振荡幅值从0.5 m/s下降到0.2 m/s，总啮合力振荡幅值从4 500 N下降到3 000 N，系统中由初始啮入引起的振荡开始逐渐衰退. 由图7（d）~（e）可知，在考虑脂膜的刚度特性后，啮合齿的总刚度下降，且作用在啮合处的法向啮合力越小，总啮合刚度相较结构刚度下降得越大. 第2齿啮合过程是系统运行的加速时期，此时脂膜状态的规律基本一致，即当量半径由啮入时的较小值（1.48 mm）增加到啮出时的较大值（17.16 mm），中心压应力随当量曲率半径增加而下降，中心膜厚随当量曲率半径增加而上升，如图7（f）~（g）所示. 齿面啮合力是高频的动态变化，且齿轮齿条瞬时速度也是变化的，因此在某次齿面啮合过程中，中心压应力与中心脂膜是动态变化的. 选取点A、B和C分别代表啮入双齿接触阶段、单齿接触阶段和啮出双齿接触阶段中脂膜最薄的点，它们对应压力、膜厚与温度的分布关系如图7（h）~（k）所示. 当第2齿刚刚进入啮合时，当量曲率半径很小，点A的卷吸速度较其他点低，该点出现较大的压应力值，中心处脂膜压应力为0.948 GPa，二次压应力峰的高度也急剧增加，达到1.285 GPa，明显大于脂膜中心压应力. 此时的脂膜厚度为第2齿啮合周期中的最小值，在颈缩处的脂膜最小值为29 nm，最大温升出现在膜厚缩颈段的初始位置，为64.8 K. 当第2齿位于啮合点B、C时，当量曲率半径相比点A时的有明显增长，卷吸速度也有较大提升，因此压应力二次峰的峰值明显低于各自的中心压应力0.713 GPa和0.509 GPa，最小膜厚也分别增加到268 、334 nm. 由于负载压力的明显下降，点B的最大温升仅有7.0 K，点C位置的温度基本没有变化（此时相滑动速度接近于零）.

系统继续运行，进入第8齿啮合期，此周期内的状态特性如图8所示. 第8齿啮合的前半段时期处于上齿条加速阶段，后半段转入减速阶段，因此啮合齿对的状态特性比第2齿啮合时期的更为复杂. 由图8（a）~（c）可知，啮合周期内3种啮合状态均出现，齿背啮合的存在使得第8齿的啮合周期明显增加. 在该齿啮合周期初期，动态误差速率和总啮合力的波动幅值较为稳定且处于较低水平，其中动态误差速率的波动幅值为0.2 m/s，总啮合力振荡幅值为2 000 N. 当系统突然转入减速阶段出现齿背啮合状态时，会引起动态误差速率和总啮合力大幅振荡. 如图8（d）所示，不同于第2齿啮合时的双齿−单齿−双齿的3段刚度，此时的结构刚度变为双齿−单齿−双齿−单齿−双齿的5段刚度，且在−0.28< $ \Delta \theta < - 0.05 $的区域，齿背啮合和齿面啮合的状态交替出现，使得结构刚度值出现明显的波动. 同理，耦合脂膜的刚度项后，啮合齿总刚度下降，且啮合处作用的啮合力越小，总啮合刚度相较结构刚度下降得越大，因此总刚度呈现高频波动的特性，如图8（e）所示. 中心压应力变化趋势与当量曲率半径的变化相反，中心膜厚变化趋势与当量曲率半径的变化相同，如图8（f）~（g）所示. 在第8啮合周期内，最小的中心膜厚为162 nm，此时的中心压应力最大，为2.299 GPa. 在此选取6个点，其中点A、B和C分别为齿面啮合状态下啮入双齿接触阶段、单齿接触阶段和啮出双齿接触阶段中脂膜最薄的点，点D、E和F分别为齿背啮合状态下啮入双齿接触阶段、单齿接触阶段和啮出双齿接触阶段中脂膜最薄的点. 由图8（h）~（n）可知，在当量曲率半径低于4 mm的点A、F，脂膜润滑状态最恶劣，点A对应的最大压应力为1.651 GPa、最小膜厚为189 nm、最大温升为123.3 K；点F对应的最大压力为1.612 GPa、最小膜厚为162 nm、最大温升分121.8 K.

啮合齿对间脂膜特性的动态变化，使摩擦力呈现复杂的波动特性，啮合齿间的摩擦作用进一步影响系统的整体特性，导致出现如图6（a）所示的上齿条未达到理论位置的情况. 通过摩擦系数表征系统运行过程中摩擦特性的变化，定义摩擦系数为 $ \mu = {F_{{\rm{f}}} }/{F_{{\rm{g}}} } $. 如图9所示，在系统运行周期内啮合齿对的摩擦系数呈现高频的波动性，在齿轮齿条运行速度较大的中间时段，摩擦系数均小于0.08. 在齿轮齿条运行的加速阶段，摩擦系数在齿对啮合周期内呈现相似性的规律. 如图10所示，以第2齿啮合周期内的摩擦系数的变化为例，摩擦系数整体上呈现先上升后下降的趋势，在啮合点靠近节点位置的时，摩擦系数的波动幅值最大，原因是只有在接近节点位置时相对滑移速度才能趋近于零. 在齿轮齿条运行的后半段（减速阶段），轮齿在齿面、齿背啮合状态间切换难以找到其一般性规律.

选取总刚度模型为验证目标，将本研究所提出方法与有限元模型（finite element modelling，FEM）的计算结果进行对比. 齿轮齿条的有限元模型的基本参数与表2、3相同. 如图11所示，在有限元模型的齿条基座面上施加Z方向的强制平动位移边界条件，并约束其他方向上的自由度. 齿轮中心处施加Z方向的强制平动位移和X轴方向的强制角位移，并约束其他方向上的自由度. 在齿轮的间隙部位，使用欧拉格式的有限元背景网格离散润滑脂可能的运动区域. 对于润滑脂的初始几何形态，同样使用欧拉网格再次离散，并赋予初始体积分数为1. 对背景网格进行布尔运算，取其余区域的体积分数为0. 在接触区域对网格进行加密处理，并采用增广拉格朗日方法防止结构体互穿.

如图12所示为第2齿的啮合刚度对比. 可以看出，本研究所提结构−脂膜耦合啮合刚度模型与FEM的计算结果基本相同. 在单齿啮合到双齿啮合的过渡阶段，线外啮合的现象在2种模型的计算结果中都得到体现，该阶段2种计算结果的误差均较大，但是均未超过3.5%，其他啮合阶段的计算结果误差均未超过2%. 在FEM计算中得到与上文分析一致的规律，即考虑脂膜影响后的齿轮齿条啮合刚度低于仅考虑结构影响的，且啮合刚度随着啮合处作用力的下降而显著下降. 因此，本研究所提方法在计算脂润滑条件下的齿轮齿条啮合刚度时，不仅具有较高的计算精度和效率，而且能够有效地反映齿轮齿条啮合过程中的动态特性变化.

（1）当在齿轮啮合过程中考虑润滑脂的瞬态热弹流效应后，轮齿间的啮合总刚度将比仅考虑结构影响时的啮合刚度低，且啮合力越低，刚度值下降越明显，局部高频刚度最大波动幅值为1.1×10⁸ N/mm，最大下降幅度为30%. 不同于准静态情况^[7,11]，啮合力不再沿着啮合线规律分配；与准静态情况下不同，此时中心压应力与中心膜厚出现高频波动现象，中心压应力的最大波动幅值为1.6 GPa，中心膜厚的最大波动为300 nm.

（2）齿轮齿条传动系统在运行过程中存在3种工作状态：齿面啮合、齿背啮合、齿轮脱啮. 系统在前期加速运动时仅出现齿面啮合和脱啮的状态，此阶段系统呈现规律性运行：随着轮齿啮合位置沿啮合线正方向移动，中心脂膜厚度增加，中心压应力下降. 系统在后期减速运动时，由于齿面啮合与齿背啮合的交替出现，当量曲率半径将出现较大的阶跃变化，使得单齿啮合周期内总啮合刚度、中心膜厚、中心压应力和摩擦系数均出现较大幅度的阶跃，系统运行无规律可循.

（3）无论是齿背啮合状态还是齿面啮合状态，最小油膜厚度、最大温升以及最大压应力即系统最恶劣润滑状态，均出现在齿轮轮齿面上靠近基圆的位置及与其啮合的齿条齿顶附近. 当每对啮合齿啮合到靠近节点位置时，脂膜温升和摩擦系数明显下降.

（4）本研究建立的在脂润滑情况下齿轮齿条动力学模型较为真实地体现了实际传动过程中的动态特性，并很好地反映了各阶段的润滑情况。期望本研究所得结论能够为后续齿轮磨损性能的研究奠定基础.