斜齿圆柱齿轮载荷分布及热弹流温度场分析

Study on Load Distribution and Thermal Elastohydrodynamic Lubrication Temperature Fields of Helical Gear System

摘要: 啮合面载荷分布是斜齿轮设计和强度校核的基础。以传动误差为基础，基于啮合面和端面刚度建立了斜齿轮单位线载荷分析模型，通过数值解法得到了斜齿轮啮合面上单位线载荷和转角误差分布。为便于设计和校核，建立了既综合考虑齿廓因素和轴向因素，又能反映啮合面载荷的特征坐标系。将斜齿轮副简化为两个反向圆锥台接触模型，完善了斜齿轮热弹流分析模型，得到了斜齿轮接触点油膜压力、厚度和温度场分布，得到了沿特征坐标分布的闪温。结果表明：斜齿轮齿廓中部承担了大部分载荷，其变化规律与转角误差相同。本文建立的特征坐标系能够反映齿轮载荷分布特征和满足设计和校核的要求。高速重载齿轮系统热弹流解闪温分布与Blok闪温分布规律基本相同。本文可以为斜齿轮设计和胶合强度校核提供理论依据。

Abstract: The load distribution of meshing surface is the basis of design and strength rating of helical gear. The unit-line load analysis model based on the transmission error, meshing face and transverse stiffness of helical gear was established. The distribution of unit-line load and transmission error was calculated by numerical method. The feature coordinate system considering both tooth profile and axial factor was established for ease of design and checking. Two reverse tapers’ contact model was proposed to simulate the helical gear pair. The more accurate thermal elastohydrodynamic lubrication analysis model was established. The oil film thickness, pressure and temperature fields of contact point and flash temperature along the feature coordinate were obtained. The results show: the middle part of tooth profile bears a large proportion of load and its variation law is the same as the angle error. The feature coordinate system proposed in this paper can reflect the characteristics of the load distribution and satisfy the requirement of helical gear strength check. The flash temperature calculated by thermal elastohydrodynamic lubrication method is almost the same with Blok flash temperature in dedendum and middle part of tooth profile. The study in the paper can provide the basis for the helical design and check of scuffing load capacity of helical gear pair.

1. 引言

平行轴斜齿圆柱齿轮是高速重载传动中的首选，其重合度高，传动平稳，振动和噪音小，已得到广泛应用。斜齿轮啮合过程中单位线载荷的计算是斜齿轮设计和强度校核的基础，Miryam B等 [1] [2] 通过弹性势能法研究了斜齿轮单位线载荷沿接触线的分布规律，Yanjun Peng等 [3] 研究了修形和错位情况下的斜齿轮载荷分布，李文良等 [4] 研究了接触线长度及螺旋角对其的影响，但载荷计算较繁琐，较难用于工程计算。

热弹流分析是校核承载能力的重要方法，杨沛然等 [5] 研究了两个反向圆锥滚子的热弹流数值解，王优强等 [6] [7] 研究了斜齿轮瞬态热弹流润滑，但其研究重点在建模和求解方法，都是直接给定单位线载荷，无法将其直接应用于斜齿轮的设计与校核中。

本文将载荷和热弹流计算引入工程应用，通过转角误差基于啮合面和端面刚度得到了斜齿轮啮合面上每一点的单位线载荷分布，提出了简化斜齿轮设计和强度校核的特征坐标系，建立了考虑润滑油粘合和密度随温度和压力变化的斜齿轮热弹流分析模型，得到了啮合面上的热弹流温度场，为斜齿轮设计与校核奠定了理论基础。

2. 斜齿轮单位线载荷计算

2.1. 斜齿轮单位线载荷分析模型

斜齿轮在传动过程中由于受载变形会产生传动误差，如式(1)所示，式中 $\Delta \theta$ 为传动误差角， $r_{b2}$ 为从动轮基圆半径。则该点单位线载荷为 $w=k\delta$ 。设该时刻接触线总长度为L，不考虑基节误差，则必满足载荷平衡方程，如式(2)，P为输入功率n₁为小轮转速，r_b1为小轮基圆半径。

