Dynamic analysis of torus involute gear including transient elastohydrodynamic effects

1. Introduction

The cylindrical involute gear has the advantages of line-contact meshing, constant working pressure angle, insensitivity to center distance variations and simplicity in use. However, in most cases, transmission performance of line-contact conjugate surfaces is not satisfactory: high sensitivity to machining or mounting error; forced deformation, deviation load and edge contact caused by variation of tooth contact stiffness or bending stiffness. Topological modification technology [1, 2] is usually utilized to produce mismatched meshing of line-contact engaged gears, but actually it is difficult to design and optimize mismatched meshing regions. Mitome [3] invented spherical gears, which geometrically have two types of tooth profiles, convex teeth and concave teeth. Insensitive to machining or mounting error, the spherical gear set allows large axial misalignments without meshing interference. For a spherical gear set, if the spherical center of tip surface or root surface deviates from the gear axis, the gear set can no longer be called a spherical gear set but a torus involute gear set [4]. This type of gear set, easily machined, can also compensate large axial misalignments. As shown in Fig. 1, the tip surface or root surface is respectively an outer or inner torus generated by a circle whose center is not in the gear axis. Each flank of cross-section along face width is still an involute profile generated from the same base circle. For the torus involute gear with convex teeth (hereinafter referred to as convex gear), tooth thickness on the reference circle gradually decreases from the middle to both ends of face width, whereas tooth thickness on the tip circle increases. For the torus involute gear with concave teeth (hereinafter referred to as concave gear), tooth thickness on the reference circle gradually increases from the middle to both ends, while tooth thickness on the tip circle decreases. Although studies on mathematical models and tooth contact analysis have been implemented [5], more sufficient proofs, especially dynamic performance, need to be provided to verify that the torus involute gear (TIG) can work with good meshing characteristics.

As the main excitations of motions along the off-line-of-action (OLOA), friction forces usually couple with backlash and time-varying meshing stiffness to make gear pairs act as nonlinear and time-varying systems. Additionally, power dissipation due to viscous shear of the lubrication oil film along the tooth contact interfaces forms the main source of mesh viscous damping [6]. On the other hand, dynamic meshing forces and slipping velocity fluctuations affect tribological behavior in terms of oil film height, lubrication viscosity and friction forces. Hence it is necessary to understand them together with acting forces and torques, so as to ensure a sufficient transmission lifetime. Currently gear dynamics is usually investigated based on vibration theory [7, 8]. For simplicity, contact for tooth surfaces is modeled through a torsional spring-damper element which acts tangentially to the base circle of each gear. The dynamic mesh impact cannot be calculated accurately in the torsional vibration model. As a numerical approach, the multi-body dynamics method (MBM) could be used to efficiently model contact of gear pairs with acceptable accuracy and considerably less computational effort compared with the finite element method (FEM) [9]. However, meshing friction is usually simplified without considering effects of lubrication oil film in most of multi-body gear models, and contact stiffness is time-independent. Fietkau [10] proposed an efficient transient elastohydrodynamic gear contact simulation method with commercial code Simpack which took into account oil films and elastic deformations directly in the multi-body simulation. But the variation of contact stiffness over mesh is not counted in, and Reynolds equation needs to be solved with the complicated multi-grid method.

In view of tooth shape's complexity of TIG and point-contact pattern, a new efficient multi-body model is in request, which can account for effects of lubrication oil film, tooth deformation and variation of contact stiffness. Therefore, the main motivation of this paper is to build a new multi-tooth dynamics model including transient elastohydrodynamic effects with low computation efforts to understand dynamic characteristics of TIG comprehensively and verify the feasibility of TIG set, which will lay theoretical foundations for the practical application of it.

Fig. 1Torus involute gears

a) Convex gear

b) Concave gear

2. Multi-tooth dynamics model