$\delta =r_{b2}\cdot \Delta \theta$ (1)

${\displaystyle {\int }_0^{L}w\text{d}l}={\displaystyle {\int }_0^{L}k}\left(l\right)r_{b2}\Delta \theta \text{d}l=\frac{9549P}{n_1{r}_{b1}}$ (2)

2.2. 斜齿轮单位线载荷计算

斜齿轮端面刚度可以由石川模型求得，斜齿轮参数如表1所示，斜齿轮端面刚度分布如图1所示。

Table 1. The parameters of the helical gear

表1. 斜齿轮参数

Figure 1. The stiffness distribution of helical gear end face

图1. 斜齿轮端面刚度分布

斜齿轮啮合面 $N_1{N}_2{N}_3{N}_4$ 为其基圆内公切面，如图2所示，接触线 $K_1{K}_2$ 为与轴线成基圆螺旋角 ${\beta }_b$ 的线段，根据端面啮合状态，可以将啮合区分为双齿啮合区 $A_1{A}_2{B}_1{B}_2$ 、 $D_1{D}_2{E}_1{E}_2$ 和单齿啮合区 $B_1{B}_2{D}_1{D}_2$ ，在双齿啮合区内，接触线总是对应成对出现的。设齿轮副从 $A_1$ 点啮入，点啮出，斜齿轮轴向重合度 ${\epsilon }_{\beta }=b/P_{ba}=b\tan {\beta }_b/\left(\pi m_t\cos {\alpha }_t\right)$ ，设 ${\epsilon }_{\beta }=n+\Delta {\epsilon }_{\beta }$ ，n为整数。 $n{P}_{ba}$ 对应的区域内，接触线总长度为常数，

由式(3)计算，其总载荷为 $W_n={\displaystyle {\int }_0^{L_n}k\left(l\right)r_{b2}\Delta \theta \text{d}l}$ ，将端面实际接触长度 $A_1{E}_1$ 离散化，设 $l_{\Delta }=A_1{E}_1/\left(N-1\right)$ ，

N为划分节点数，则 $n{P}_{ba}$ 对应载荷如式(4)。

Figure 2. The meshing model and contact face of helical gear

图2. 斜齿轮啮合模型及啮合面

$L_n=n{\epsilon }_{\alpha }{P}_{bt}/\sin {\beta }_b$ (3)

$W_n=n{\displaystyle \underset{i=1}{\overset{N}{\sum }}{k}_i{r}_{b2}{l}_{\Delta }\Delta \theta /\sin {\beta }_b}$ (4)

$\Delta b=\Delta {\epsilon }_{\beta }{P}_{ba}$ 对应的接触线长度和分布规律是随时间变化的，与 $\Delta B$ 和 ${\beta }_b$ 大小有关，可分四种情况，如图3所示，在s1区域，接触线长度最长，s2区域最短，其余区域是渐变的，啮合线长度 $L_{\Delta }$ 为该时刻各段接触线长度之和， $W_{\Delta }$ 为各点单位线载荷之和，设某时刻共有m条线段参与啮合，每条线段上离散后点数为 $num$ ，则有：

$W_{\Delta }={\displaystyle \underset{i=1}{\overset{m}{\sum }}{\displaystyle \underset{j=1}{\overset{num-1}{\sum }}{k}_j{r}_{b2}{l}_{\Delta }\Delta \theta /\sin {\beta }_b}}$ (5)

$L=L_n+L_{\Delta }$ (6)

$W=W_n+W_{\Delta }$ (7)

Figure 3. The analysis model of length of the contact line of helical gear

图3. 斜齿轮单位线载荷分析模型

斜齿轮接触线总长度如式(6)，总载荷如式(7)。采用数值计算方法，取初值 $L=L_n$ 和 $k_0$ ，求得 $\Delta {\theta }_0$ ，由式(7)得到总载荷 $W_0$ ，根据总载荷与实际载荷差值修正 $\Delta \theta$ ，直到载荷精度满足要求，即得到该时刻的转角误差 $\Delta \theta$ ，根据 $w=k\delta$ 得到接触线上该点的单位线载荷。其啮合面单位线载荷如图4所示，啮合过程中转角误差分布如图5所示，斜齿轮单位线载荷分布趋势与转角误差规律相同，呈方波状分布，方向与齿宽方向成 ${\beta }_b$ 角度，只是由于各点啮合刚度不同，齿廓中部单位线载荷比较大，齿顶和齿根单位线载荷较小。

Figure 4. The unite-linear load distribution of helical gear

图4. 斜齿轮三维单位线载荷分布

Figure 5. The 3D transmission error distribution of helical gear

图5. 斜齿轮三维传动误差分布

斜齿轮啮合面上每一点单位线载荷、滑动速度、曲率半径等不同，其载荷分布、应力分布、温度分布均应为三维模型，但三维模型计算量大，设计或校核运算周期长，而实际设计和强度校核过程中，往往只关心危险点的强度，故为了既能承载斜齿轮相关信息，又能体现斜齿轮本质，建立斜齿轮特征坐标系，用线段 $\overline{A_1{E}_2}$ 上各点来表征整个啮合面信息，坐标采用各点齿廓坐标，为示与端面坐标区别，用无量纲坐标 $\Psi$ 表示，定义与 $\Gamma$ 相同，特征坐标经过啮入点和啮出点，啮入点和啮出点接触线长度最短，滑动速度大，曲率半径小，故为斜齿轮最危险点之一，线段 $\overline{A_1{E}_2}$ 既包含了齿廓方向信息，也包括了齿宽方向信息。特征坐标下单位线载荷和转角误差如图6所示，显然，啮入点载荷和对应轴向最大值基本相同，转角误差也相同，因此特征坐标能反映整个啮合面的变化情况。

Figure 6. Unite-linear load and transmission error distribution of helical gear

图6. 斜齿轮单位线载荷和传动误差分布

3. 斜齿轮热弹流计算

3.1. 斜齿轮分析基本模型

齿轮系统在运转过程中会在接触处形成润滑油膜，由于挤压和粘性剪切，油膜内能增加，温度升高，通过对流换热和热传导使齿轮温度升高。热稳态时，轮齿本体温度高于润滑油温度，因此在进入啮合瞬间，轮齿作为热源对润滑油加热，使润滑油温度升高，由于油膜极薄，几乎接触瞬时润滑油温度即与本体温度一样，以往研究均以热稳态润滑油温度为入口油温，这在接触区前半段，润滑油界面温度低于轮齿温度，这与实际不符，本文以热平衡时本体温度作为初始温度 $t_0$ ，使热弹流温度更符合工程实际。

3.1.1. 斜齿轮热弹流基本方程

斜齿轮接触线上曲率半径沿齿宽是逐渐变化的，接触线上K点到 $N_1{N}_2$ 距离即为主动轮曲率半径 $r_{K1}$ ，到 $N_3{N}_4$ 距离为从动轮曲率半径 $r_{K2}$ ，可简化为两个个反向圆锥台相对滚动，如图7所示，在xoz平面内为两个椭圆，其长短轴如式(8)，在接触处曲率半径如式(9)。

${a_{Ki}|}_{i=1,2}=\frac{r_{Ki}{\cos }^2{\beta }_b}{1-2{\sin }^2{\beta }_b}$ ， ${b_{Ki}|}_{i=1,2}=\frac{r_{Ki}\cos {\beta }_b}{1-2{\sin }^2{\beta }_b}$ (8)

${R_i|}_{i=1,2}=\frac{a_{Ki}^2}{b_{Ki}}=\frac{r_{Ki}{\cos }^3{\beta }_b}{1-2{\sin }^2{\beta }_b}$ (9)

${u_i|}_{i=1,2}=\pi \cdot n_i{R}_i/30$ (10)

3.1.2. 热弹流基本方程

流体润滑雷诺方程为：

Figure 7. Contact analysis model of helical gear

图7. 斜齿轮齿面接触分析模型

$\frac{\partial }{\partial x}\left(F_{\text{2}}\frac{\partial p}{\partial x}\right)=u_2\frac{\partial \left(\rho h\right)}{\partial x}-\left(u_2-u_1\right)\frac{\partial }{\partial x}\left(\frac{\rho F_1}{F_0}\right)$ (11)