An approach of Ebrahimi [9] with circumferentially movable teeth is adopted. Fig. 2 shows a model of meshing multi-tooth. Each rigid tooth has one degree of freedom (DOF): the rotation about the rigid gearwheel center whose angular displacement is ${\theta }_i$. The gearwheel has one rotational DOF with respect to the ground, which is denoted by $\gamma$. All teeth and their respective contacts are treated separately. In order to count in elastic deflection of the tooth, each tooth is connected with the gearwheel by a rotatable spring-damper element. The forces resulted from the flank contact are calculated based on the contact force model including effects of mixed elastohydrodynamic lubrication (EHL).

2.1. Calculation of rotatable spring element stiffness

In order to compute the stiffness of a rotatable spring-damper element, its relation with tooth deformation needs to be figured out. As shown in Fig. 3, for simplicity, misalignment of the LOA induced by flank contact deformation is neglected. i.e. contact points only translate along the LOA over meshing.

Fig. 2Dynamic model of meshing multi-tooth

Fig. 3Meshing point translating along LOA

In terms of geometric characteristics of involutes, the stiffness of a rotatable spring-damper element is calculated as:

1

$k_w=\frac{T}{\mathrm{\Delta }\theta }=\frac{F_n{r}_b}{\mathrm{\Delta }\theta }=\frac{F_n}{\overline{CD}}{r}_b^2=k{r}_b^2.$

Fig. 4Approximation of TIGs

a) Convex gear

b) Concave gear

Fig. 5Model for calculating tooth deformation

As is known, a gear with continuous linear shifting is called the beveloid gear or conical involute gear. In essence, the TIG is a special spur gear with continuous shifting in the second order. To calculate the tooth deformation, each TIG is approximately viewed as a combination of two conical involute gears as depicted in Fig. 4. In the middle cross-section of a TIG, a coordinate system is established as Fig. 5. Then based on elastic mechanics, the relations between Hertzian contact stiffness $k_h$, bending stiffness $k_b$, shear stiffness $k_s$, axial compressive stiffness $k_a$, fillet foundation stiffness $k_f$and Hertzian energy $U_h$, bending energy $U_b$, shear energy $U_s$, axial compressive energy $U_a$, fillet foundation energy $U_f$are respectively gained, which are written as:

2

$U_a=\frac{F_n^2}{2k_a}={\int }_{y_{{}_D}}^{y_{{}_{\beta }}}\frac{F_a^2}{2E{A}_k}d{y}_k\approx {\int }_{r_{{}_f}}^{r_{{}_{\beta }}}\frac{F_a^2}{2E{A}_k}d{r}_k,$

where:

$F_a=F_n\mathrm{s}\mathrm{i}\mathrm{n}\beta ,\mathrm{}{A}_k=\frac{B\left(s_k+s_{kl}\right)}{2\mathrm{c}\mathrm{o}\mathrm{s}\sigma },$