式中： $F_0={\displaystyle {\int }_0^{h}1/\eta \text{d}z}$ ， $F_1={\displaystyle {\int }_0^{h}z/\eta \text{d}z}$ ， $F_2={\displaystyle {\int }_0^h\rho \cdot z/\eta \left(z-F_1/F_0\right)\text{d}z}$ ， $\rho ,\eta$ 分别为考虑压力和温度影响的润滑油密度和粘度，如式(12)、(13)，h为油膜厚度，如式(14)，其中p必须满足载荷平衡方程，如式(15)。

$\rho \left(p,t\right)={\rho }_0\left(1+\frac{0.6\times {10}^{-9}p}{1+1.7\times {10}^{-9}p}-0.00065\left(t-t_0\right)\right)$ (12)

$\eta \left(p,t\right)={\eta }_0\exp \left\{A_1\left[-1+{\left(1+A_{2}p\right)}^z{\left(\frac{t-138}{t_0-138}\right)}^{-s_0}\right]\right\}$ (13)

式中： $A_1=\ln {\eta }_0+9.67$ ， $A_2=5.1\times {10}^{\text{-}9}$ ， $z=\alpha /\left(A_1{A}_2\right)$ ， $s_0=\beta /\left(A_1\left(t_0-138\right)\right)$ 。

$h=h_0+\frac{x^2}{2R}-\frac{4}{\pi E}{\displaystyle {\int }_{x_a}^{x_b}p\left(x\right)\ln \left|x-s\right|\text{d}s}$ (14)

式中， $h_0$ 为油膜中心处膜厚，R为综合曲率半径，E为当量弹性模量。

$R=\frac{R_1{R}_2}{R_1+R_2}$ ， $\frac{1}{E}=\frac{1}{2}\left(\frac{1-{\nu }_1^2}{E_1}+\frac{1-{\nu }_2^2}{E_2}\right)$

$w={\displaystyle {\int }_{x_a}^{x_b}p\left(x\right)\text{d}x}$ (15)

能量方程：

$\stackrel{对流能量}{\stackrel{︷}{\rho \cdot c_p\left(u\frac{\partial t}{\partial x}\right)}}-\stackrel{传导能量}{\stackrel{︷}{\lambda \frac{{\partial }^{2}t}{\partial z^2}}}+\stackrel{压缩发热}{\stackrel{︷}{\frac{t}{\rho }\frac{\partial \rho }{\partial t}u\left(\frac{\partial p}{\partial x}\right)}}=\stackrel{粘性发热}{\stackrel{︷}{\eta {\left(\frac{\partial u}{\partial z}\right)}^2}}$ (16)

运动方程：

$\frac{\partial p}{\partial x}=\frac{\partial }{\partial z}\left(\eta \frac{\partial u}{\partial z}\right)$ (17)

能量方程边界条件：

$\{\begin{array}{l}t\left(x,0\right)=\frac{\lambda }{\sqrt{\pi {\rho }_1{c}_1{u}_1{k}_1}}{\displaystyle {\int }_{-\infty }^x{\frac{\partial t}{\partial z}|}_{z=0}\frac{\text{d}s}{\sqrt{x-s}}+t_a}\\ t\left(x,h\right)=\frac{\lambda }{\sqrt{\pi {\rho }_2{c}_2{u}_2{k}_2}}{\displaystyle {\int }_{-\infty }^x{\frac{\partial t}{\partial z}|}_{z=h}\frac{\text{d}s}{\sqrt{x-s}}+t_b}\end{array}$ (18)

3.2. 热弹流方程的求解及结果分析

3.2.1. 热弹流方程的求解

将方程(11)~(18)进行无量纲化，采用有限差分法进行离散，x方向离散成81个不等距点，z方向离散成21个等距点，在给定温度场情况下，求解雷诺方程，得到油膜压力和厚度分布，在给定压力和油膜厚度情况下求解能量方程，采用逐行扫描法求解温度场，将两者联立循环求解，直至得到稳定温度场。采用有限元方法求解本体温度，热平衡时润滑油温度 $t_{oil}=60˚\text{C}$ ，通过有限元法 [8] [9] 求解得到本体温度 $t_0=396\text{K}\left(123˚\text{C}\right)$ ，如图8，润滑油参数如表2所示。

Figure 8. The bulk temperature fields solved by FEM method

图8. 斜齿轮有限元本体温度场分布

Table 2. The parameters of the oil

表2. 润滑油参数

3.2.2. 斜齿轮热弹流温度场计算结果

图9给出了啮入点的油膜厚度和压力分布，对于高速重载齿轮，油膜中心压力往往比较大，本例中中心压力已经达到1.45 GPa，油膜厚度约为0.4 μm，与齿面粗糙度数量级相当。由于压力较大，故压力分布和油膜分布不再具有典型的弹流特征，二次压力峰和颈缩靠近出口且很小，压力分布接近赫兹压力分布。

接触区的三维温度场如图10所示，(a)为三维温度场分布，(b)为中心油层和界面温度分布，可见油膜中层温度远大于界面温度，油膜中层温度分布与压力分布趋势相同，在压力最大处取得最大值，而界面温度则略有不同，温度一直升高在出口附近略有下降，取其最大值作为最高温度，其与本体温度之差为该点闪温。

Figure 9. The thickness and pressure of the oil film

图9. 斜齿轮油膜厚度和压力分布

(a) 三维温度场分布 (b) 油层中心及界面温度分布

Figure 10. The temperature fields of the oil film

图10. 斜齿轮油膜温度场分布

3.2.3. 热弹流闪温与Blok闪温的比较

闪温是引起齿轮热胶合破坏的主要原因，齿面高温致使油膜破裂，两齿面直接接触，摩擦系数增大，从而导致更高温度，使齿面撕裂，引起齿轮热胶合破坏。Blok根据摩擦理论推导了齿面接触闪温公式，如式(19)。

$t_f=0.7858\frac{fw\left|u_1-u_2\right|}{\left(\sqrt{{\lambda }_1{\rho }_1{c}_1{u}_1}+\sqrt{{\lambda }_2{\rho }_2{c}_2{u}_2}\right)\sqrt{b_0}}$ (19)

式中， $f=0.06$ 为摩擦系数，w为单位线载荷， ${\lambda }_i,{\rho }_i,c_i,u_i$ 为轮齿热传导系数、密度、比热容和切向滑动速度，其中1、2分别表示主从动轮， $b_0$ 为赫兹接触半宽。

图11为特征坐标下的热弹流闪温和Blok闪温的比较，两者分布趋势基本相同，在主动轮齿根和齿廓节点附件，热弹流闪温大于Blok闪温，而在齿顶附近，热弹流闪温小于Blok闪温。

Figure 11. The flash temperature of helical gear

图11. 斜齿轮闪温分布

4. 结论

1) 以端面刚度和啮合面为基础，建立了斜齿轮啮合面单位线载荷分析模型，得到了斜齿轮啮合面上每一点的单位线载荷。啮合面单位线载荷分布与转角误差趋势相同，呈方波状分布，方向与齿宽方向成 ${\beta }_b$ 角度，齿廓中部载荷比较大，齿顶和齿根载荷较小。

2) 建立了便于斜齿轮设计和校核的特征坐标系，既包含了齿廓信息，又包含齿宽信息，能够反映斜齿轮危险点信息，特征坐标系能够代表斜齿轮啮合面情况，可以有效地减小斜齿轮设计和校核工作量。

3) 建立了更符合实际工况的斜齿轮热弹流分析模型，得到了斜齿轮接触点的油膜厚度、压力和温度场分布，以及沿特征坐标的闪温分布，高速重载斜齿压力和颈缩靠近出口，沿特征坐标闪温大致呈V形分布，在主动轮啮入点最大，计算结果与Blok闪温分布规律一致，在主动轮齿根和齿廓中部，热弹流闪温大于Blok闪温，在齿顶处，热弹流闪温小于Blok闪温。