$s_k=s\frac{r_k}{r\mathrm{\text{'}}}-2r_k\left(\mathrm{i}\mathrm{n}\mathrm{v}{\alpha }_k-\mathrm{i}\mathrm{n}\mathrm{v}\alpha \right),$

$s_{kl}=\left(s+2{\epsilon }_l{m}_n\mathrm{t}\mathrm{g}\alpha \right)\frac{r_k+{\epsilon }_l{m}_n}{r\text{'}}-2\left(r_k+{\epsilon }_l{m}_n\right)(\mathrm{i}\mathrm{n}\mathrm{v}{\alpha }_{kl}-\mathrm{i}\mathrm{n}\mathrm{v}\alpha )$

${\alpha }_k=arccos\frac{r_b}{r_k},\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}{\alpha }_{kl}=\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{c}\mathrm{o}\mathrm{s}\frac{r_b}{r_k+{\epsilon }_l{m}_n},\sigma =\mathrm{a}\mathrm{r}\mathrm{c}\mathrm{t}\mathrm{g}\left(\frac{{\epsilon }_l{m}_n}{B}\right),\mathrm{}{\epsilon }_l=\frac{r_{al}-r_a}{m_n},$

$U_b=\frac{F_n^2}{2k_b}={\int }_{y_{{}_D}}^{y_{{}_{\beta }}}\frac{M_k^2}{2E{I}_k}d{y}_k\approx {\int }_{r_{{}_f}}^{r_{{}_{\beta }}}\frac{M_k^2}{2E{I}_k}d{r}_k,$

$M_k=F_n\left(y_{\beta }-r_k\right)\mathrm{c}\mathrm{o}\mathrm{s}\beta -F_n{b}_{\beta }\mathrm{s}\mathrm{i}\mathrm{n}\beta ,\mathrm{}{I}_k=\frac{B}{\left(s_{kl}-s_k\right)\mathrm{c}\mathrm{o}\mathrm{s}\sigma }\left(\frac{s_{kl}^4}{48}-\frac{s_k^4}{48}\right),$

$U_s=\frac{F_n^2}{2k_s}={\int }_{y_{{}_D}}^{y_{{}_{\beta }}}\frac{1.2F_b^2}{2G{A}_k}d{y}_k\approx {\int }_{r_{{}_f}}^{r_{{}_{\beta }}}\frac{1.2F_b^2}{2G{A}_k}d{r}_k,\mathrm{}{F}_b=F_n\mathrm{c}\mathrm{o}\mathrm{s}\beta .$

Stiffness related to the elasticity of root foundation is as follows [11]:

3

$\frac{1}{k_f}=\frac{\mathrm{c}\mathrm{o}\mathrm{s}\sigma \mathrm{c}\mathrm{o}{\mathrm{s}}^2\beta }{EB}\left\{5.306{\left(\frac{L_f}{H_f}\right)}^2+2\left(1-\nu \right)\left(\frac{L_f}{H_f}\right)+1.534\left(1+\frac{0.4167\mathrm{t}\mathrm{a}{\mathrm{n}}^2\beta }{1+\nu }\right)\right\},$

where $L_f=y_{\beta }-y_D-b_{\beta }\mathrm{t}\mathrm{a}\mathrm{n}\beta$; $H_f=2b_D$.

The total tooth stiffness is as follows:

4

$k=\frac{1}{\frac{1}{k_b}+\frac{1}{k_s}+\frac{1}{k_a}+\frac{1}{k_f}}.$

Over a mesh cycle the stiffness of a single tooth can be calculated from the effective root to addendum by treating the pressure angle of each contact point as a variable. Then Eq. (1) is employed to give the time-varying stiffness of a rotatable spring-damper element.

3. Model verification

As the measurement of contact forces is very difficult or even impossible for many technical problems, so far there is no direct way of measuring dynamic contact forces of teeth. To validate the simulation model, the basic experimental impact investigation is performed using one gear and one simple impact body. For the simplicity of the experimental setup, the investigated gear is not rotating. A cuboid is used as the impact body.

The recoil velocities of the cuboid are indirectly measured using the oscilloscope, and these velocities are compared with simulations. According to the impulse theorem, when there is no friction force, only the gravity force and contact force contribute to the velocity change of the impact cuboid. An assumption can therefore be made that if recoil velocities of the cuboid agree, the contact forces must agree too. The air supplied slide rail is used to guarantee zero friction. As illustrated in Fig. 7, the cuboid translates on an air supplied slide rail in direction of the LOA and collides with a tooth. The slide rail is mounted on a precision rotary stage that allows to adjust the alignment angle $\alpha$. The rotary stage itself is mounted on a frame that allows to adjust the alignment angle $\varphi$. Electrical signals can be acquired to determine the time interval $t_n$ between two successive collisions. So the initial velocity for the subsequent collision is $V_{in}=0.5g{t}_n\mathrm{s}\mathrm{i}\mathrm{n}\varphi$, and its recoil velocity is $V_{out}=0.5g{t}_{n+1}\mathrm{s}\mathrm{i}\mathrm{n}\varphi$. The investigated gears have the parameters shown in Table 1. Table 2 demonstrates the measured and simulated results. From this table we note that the simulated results agree well with the measured ones.

Fig. 6Flow chat of friction coefficient calculation

Fig. 7Impact experiment set up (1 – air supplied rail; 2 – inclined bracket; 3 – inclination angle adjustment; 4 – vertical bracket; 5 – rotary stage; 6 – fixture; 7 – TIG; 8 – impact cuboid)

4. Results and discussion

A numerical example of a TIG pair is provided to validate the above approach. From parameters in Table. 1, the tooth surfaces for engaged gears are constructed from the mathematical models, and then 3-D models can be built. Fig. 8 shows the tooth surfaces of engaged TIGs. Then TCA is conducted for each rotating angle of the convex gear (i.e. gear 1) based on theory of differential geometry. As displayed in Table 3, principle curvatures, comprehensive curvature radius, semi-major axis and semi-minor axis are obtained for each contact ellipse. Using Fourier fitting, they can be represented as functions of angular displacement of the gear wheel.

Through polynomial fitting, spring-damper stiffness and contact stiffness can be represented as functions of pressure angles. Fig. 9 displays spring-damper stiffness varying with pressure angle at each contact point, and contact stiffness of the TIG pair is shown in Fig. 10. From above stiffness, the mesh stiffness can be calculated.

To verify the analytical result of mesh stiffness, an 8-9-tooth finite element model is built as follows (See Fig. 11). The preprocess of the finite element model is completed in Hypermesh; Either of the engaged gears can only rotate about its own axis; Teeth are connected with a central node by massless rigid elements; The driving gear has a certain angular velocity, and a torque 20 N·m is exerted on the other gear. Contact algorithm is explicit surface to surface contact. The dynamic finite element analysis (FEA) is performed in ABAQUS to obtain the comprehensive deformation of TIGs. The mesh stiffness is then computed, and the comparison of two methods is made in Fig. 12.

Fig. 8Tooth surfaces of TIGs

a) Convex gear

b) Concave gear

On the basis of above work, a 3-D multi-rigid-body dynamic model is built in ADAMS as follows. Either of the gearwheels has one rotational DOF with respect to the ground; Each tooth can rotate about its gearwheel center, which is connected with the gearwheel by a rotatable spring-damper element; Contact pairs are defined according to the mating relationship of teeth; The pinion (gear 1) is driven by a constant angular velocity, and a constant torque $T_2=\mathrm{}$20 N·m is simultaneously exerted on the gear (gear 2).

Fig. 9Stiffness of rotatable spring-damper

Fig. 10Contact stiffness of TIG pair

Fig. 11Finite element model of TIG pair

Fig. 12Mesh stiffness of TIG pair

4.2. Contact force

Fig. 15 shows the normal forces for 3 successive meshing tooth pairs. From this figure the transmission continuity of the gear set can be checked. The duration for double pair meshing over one mesh cycle is $\mathrm{\Delta }t\text{,}$ and the mesh cycle period is $\mathrm{\Delta }T\text{.}$ So the contact ratio $\mathrm{\Omega }=1+\mathrm{\Delta }t/\mathrm{\Delta }T=\text{1+}\text{0.0188}/\text{0.044}=\text{1.4273,}$ which shows that transmission continuity is achieved. Moreover, it can also be seen from Fig. 15 that the normal contact force is nearly independent of lubrication viscosity when considering effects of EHL.

Fig. 16 is the result obtained by superposing the normal contact forces of all tooth pairs. The magnitude of normal contact forces can be verified preliminarily. Utilizing the conventional gear design approach, the normal contact force is $F_n=T_2/r_{b2}=\text{20000}/\text{57.56}=$ 347.49 N, and the average value of our model is 353.35 N. The increase in average normal contact force is induced by the effect of impacts during meshing in and out as well as friction forces. Besides, Fig. 16 indicates that when the preceding tooth pair meshes out, at the pitch or the subsequent one meshes in, there exists a major fluctuation of normal contact force.

Fig. 13DTE for different original viscosity (n= 60 r/min)

Fig. 14DTE for different angular velocity

Fig. 15Normal contact force of Tooth pair 3-5 (n= 60 r/min)

Fig. 16Total normal contact force (η0= 0.012 Pa·s, n= 60 r/min)

Fig. 17Friction force of 4th tooth pair (η0= 0.012 Pa·s, n= 60 r/min)

Fig. 18Friction force of 4th tooth pair (η0= 0.012 Pa·s)

Fig. 17 displays the friction force of the 4th tooth pair for different original viscosity (The friction force of the 3rd pair and the 5th pair are also plotted in black thin line). As depicted in Fig. 17, at the pitch point the relative slip velocity between tooth flanks changes its sign and therefore the friction force does as well. As a result, the normal contact force tends to decline (See Fig. 16). Moreover, original viscosity slightly affects the friction force. The friction force rises with the falling original viscosity because the asperity contact force rises.

Fig. 18 displays the friction force of 4th tooth pair along LOA for different angular velocity. It can be seen from Fig. 18 that the angular velocity has great effect on friction forces. As angular velocity rises, the relative slipping velocity and entraining velocity of engaged flanks increase, and consequently the fluctuation and peak value of friction force tend to increase. However, the general tendency of all plots are similar.

Frequency-domain analysis for normal contact force and angular acceleration of the driven gear are shown in Fig. 19. From the figures the primary resonances are both found at 40.05 Hz, nearly equal to gear mesh frequency ($z_1{n}_1/\text{60}$ where $z_1$ is the tooth number for gear 1, and $n_1$ is the rotational speed for gear 1), which is associated with a mode having the LOA transverse motions and the torsional motions of the wheels.

Fig. 19Results in frequency domain (η0= 0.012 Pa·s, n= 60 r/min)

a) Normal contact force

b) Angular acceleration of the driven gear

4.3. Lubrication oil film height

Fig. 20 displays the average lubrication oil film height of the 4th tooth pair for different original viscosity (Those of the 3rd pair and the 5th pair are also plotted in black thin line). As shown in Fig. 20, original viscosity affects the lubrication oil film height markedly. With the increasing of original viscosity, the average lubrication oil film height increases. For a generic tooth pair, as the normal contact force rises, it begins to mesh in, and the average lubrication oil film height declines gradually due to the increasing normal contact force. When the subsequent tooth mesh in, the contact normal force decreases, and therefore the average lubrication oil film height rises.

Fig. 20Average lubrication oil film height of 4th tooth pair (n= 60 r/min)

Fig. 21Average lubrication oil film height of 4th tooth pair (η0= 0.012 Pa·s)

Fig. 21 illustrates the average lubrication oil film height for different angular velocity. As depicted in Fig. 21, angular velocity influences the lubrication oil film height evidently. Obviously, the increasing of angular velocity causes the entraining velocity of engaged flanks to increase. As a result, the average lubrication oil film height tends to rise and fluctuate clearly. All in all, either angular velocity or original viscosity affects lubrication oil film height more strongly than it does to contact forces.

5. Conclusions

A new efficient model for gear contact simulation in mixed EHL is presented and integrated in multi-body simulation environment. Lubrication oil film height and tooth deformations are directly counted in this model. Time-varying contact stiffness and rotatable spring stiffness are both taken into account. Combined with TCA, the computation of friction coefficient is implemented with high efficiency. Subroutines for the computation of normal tooth contact forces and friction coefficient are developed and called by the numerical integrator in commercial code ADAMS. A 3-D multi-body model of a TIG set is employed and verified by an impact experiment. The results indicate that the continuity and precision of transmission can be guaranteed using this type of gear set. Important values like dynamic transmission error, normal and tangential contact forces and average film heights are calculated with high efficiency for different conditions, which lay foundations for application of TIG. Besides, it is feasible to analyze the tribo-dynamic behavior of other point-contact gears using the approach in this paper.