Comparison of meshing characteristics and dynamic behavior of hypoid gear drive under different load conditions

Abstract

The study of loaded meshing characteristics and dynamic behavior is significant for evaluating the stability of hypoid gear under actual working conditions. Loaded tooth contact analysis (LTCA) of a hypoid gear drive under different loads is carried out. The equivalent meshing force, loaded transmission error (LTE), comprehensive elastic deformation, and time-varying meshing stiffness (TVMS) are investigated. Results show that the LTE and TVMS are sensitive to load, and the LTE and TVMS curves show obvious asymmetric as the load increases. Under low load conditions, the actual contact ratio increases significantly with the increase of load, but tends to be flat when the load reaches a certain range. The obtained TVMS and LTE are fitted by higher Fourier series and introduced into an 8-DOF dynamic model. By comparing the time history diagram, phase diagram, Poincaré section, and spectrum, the dynamic characteristics of the system under different loads are contrasted. On low load condition, the system may exhibit chaotic motion in certain meshing frequency bands and damping ratio intervals. While the tooth backlash increases to a critical value, the system may change to chaotic motion from periodic motion. The results provided a basis for improving the meshing performance of hypoid gear.

Keywords

Hypoid gear drive LTCA meshing characteristics dynamics behavior load condition

Introduction

At present, the study of hypoid gears can be roughly divided into two directions. One focuses on optimizing the machining process to obtain high-precision tooth profiles.^1–8 Through the static geometric contact analysis, the gear blank parameters and machine-tool settings can be optimized. Tooth surface stress, load distribution, transmission error can also be obtained by the quasi-static finite element analysis.^9–15 However, the results achieved by this “static analysis” cannot fully reflect the vibration and noise levels of the gear drive during entire meshing process. Another direction focuses on establishing a multi-freedom dynamics model of the gear system. With the consideration of the nonlinear factors, such as the excitation of the time-varying stiffness, dynamic transmission error, and backlash nonlinear, the dynamics response and vibration features of the gear pair can be researched.^16–25 Results obtained by this “dynamics analysis” can reflect the dynamic behavior of the gear system and analyze the causes of vibration and noise. Nevertheless, some of the nonlinear factors are generally expressed using approximate formulas, like the meshing stiffness, which may lead to inaccurate results.^26,27 So far, the above two directions have usually been studied independently, so the static analysis and dynamics modeling of hypoid gears have not been well connected.

Stiffness fluctuation is one of the main excitations of the noise and vibration of gear system.^28–30 Accurate determination of meshing stiffness has always been a major problem in the dynamics research of gears. For hypoid gears, dynamics studies generally established a mathematical model of the gear drive to investigate the effects of the stiffness, damping, and inertial mass on the dynamic characteristics. Among the three types of parameters, the most significant difference between the hypoid gears and the spur gears is the meshing stiffness. Main reason is that the magnitude and direction of the meshing force of spur gears can be determined easily, and there are reliable methods and empirical formulas, while the tooth profiles and meshing process of hypoid gears are both complicated.¹ During the meshing process of a hypoid gear pair, the direction of the meshing force and the position of the meshing point vary with the gear’s rotation angle and the magnitude of the load. Compared with the spur gears, it is relatively rare to study how to precisely express the time-varying meshing stiffness of hypoid gears. Up to now, the calculation of the cylindrical gears’ meshing stiffness has been generally solved. Kiyohiko derived the formulas for calculating the bending stiffness and contact stiffness of involute cylindrical gears.³¹ An empirical formula for calculating the single-tooth stiffness for the stand full-depth or addendum modified involute gears was proposed by Kuang and Yang,³² which has been widely used to calculate the meshing stiffness of spur gears.

Today, it is still difficult to measure the meshing stiffness of hypoid gears by the experimental methods, while the finite element analysis method can precisely simulate the instantaneous contact state of the gear drive. For computing the meshing stiffness of spiral bevel gears, Gosselin et al. presented a calculation method based on the finite strip theory, which obtained the displacement curves along the tooth height and tooth length direction, while no stiffness curve was achieved.³³ Mennem used the finite element method to compute the tooth contact softness under different loads and finally worked out the time-varying stiffness.³⁴ Wang et al. established a generalized nonlinear time-varying (NLTV) dynamic model of a hypoid gear pair and approximated the gear meshing stiffness using a constant and a triangular series.²⁶ Mohammadpour et al. and Karagiannis et al. calculated the meshing stiffness of hypoid gears by the finite element software CALYX. Then they simplified the meshing stiffness to a triangular series form and finally obtained the empirical formula of the meshing stiffness.^35,36 Fang established a meshing model of a hypoid gear pair based on the static force state of the gears and calculated the dynamic load, while this model did not consider the time-varying characteristics of the meshing stiffness.³⁷ According to the basic theory of meshing stiffness, Yang et al. calculated the stiffness of a single pair of teeth, then worked out the integrated meshing stiffness of the gear pair by overlaying the stiffness of the simultaneously meshing teeth.³⁸

Existing finite element methods for calculating meshing stiffness of hypoid gears are generally based on measuring the elastic deformation at the calculating nodes on tooth surfaces. The bending deformation of the tooth profile and the shear deformation of the tooth root are not considered, which may result in less precise stiffness data. Obviously, the deformation of the tooth surface is called the contact deformation more suitable than the meshing deformation. Wang presented a more accurate method to calculate the meshing stiffness of spiral bevel and hypoid gears by computing the dynamic loaded transmission error and the quasi-static unloaded transmission error.¹⁸ Since the rotation angle can be measured more precisely than the elastic deformation by the finite element method, this calculation principle is more accurate than the traditional method. In his PhD dissertation, Wang stated that this method was implemented by the Contact Analysis Program Package (CAPP), so there is no detailed description of the actual measuring process in the dissertation. This calculation method and the CAPP were also used in other scholars’ research.^19,38 Among the existing literatures, the comparison of the meshing stiffness measured by the finite element method under different load conditions was rare. Besides, the reports on the accurate fitting of hypoid gears’ meshing stiffness data are also scarce.

Dynamics study can not only analyze and predict the vibration and noise levels of the gear system, but also put forward new requirements for the design and processing of the gears. The vibration performance of a hypoid gear pair will directly affect the stability and reliability of the transmission system, so it is very valuable to analyze its dynamics. Most dynamics researches usually study the effects of different nonlinear factors on the dynamic responses by establishing a multi-freedom dynamics model. Since hypoid gears’ meshing stiffness is difficult to measure and calculate, the first or second order form of the Fourier series is usually used to represent it, which may result in inaccurate analysis results. Due to the lack of research on precisely measuring the meshing stiffness of hypoid gears, it is rare to import the accurate stiffness data into the dynamics model. In addition, the analysis of comparing the influence of the meshing stiffness under different loads on the dynamics characteristics of hypoid gears is also rare now.

In this research, we give a detailed mathematical model for computing the meshing stiffness by measuring the contact forces and the dynamic and static transmission error of a hypoid gear pair. The static transmission error is calculated by the TCA method. The actual rotational displacement of the gear and the pinion under different loads can be obtained based on the post-processing results of ABAQUS software, then the dynamic transmission error under different can be worked out based on the difference value of the actual rotation angle and the theorical angle of the pinion. By using the same method to calculate the spur gears’ meshing stiffness, the reliability and accuracy of the used method are verified. In order to compare the influence of different meshing stiffness on the dynamics of the gear system, the stiffness data is fitted by an eighth-order Fourier series and imported into an 8-dof dynamics model of the gear pair. Simultaneously, the accuracy and goodness of the results fitted by the second-order and eighth-order Fourier series are also compared. By contrasting the dynamics characteristics under different loads, the effects of different meshing stiffness on the motion features of the system are compared. Through the meshing stiffness, the connection between the static analysis of the meshing teeth and the dynamics analysis of the gear system has been established in this research.

Mathematical description of the meshing parameters

Equivalent meshing force (EMF)

Due to the high transmission ratio of hypoid gear drive, there are multiple pairs of teeth generally in contact at the same time during the meshing process. As shown in Figure 1, ${\vec{r}}_{i} (r_{x}^{i}, r_{y}^{i}, r_{z}^{i})$ is the radial vector of the application point of the meshing force on the i-th tooth pair, which direction is from the origin of the meshing coordinate system to the midpoint of the instantaneous contacting ellipse. ${\vec{f}}_{i} (f_{x}^{i}, f_{y}^{i}, f_{z}^{i})$ is the vector of the equivalent meshing force, which application point is the midpoint of the contacting ellipse. The direction of ${\vec{f}}_{i}$ points to the outside of the tooth surface along the normal direction at the application point. The magnitude $| {\vec{f}}_{i} | (| f_{x}^{i} |, | f_{y}^{i} |, | f_{z}^{i} |)$ of ${\vec{f}}_{i}$ can be computed based on the results of the meshing force in LTCA. The magnitude of the resultant force ${\vec{F}}_{Q} (F_{x}, F_{y}, F_{z})$ of all equivalent meshing force on the contacting surfaces can be computed based on

| {\vec{F}}_{Q} | = \sqrt{{F_{x}}^{2} + {F_{y}}^{2} + {F_{z}}^{2}}

(1)

where F_x, F_y, F_z are the component of ${\vec{F}}_{Q}$ on the x-, y-, and z-axis, respectively, which can be calculated based on

F_{k} = \sum_{i = 1}^{q} | f_{k}^{i} |, k = x, y, z

(2)

where q is the number of the contacting tooth surfaces. The unit vector ${\vec{L}}_{Q} (n_{x}, n_{y}, n_{z})$ of the direction of ${\vec{F}}_{Q}$ can be worked out according to

n_{k} = \frac{F_{k}}{| {\vec{F}}_{Q} |}, k = x, y, z

(3)

Figure 1.

Equivalent meshing force and application points.

Equivalent meshing point (EMP)

Based on the radial vector ${\vec{r}}_{i} (r_{x}^{i}, r_{y}^{i}, r_{z}^{i})$ of the meshing force and the resultant force ${\vec{F}}_{Q} (F_{x}, F_{y}, F_{z})$ of the gear pair, according to the principles of mechanical balance, the radial vector of the application point P_Q[x_Q, y_Q, z_Q] of ${\vec{F}}_{Q}$ in the meshing coordinate system can be worked out based on

x_{Q} = \frac{\sum_{i = 1}^{q} r_{x}^{i} | {\vec{f}}_{i} |}{\sum_{i = 1}^{q} | {\vec{f}}_{i} |}, y_{Q} = \frac{M_{x} + F_{z} x_{Q}}{F_{z}} z_{Q} = \frac{M_{y} + F_{z} x_{Q}}{F_{x}}

(4)

herein, M_x and M_y are the torque of ${\vec{F}}_{Q}$ around the x- and y-axis, which can be obtained by LTCA. P_Q can be regarded as the instantaneous equivalent meshing point of the gear pair, which does not exist actually, but is calculated by the equivalent meshing force and the mechanical principles, and can reflect the meshing state of the gear pair. Considering the effects of the errors from the meshing state, as the gear pair rotates, P_Q moves in a small space near the contacting area theoretically. Due to the periodic characteristics of the meshing process, the actual trajectory of P_Q is a closed curve.¹⁸

Transmission error (TE)

Since the mating surfaces of hypoid gear pair are in local conjugate point contact state, the TE of the gear drive occurs inevitably during the meshing process, which is generally defined as the difference between the actual and the theoretical rotation angle of the driven gear, as follows

δ (φ_{P}) = (φ_{G} - φ_{G}^{0}) - \frac{z_{P}}{z_{G}} (φ_{P} - φ_{P}^{0})

(5)

where z_P and z_G are the teeth number of the pinion and gear, respectively. $φ_{P}^{0}$ and $φ_{G}^{0}$ are the initial rotation angle of the pinion and gear, respectively. $φ_{P}$ and $φ_{G}$ are the actual rotation angle of the pinion and gear, respectively. In practice, the magnitude of the load affects the position of the contact area and the instantaneous meshing state. Therefore, the loaded TE under different loads is different from the theoretical TE, which should be calculated by LTCA.

Comprehensive elastic deformation (CED)

The CED u_n of the contact area of a hypoid gear pair consists of the contact deformation u_c of the tooth surface, the bending deformation u_b of the tooth profile, the shear deformation u_s of the tooth, and the torsional deformation u_t of the basement of the gear, that is

u_{n} = \sum_{i = p, g} (u_{c}^{i} + u_{b}^{i} + u_{s}^{i} + u_{t}^{i})

(6)

where the subscript p and g represent the pinion and the gear respectively. It is difficult to measure the above four kinds of deformations actually, and usually only the tooth surface contact deformation u_c can be directly extracted from the LTCA results, which is generally used as the elastic deformation of the gear pair in some previous studies.

In order to obtain u_n more precisely, a novel calculation principle of CED by extracting the coordinates of the instantaneous contacting area and the rotation angle of the gear pair is proposed in this paper. As shown in Figure 2, during the rotation process of the gear pair, the actual displacements of the meshing points at one moment are composed of the CED and the rotational displacement of the tooth. In Figure 2, the tooth profiles of the pinion and the gear are the projections of the gears in the cross-section of the pinion and the gear, which are perpendicular to the rotation shafts. In fact, the projections are perpendicular to each other, so the points in Figure 2 should be described by the three-dimensional coordinates [x_i, y_i, z_i] (i = P, G). The points $P_{0} [x_{0}^{P}, y_{0}^{P}, z_{0}^{P}]$ and $G_{0} [x_{0}^{G}, y_{0}^{G}, z_{0}^{G}]$ correspond to the initial position of the meshing points at the initial moment t₀. $P_{t} [x_{t}^{P}, y_{t}^{P}, z_{t}^{P}]$ and $G_{t} [x_{t}^{G}, y_{t}^{G}, z_{t}^{G}]$ correspond to the actual position of the meshing points at the contacting moment t. ${P'}_{t} [x'_{t}^{P}, y'_{t}^{P}, z'_{t}^{P}]$ and ${G'}_{t} [x'_{t}^{G}, y'_{t}^{G}, z'_{t}^{G}]$ correspond to the theoretical position of the meshing points at moment t. The pinion and the gear respectively rotate around the X_H and Y_H axes of the meshing coordinate S_H[X_H, Y_H, Z_H]. The free rotation radius of P_t and G_t can be calculated based on

r_{0}^{P} = \sqrt{{(y_{0}^{P})}^{2} + {(z_{0}^{P})}^{2}}, r_{0}^{G} = \sqrt{{(x_{0}^{G})}^{2} + {(z_{0}^{G})}^{2}}

(7)

Figure 2.

Actual and theoretical displacements of the meshing points on tooth surfaces.

The theoretical displacement s_P(t) and s_G(t) of ${P'}_{t}$ and ${G'}_{t}$ at moment t can be calculated by

s_{i} (t) = 2 r_{0}^{i} (t) \cdot \sin (\frac{φ_{i} (t)}{2}), (i = P, G)

(8)

where the φ_P(t) and φ_G(t) are the actual rotation angle of the pinion and the gear respectively. The actual displacement S_P(t) and S_G(t) of $P_{t}$ and $G_{t}$ at the time t can be calculated by

S_{i} (t) = \sqrt{{(x'_{t}^{i} - x_{0}^{i})}^{2} + {(y'_{t}^{i} - y_{0}^{i})}^{2} + {(z'_{t}^{i} - z_{0}^{i})}^{2}}, (i = P, G)

(9)

At the time t, it is supposed that the numbers of the contacting points on the tooth surfaces of the pinion and the gear are m and l respectively, thus the comprehensive elastic deformation u_n of the gear pair can be computed based on

\begin{matrix} u_{n} (t) = u_{n}^{P} (t) + u_{n}^{G} (t) = \frac{1}{m} \cdot \sum_{i = 1}^{m} | (S_{i}^{P} (t) - s_{i}^{P} (t)) | \\ + \frac{1}{l} \cdot \sum_{j = 1}^{l} | (S_{j}^{G} (t) - s_{j}^{G} (t)) | \end{matrix}

(10)

The calculation process of CED can be divided into the following steps,

Step 1. Select the j-th (j = 1, 2, …, n, n is the maximum number of teeth engaged simultaneously) teeth pair after the gear drive engages stably as the studied teeth. Extract the rotation angle of the gear and the pinion at each analysis step from the LTCA outputs.

Step 2. Select the time t₀ after the gear drive meshes stably as the initial time, and the j-th teeth pair must not be engaged at t₀. Extract the coordinates $[x_{0}^{P}, y_{0}^{P}, z_{0}^{P}]$ and $[x_{0}^{G}, y_{0}^{G}, z_{0}^{G}]$ of the nodes on the surfaces of the j-th teeth pair. The rotation angles of the gear and the pinion at t₀ are recorded as φ_G(t₀) and φ_P(t₀).

Step 3. During the meshing process of the j-th teeth pair, at each analysis step, take all contacting nodes as a set and extract their coordinates $[x_{t}^{P}, y_{t}^{P}, z_{t}^{P}]$ and $[x_{t}^{G}, y_{t}^{G}, z_{t}^{G}]$ , and record the rotation angle φ_G(t) and φ_P(t), respectively.

Step 4. During the meshing process of the j-th teeth pair, at each analysis step, based on the above parameters and equations (7)–(9), the actual and theoretical displacements Sⁱ(t) and sⁱ(t) (i = P, G) of the contacting nodes can be worked out. The CED u_n of the gear pair at can be calculated based on equation (10).

Time-varying meshing stiffness (TVMF)

According to the definition of the meshing stiffness of the gear pair,³² the time-varying meshing stiffness k_n(t) of the hypoid gear drive can be calculated based on

k_{n} (t) = \frac{| {\vec{F}}_{Q} (t) |}{u_{n} (t)}

(11)

Verification of the LTCA method

The exact formulas for computing the meshing stiffness of the hypoid gear drive are rarely reported in current literatures, while those of the spur gear pair have been widely adopted. To verify the proposed principle, the same method is used to work out the meshing stiffness of a spur gear drive, and the results are compared with those computed by Kuang and Yang.³² The parameters of the studied spur gear pair are given in Table 1. The instantaneous meshing force on the contact area obtained by LTCA are shown in Figure 3.

Table 1.

Parameters of the studied spur gear pair and the settings of LTCA.

Parameter	Pinion	Gear	Parameter	Pinion	Gear
Number of teeth	18	31	Elastic modulus (N mm⁻²)	2.06e5
Module (mm)	6		Poisson’s ratio	0.3
Pressure angle (deg)	20		Density (kg m⁻³)	7.85e3
Actual center distance (mm)	152.345		Friction coefficient	0.05
Addendum modification coefficient (mm)	0.5	0.5	Contact type	Hard contact
Tooth width (mm)	50	50	Element type	C3D8R
Addendum coefficient	1	1	Analysis type	Implicit, statics analysis
Tip clearance coefficient	0.25	0.25	Slippage	Finite sliding
Load torque (Nm)	0	30	Analysis step size (s)	0.005
Rotation speed (rad/s)	0.5	0	Number of elements	37,980	106,950

Figure 3.

Meshing force on the tooth surfaces of a pair of spur gears.

The meshing stiffness of a single tooth pair and the comprehensive meshing stiffness of the gear drive calculated according to the principle of this paper and the formulas proposed by Kuang and Yang are shown in Figure 4(a) and (b) respectively, where “Km1” and “Km2” represent the meshing stiffness of the first and second tooth pair computed by the classic formulas and “FEM” represents that worked out by the proposed principle in this paper. It is observed that the results obtained by the two methods are essentially same, thus the accuracy of the proposed principle to compute the meshing stiffness can be verified.

Figure 4.

Meshing stiffness of the studied spur gear pair: (a) meshing stiffness of single tooth pair and (b) comprehensive meshing stiffness of spur gear drive.

In order to verify the accuracy of the modeling method and the LTCA process, the LTCA simulation under 30 Nm and the rolling test are respectively implemented. Due to the small load required for the rolling test, the LTCA simulation load is set to 30 Nm. As shown in Figure 5, it can be seen that the contact patterns obtained by LTCA are consistent with the results of the rolling test, and due to the small load condition, the contact pattern is close to the toe of the tooth, which will move toward the heel under general load. Therefore, the finite element model of the gear pair and the LTCA principle are proven to be effective, and the meshing parameters worked out by LTCA can be used in the proposed calculation principle.

Figure 5.

Contact patterns on gear tooth surface obtained by LTCA and the rolling test: (a) contact pattern obtained by LTCA and (b) contact pattern obtained by the rolling test.

Time-varying meshing characteristics

A hypoid gear pair processed by the face-milled method is taken as the case study in this paper, which geometric parameters are given in Table 2. A series of LTCA cases in which the gear pair is driven at 100–6000 Nm load torque and 60 r/min rotate speed are carried out.

Table 2.

Parameters of the studied hypoid gear pair.

Parameters	Pinion	Gear	Parameters	Pinion	Gear
Number of teeth	10	41	Face angle (deg)	20.6584	74.8365
Module (mm)	4.741		Root angle (deg)	14.1923	68.0107
Shaft angle (deg)	90		Pitch apex beyond crossing point (mm)	23.1139	−3.1263
Offset distance (mm)	−31.8		Face cone apex beyond crossing point (mm)	11.9258	−3.6982
Tooth width (mm)	33.7277	28	Root cone apex beyond crossing point (mm)	21.8026	−1.4089
Mean spiral angle (deg)	49.9833	29.3980	Mean addendum (mm)	5.9348	1.2156
Pressure angle (deg)	21		Mean dedendum (mm)	2.3381	7.0574
Pitch angle (deg)	15.2889	73.6761	Working depth (mm)	14.84
Mass (kg)	0.3303	1.7743	Mean backlash between contacting teeth	0.15
Radius of gyration (mm)	27.6	80.7	Moment of inertia (kg mm²)	251.8	11,570

Taking the results under 4000 Nm load torque as an example, the position of the equivalent meshing points is shown in Figure 6, and the EMF, loaded TE, CED, and TVMF are shown in Figure 7.

Figure 6.

Equivalent meshing points in the meshing coordinate system under 4000 Nm condition: (a) sight 1 of equivalent meshing points and (b) sight 2 of equivalent meshing points.

Figure 7.

Time-varying meshing parameters of the hypoid gear drive under 4000 Nm: (a) equivalent meshing force, (b) loaded transmission error, (c) comprehensive elastic deformation, and (d) time-varying meshing stiffness.

Because of the periodicity of the meshing state, the ideal trajectory of the equivalent meshing points is a closed curve. Due to the errors in the modeling and analyzing process in LTCA, the equivalent meshing points obtained in this research are located in a small plane in the meshing coordinate system, as shown in Figure 6. It can be seen from Figure 7 that the EMF, loaded TE, CED, and TVMF all change periodically, and there is no obvious sudden change during the meshing process. As shown in Figure 4, there is a significant abrupt change between 620 and 1160 MN/m of the meshing stiffness of the gear drive. Therefore, the meshing process of the hypoid gear pair is relatively stable. The fluctuation frequency of the meshing parameters is consistent with the meshing frequency of the gear drive, which is 10 Hz in this case.

Loaded TE and actual contact ratio

The obtained loaded TE of the gear pair under different loads is shown in Figure 8(a). It can be seen that the absolute value of the loaded TE increases as the load rises, while the increment decreases gradually. The reason is that the increase of the load makes the actual contact pattern on tooth surface deviate from the ideal contact area, which moves toward the toe and the heel of the tooth, thus the loaded TE increases. However, as the load rises, the increment of the deviation of the contact pattern decreases gradually, while the actual contact ratio increases (as shown in Figure 8(b)), which enhances the stability of the transmission state, so the increment of the loaded TE decreases as the load rises. The influence of the variation of the load on the transmission error of the hypoid gear pair is significant, and it is necessary to consider this effect on the dynamic behavior of the gear drive.

Figure 8.

Loaded TE and actual contact ratio of the studied gear pair: (a) loaded transmission error, (b) actual contact ratio, and (c) contact patterns on gear tooth surface under different loads.

The actual contact ratio in Figure 8(b) varies from 1.38 to 3.01, and the relatively large range of the variation is due to the load torque of the finite element analysis of the gear drive set from 100 to 6000 Nm. From Figure 8(c) it can be seen that as the load increases, the contact pattern extends from the middle of the tooth surface toward the toe and the heel of the tooth. So that the axial contact ratio increases significantly, which leads to the actual contact ratio enhanced obviously. In practice, the load torque is usually in the range of tens to hundreds of Nm, so the variation range of the actual contact ratio will be smaller.

Considering the periodic feature of the meshing state, the Fourier series expansion is used to fit the loaded TE under different loads in this study, which can be described as

{\bar{e}}_{n} (t) = e_{m} + \sum_{i = 1}^{N} (e_{a i} \cos (i Ω_{h} t + ϕ_{e}) + e_{b i} \sin (i Ω_{h} t + ϕ_{e}))

(12)

where e_m is the average of the loaded TE, e_ai and e_bi are the magnitude of the i-th cosine and sine harmonic component of the loaded TE, ϕ_e is the initial phase of the harmonic component, Ω_h is the meshing frequency. The “R-square” is used to evaluate the goodness of the fitting, which ranges from 0 to 1, and the closer to 1, the better the fitting is. It is found by comparison that the six-order Fourier series expansion is best for the fitting of the loaded TE in this study, and the coefficients are given in Table A1.

EMF and TVMS

The obtained EMF and TVMS of the studied gear pair under different loads are shown in Figure 9(a) and (b), respectively. The increase of the EMF is consistent with that of the load torque, and the fluctuating frequency of the EMF equals to the meshing frequency. When the load torque is lower than 2000 Nm, the TVMS rises significantly as the load increases. The reason is that the actual contact ratio of the gear pair increases rapidly as the load rises under lower load torque (as shown in Figure 8(b)), that is, the number of the teeth contacting at the same time increases rapidly, which leads to the TVMS rise obviously. After the load reaches a certain amplitude, as the load increases, the TVMS rises relatively slower.

Figure 9.

(a) Equivalent meshing force and (b) time-varying meshing stiffness of the studied gear pair.

Furthermore, it is worth noting that as the TVMS increases, the maximum of TVMS appears gradually earlier during one meshing cycle. The reason is that as the load increases, the actual contact area expands to the toe and the heel of the tooth. The trajectory of the contact pattern on hypoid gear tooth surface is from the heel to the toe. Therefore, as the load rises, the contact at the heel appears earlier. Because the tooth profile at the heel is thicker, the bending and shear deformation of the tooth is not easy to be generated, so the TVMS at the heel is greater than that at other position of the tooth. Therefore, as the load rises, the first contact pattern moves toward the heel of the tooth, and the maximum of TVMS appears gradually earlier.

As the load changes, the TVMS exhibits strong nonlinear characteristics, which should be concerned in the study of the dynamic behavior of the gear drive. In this paper, the Fourier series expansion is also used to fit the TVMS under different loads, which can be described as

\bar{k} (t) = k_{m} + \sum_{i = 1}^{N} (k_{a i} \cos (i Ω_{h} t + ϕ_{k}) + k_{b i} \sin (i Ω_{h} t + ϕ_{k}))

(13)

where k_m is the average of the TVMS, k_ai and k_bi are the magnitude of the i-th cosine and sine harmonic component of the TVMS, ϕ_k is the initial phase of the harmonic component. It is found by comparison that the three-order Fourier series expansion is best for the fitting of the TVMS in this study, and the coefficients are given in Table A2.

Dynamic model of hypoid gear system

Dynamic model of the gear pair

The simplified dynamics model and contact force model of a hypoid gear pair with the elastic supports is illustrated in Figure 10. p and g indicate the pinion and the gear, respectively, which are simulated as lumped mass and lumped rotational inertia. The shafts are simulated as rigid bodies without mass and the bearings are simulated as springs and damper. T_p and T_g are torques of the pinion and gear, respectively. The X-axis is set as the pinion shaft and the Y-axis is set as the gear shaft. k(τ) is the time-varying meshing stiffness where τ is time. e(τ) is the comprehensive TE along the line of action. c_m is the meshing damping coefficient.

Figure 10.

(a) Dynamics model and (b) contact force model of a hypoid gear pair.

Thus, the entire dynamic system has 12 degrees of freedom, including the bending vibration, axial vibration, torsional vibration, and shimmy of the pinion and the gear along the X-axis, Y-axis, and Z-axis. In Figure 10, “k” represents the stiffness coefficient, “c” represents the damping coefficient, and the subscripts “1” and “2” represent the pinion and the gear respectively. Since the influence of the shimmy of the gears is not significant to the dynamics of the system, in order to simplify the dynamics calculation, the shimmy is ignored in this study. Therefore, the transmission system of the hypoid gear pair is finally built as an 8-dof nonlinear dynamics model with the time-varying meshing stiffness and backlash. The generalized displacements of the system can be expressed as

q = {X_{1}, Y_{1}, Z_{1}, θ_{1 x}, X_{2}, Y_{2}, Z_{2}, θ_{2 y}}^{T}

(14)

where, $X_{i}, Y_{i}, Z_{i} (i = 1, 2)$ represent the linear displacement of the two gears along the X-axis, Y-axis, and Z-axis respectively, $θ_{1 x}, θ_{2 y}$ represent the angular displacement of the pinion along the X-axis and the angular displacement of the gear along the Y-axis.

Mathematical model

The contact force at the meshing point can be worked out based on the relative position parameters of the mating tooth surfaces. After establishing the positional relationship between the pinion and the gear shown in Figure 10, the relative displacement of the tooth surfaces along the direction of the normal force at the contact point can be represented as

\begin{matrix} X_{n} = X_{1} \cos δ_{1} \sin α_{n} - (X_{2} + r_{p 2} θ_{2 y}) \cos δ_{2} \sin α_{n} \\ + ((Y_{1} + r_{p 1} θ_{1 x}) (\cos δ_{1} \cos α_{n} \sin β_{1} + . . . \sin δ_{1} \sin α_{n}) \\ - Y_{2} (\cos δ_{2} \cos α_{n} \sin β_{2} + \sin δ_{2} \sin α_{n})) \\ - (Z_{1} (\sin δ_{1} \cos α_{n} \sin β_{1} + . . . \cos δ_{1} \sin α_{n}) \\ - Z_{2} (\sin δ_{2} \cos α_{n} \sin β_{2} + \cos δ_{2} \sin α_{n})) - e_{n} (t) \end{matrix}

(15)

e_{n} (t) = \sum_{i = 1}^{N} A_{ei} \cos (i Ω_{h} t + ϕ_{e})

(16)

where, $δ_{1}$ , $δ_{2}$ are the pitch angle of the pinion and the gear respectively; $α_{n}$ is the pressure angle on the normal plane; $β_{1}$ , $β_{2}$ are the helical angle at the midpoint on the pitch cone of the pinion and the gear respectively; $r_{p 1}$ , $r_{p 2}$ are the radius of the contact point on the pinion surface and the gear surface respectively; $e_{n} (t)$ is the static transmission error on the normal direction; $Ω_{h}$ is the frequency of the internal excitation from the gear system; $A_{ei}$ is the amplitude of the ith-order harmonic wave of the static transmission error; $ϕ_{e}$ is the initial phase of the static transmission error.

The normal force at the contact point and its component forces along the coordinate axes can be respectively expressed as

{\begin{matrix} F_{n} = \bar{k} (t) \bar{f} (X_{n}) + c_{m} {\overset{\cdot}{X}}_{n} \\ F_{x} = - F_{n} (\sin α_{n} \cos δ_{2} + \cos α_{n} \sin δ_{2} \sin β_{2}) \\ F_{y} = F_{n} (\sin α_{n} \sin δ_{2} - \cos α_{n} \cos δ_{2} \sin β_{2}) \\ F_{z} = - F_{n} \cos α_{n} \cos β_{2} \end{matrix}

(17)

\bar{k} (t) = k_{m} + \sum_{i = 1}^{N} A_{ki} \cos (i Ω_{h} t + ϕ_{k})

(18)

\bar{f} (X_{n} (t)) = {\begin{matrix} X_{n} (t) - b X_{n} (t) > b \\ 0 - b \leq X_{n} (t) \leq b \\ X_{n} (t) + b X_{n} (t) < - b \end{matrix}

(19)

where, ${\bar{k}}_{t}$ is the time-varying meshing stiffness of the gear pair; $k_{m}$ , $c_{m}$ are the mean meshing stiffness and mean meshing damping; $A_{ki}$ is the ith-order harmonic wave of the meshing stiffness; $ϕ_{k}$ is the initial phase of the meshing stiffness; $\bar{f} (X_{n})$ is the nonlinear function of the backlash of the gear pair; $b$ is half of the average meshing clearance. The differential equations of the motion of the gear pair can be represented as

{\begin{matrix} m_{1 x} {\overset{\cdot\cdot}{X}}_{1} + c_{x 1} {\overset{\cdot}{X}}_{1} + k_{x 1} X_{1} = - F_{x} \\ m_{1 y} {\overset{\cdot\cdot}{Y}}_{1} + c_{y 1} {\overset{\cdot}{Y}}_{1} + k_{y 1} Y_{1} = - F_{y} \\ m_{1 z} {\overset{\cdot\cdot}{Z}}_{1} + c_{z 1} {\overset{\cdot}{Z}}_{1} + k_{z 1} Z_{1} = - F_{z} \\ I_{1 x} {\overset{\cdot\cdot}{θ}}_{1 x} = T_{1} - F_{n} r_{b 1} \\ m_{2 x} {\overset{\cdot\cdot}{X}}_{2} + c_{x 2} {\overset{\cdot}{X}}_{2} + k_{x 2} X_{2} = F_{x} \\ m_{2 y} {\overset{\cdot\cdot}{Y}}_{2} + c_{y 2} {\overset{\cdot}{Y}}_{2} + k_{y 2} Y_{2} = F_{y} \\ m_{2 z} {\overset{\cdot\cdot}{Z}}_{2} + c_{z 2} {\overset{\cdot}{Z}}_{2} + k_{z 2} Z_{2} = F_{z} \\ I_{2 y} {\overset{\cdot\cdot}{θ}}_{2 y} = - T_{2} + F_{n} r_{b 2} \\ T_{1} = T_{1 m} + T_{1 v}, T_{2} = T_{2 m} \end{matrix}

(20)

where, $m_{ij}, I_{ij} (i = 1, 2; j = x, y, z)$ are the lumped mass and lumped rotational inertia of the gears; $c_{ij}, k_{ij} (i = 1, 2; j = x, y, z)$ are the coefficients of the damping and the stiffness of the gears along X-axis, Y-axis, and Z-axis; $T_{1}$ , $T_{2}$ are the drive torque and the load torque of the gear pair; $r_{b 1}$ , $r_{b 2}$ are the base radius of the pinion and the gear.

The drive torque $T_{1}$ on the pinion is consisted of the mean load torque $T_{1 m}$ and the time-varying torque $T_{1 v}$ , while the load torque $T_{2}$ on the gear is a constant torque $T_{2 m}$ . Obviously, equation (20) is a positive semi-definite nonlinear second-order differential equation set of an 8-dof system. Introducing the relative displacement $X_{n}$ as a new degree of freedom in equation (20), thereby the two torsional vibration displacements $θ_{1 x}, θ_{2 y}$ can be eliminated, and the freedom of the equation set can be reduced to seven. In order to simplify the equations, the following nine intermediate variables are introduced into the equation set.

{\begin{matrix} a_{1} = \cos δ_{1} \sin α_{n} \\ a_{2} = \cos δ_{2} \sin α_{n} \\ a_{3} = \cos δ_{1} \cos α_{n} \sin β_{1} + \sin δ_{1} \sin α_{n} \\ a_{4} = \cos δ_{2} \cos α_{n} \sin β_{2} + \sin δ_{2} \sin α_{n} \\ a_{5} = \sin δ_{1} \cos α_{n} \sin β_{1} + \cos δ_{1} \sin α_{n} \\ a_{6} = \sin δ_{2} \cos α_{n} \sin β_{2} + \cos δ_{2} \sin α_{n} \\ a_{7} = \sin δ_{2} \cos α_{n} \sin β_{2} + \cos δ_{2} \sin α_{n} = a_{6} \\ a_{8} = - \cos δ_{2} \cos α_{n} \sin β_{2} + \sin δ_{2} \sin α_{n} \\ a_{9} = \cos β_{2} \cos α_{n} \end{matrix}

(21)

The intermediate variables a₁–a₉ in equation (21) are all derived from the stress analysis in Figure 10. The derivation of damping and stiffness matrix and load vector using variables a₁–a₉ is shown in the appendix. Combining the two torsional vibration equations in equation (21), a new differential equation can be obtained.

\begin{matrix} - a_{1} m_{e} {\overset{\cdot\cdot}{X}}_{1} - a_{3} m_{e} {\overset{\cdot\cdot}{Y}}_{1} + a_{5} m_{e} {\overset{\cdot\cdot}{Z}}_{1} + a_{2} m_{e} {\overset{\cdot\cdot}{X}}_{2} + a_{4} m_{e} {\overset{\cdot\cdot}{Y}}_{2} \\ - a_{6} m_{e} {\overset{\cdot\cdot}{Z}}_{2} + m_{e} {\overset{\cdot\cdot}{X}}_{n} + b_{2} c_{m} {\overset{\cdot}{X}}_{n} + b_{2} \bar{k} (t) \bar{f} (X_{n}) = b_{1} F - {\overset{\cdot\cdot}{e}}_{n} (t) \end{matrix}

b_{1} = (\frac{a_{3}}{m_{1}} + \frac{a_{2}}{m_{2}}) m_{e}, b_{2} = (\frac{a_{3} r_{p 1}}{r_{b 1} m_{1}} + \frac{a_{2} r_{p 2}}{r_{b 2} m_{2}}) m_{e}

(22)

where, $m_{e}$ is the equivalent mass of gear pair, $m_{e} = m_{1} \cdot m_{2} / (m_{1} + m_{2})$ . $m_{1}$ , $m_{2}$ are the equivalent mass of the pinion and the gear respectively, $m_{1} = I_{1 x} / r_{b 1}^{2}$ , $m_{2} = I_{2 y} / r_{b 2}^{2}$ .

F = F_{1 m} + F_{1 v}

(23)

where, $F_{1 m}$ , $F_{1 v}$ are the amount of the average and variation of the circumferential force of the pinion respectively, $F_{1 m} = T_{1 m} / r_{b 1} = T_{2 m} / r_{b 2}$ , $F_{1 v} = a_{3} r_{p 1} m_{e} T_{1 v} / I_{1 x} = \sum_{i = 1}^{N} A_{Fi} \cos (i Ω_{e} t + ϕ_{e})$ . $Ω_{e}$ is the frequency of the external excitation, $A_{Fi}$ is the amplitude of the ith-order harmonic wave of the external load; $ϕ_{e}$ is the initial phase of the external load. Then the dimensionless form governing equations of gear motion are concluded as

{\begin{matrix} {\overset{\cdot\cdot}{x}}_{1} + 2 ζ_{x_{1}} {\overset{\cdot}{x}}_{1} - 2 ζ_{m x_{1}} a_{7} {\overset{\cdot}{x}}_{n} + k_{x_{1}} x_{1} - a_{7} k_{m x_{1}} (t) f (x_{n}) = 0 \\ {\overset{\cdot\cdot}{y}}_{1} + 2 ζ_{y_{1}} {\overset{\cdot}{y}}_{1} + 2 ζ_{m y_{1}} a_{8} {\overset{\cdot}{x}}_{n} + k_{y_{1}} y_{1} + a_{8} k_{m y_{1}} (t) f (x_{n}) = 0 \\ {\overset{\cdot\cdot}{z}}_{1} + 2 ζ_{z_{1}} {\overset{\cdot}{z}}_{1} - 2 ζ_{m z_{1}} a_{9} {\overset{\cdot}{x}}_{n} + k_{z_{1}} z_{1} - a_{9} k_{m y_{1}} (t) f (x_{n}) = 0 \\ {\overset{\cdot\cdot}{x}}_{2} + 2 ζ_{x_{2}} {\overset{\cdot}{x}}_{2} + 2 ζ_{m x_{2}} a_{7} {\overset{\cdot}{x}}_{n} + k_{x_{2}} x_{2} + a_{7} k_{m x_{2}} (t) f (x_{n}) = 0 \\ {\overset{\cdot\cdot}{y}}_{2} + 2 ζ_{y_{2}} {\overset{\cdot}{y}}_{2} - 2 ζ_{m y_{2}} a_{8} {\overset{\cdot}{x}}_{n} + k_{y_{2}} y_{2} - a_{8} k_{m y_{2}} (t) f (x_{n}) = 0 \\ {\overset{\cdot\cdot}{z}}_{2} + 2 ζ_{z_{2}} {\overset{\cdot}{z}}_{2} + 2 ζ_{m z_{2}} a_{9} {\overset{\cdot}{x}}_{n} + k_{z_{2}} z_{2} + a_{9} k_{m z_{2}} (t) f (x_{n}) = 0 \\ - a_{1} {\overset{\cdot\cdot}{x}}_{1} - a_{3} {\overset{\cdot\cdot}{y}}_{1} + a_{5} {\overset{\cdot\cdot}{z}}_{1} + a_{2} {\overset{\cdot\cdot}{x}}_{2} + a_{4} {\overset{\cdot\cdot}{y}}_{2} - a_{6} {\overset{\cdot\cdot}{z}}_{2} + {\overset{\cdot\cdot}{x}}_{n} + 2 ζ_{m} b_{2} {\overset{\cdot}{x}}_{n} \\ + k (t) b_{2} f (x_{n}) = f_{1 m} + f_{1 v} + f_{e} \end{matrix}

(24)

where, $f_{1 m}$ , $f_{1 v}$ are the amount of the average and variation of the dimensionless external force, $f_{e}$ is the dimensionless form of the force from the error excitation.

\begin{matrix} x_{j} = X_{j} / b, y_{j} = Y_{j} / b, z_{j} = Z_{j} / b, x_{n} = X_{n} / b, \\ ω_{n} = \sqrt{k_{m} / m_{e}}, ω_{ij} = \sqrt{k_{ij} / m_{ji}}, ω_{e} = Ω_{e} / ω_{n}, \\ ω_{h} = Ω_{h} / ω_{n}, ζ_{ij} = c_{ij} / (2 m_{ji} ω_{n}), ζ_{mij} = c_{m} / (2 m_{ji} ω_{n}), \\ ζ_{m} = c_{m} / (2 m_{e} ω_{n}), k_{ij} = ω_{ij}^{2} / ω_{n}^{2}, k_{mij} = k (t) / (m_{ji} ω_{n}^{2}), \\ f_{1 m} = b_{1} F_{1 m} / (m_{e} b ω_{n}^{2}), f_{1 v} = b_{1} F_{1 v} / (m_{e} b ω_{n}^{2}), τ = ω_{n} t \end{matrix}

\begin{matrix} k (t) = \frac{\bar{k} (t)}{k_{m}} = 1 + \sum_{i = 1}^{A_{kj}} \cos (i ω_{h} t + ϕ_{k}), \\ f_{e} = \sum_{i = 1}^{N} \frac{A_{ei}}{b} \cos (i ω_{h} τ + ϕ_{e}), i = x, y, z; j = 1, 2 . \end{matrix}

\bar{f} (x_{n}) = \frac{\bar{f} (X_{n})}{b} = {\begin{matrix} x_{n} - 1 x_{n} > 1 \\ 0 - 1 \leq x_{n} \leq 1 \\ x_{n} + 1 x_{n} < - 1 \end{matrix}

Simplifications, assumptions, and parameter settings

Before solving the nonlinear differential equation set, it is necessary to make appropriate simplifications and assumptions for some of the excitation terms in equation (24) to simplify the calculation process. Firstly, the driving torque on the pinion is assumed relatively stable, and its fluctuation is ignored. Therefore, $T_{1 v} = 0$ and $f_{1 v} = 0$ . For the dimensionless error excitation force, only the first order harmonic component is considered into calculation

f_{e} = γ ω_{h}^{2} \cos (ω_{h} τ + ϕ_{e})

(25)

where, $γ$ is the coefficient of the static transmission error, $γ = A_{e 1} / b$ .

For the time-varying meshing stiffness, according to the above fitting results by the Fourier series, the first to the eighth harmonic components are considered into calculation

k (τ) = 1 + \sum_{i = 1}^{8} ε_{i} \cos (ω_{h} τ + ϕ_{k_{i}})

(26)

where, $ε_{i}$ is the ith-order coefficient of the static transmission error, $ε_{i} = A_{ei} / k_{m}$ . The mean meshing stiffness and the coefficients of each order meshing stiffness are given in Table A2.

Since equation (24) is a strongly nonlinear second-order differential equation set, it is very difficult to solve it by analytical method. In this study, the fifth-order variable step size adaptive Runge-Kutta numerical integration method is used to solve it. Therefore, the seven second-order differential equations in equation (24) should be transferred to fourteen first-order differential equations. The state variables are introduced into the equation.

\begin{matrix} u (τ) = {u_{1}, u_{2}, u_{3}, u_{4}, u_{5}, u_{6}, u_{7}, u_{8}, u_{9}, u_{10}, u_{11}, u_{12}, u_{13}, {u_{14}}}^{T} \\ = {x_{1}, {\overset{\cdot}{x}}_{1}, y_{1}, {\overset{\cdot}{y}}_{1}, z_{1}, {\overset{\cdot}{z}}_{1}, x_{2}, {\overset{\cdot}{x}}_{2}, y_{2}, {\overset{\cdot}{y}}_{2}, z_{2}, {\overset{\cdot}{z}}_{2}, {x_{n}, {\overset{\cdot}{x}}_{n}}}^{T} \end{matrix}

(27)

Therefore, equation (24) can be transferred to a matrix form as

\overset{\cdot}{u} (τ) = Hu (τ) + f (u_{13}) A + (f_{1 m} + f_{1 v} + f_{e}) B

(28)

where, $H$ is the matrix of the linear terms; $A$ is the column vector of the coefficients of the nonlinear terms; $B$ is the column vector of the coefficients of the excitation terms.

In this research, the “ode45()” function in MATLAB R2016b software is employed to solve equation (24). The backlash of the gear pair is set as 72 $μ m$ , so $b = 36$ . The relevant dimensionless parameters in the nonlinear differential equation (24) are set as

\begin{matrix} f_{1 v} = 0, f_{1 m} = 0.5, Ω_{e} = 0, ζ_{ij} = 0.01, \\ ζ_{mij} = 0.02, (i = x, y, z; j = 1, 2), ε = 0.8, \\ ζ_{m} = 0.05, ϕ_{e} = 0, ϕ_{k} = 0 . \end{matrix}

The meshing stiffness under 500, 1000, and 2000 Nm are given in Appendix equation (A1), and the initial condition of the differential equation set is set as

\begin{array}{l} x_{n} (0) = 0.1, {\dot{x}}_{n} (0) = 0.1, x_{i} (0) = y_{i} (0) = z_{i} (0) = 0, \\ {\dot{x}}_{i} (0) = {\dot{y}}_{i} (0) = {\dot{z}}_{i} (0) = 0, (i = 1, 2) \end{array}

Comparison of the effects of different meshing stiffness on the dynamics

In order to investigate the effects of different meshing stiffness on the dynamics characteristics of the hypoid gear system, the above three expressions of dimensionless stiffness are introduced into the dynamics model. After completing the simulation, the time history curve, phase diagram, Poincaré section, and the bifurcation diagram of the system are used to analysis the influence. When the normalized excitation frequency $ω_{h}$ is set as 0.23, the time history map, trajectory on phase plane, and Poincaré plot of the normalized relative displacement $x_{n}$ under different load conditions are illustrated in Figure 11.

Figure 11.

As can be seen from Figure 11, the vibration amplitude of the normalized relative displacement $x_{n}$ under different loads is significantly different. Moreover, the vibration of $x_{n}$ under 500 N m load condition is a two-cycle motion, and that under 2000 N m load condition is a four-cycle motion, while that under 1000 N m load condition is a chaotic motion. Apparently, at the same excitation frequency, there are big differences in the steady-state dynamics features of the hypoid gear system corresponding to different meshing stiffness.

When $ω_{h}$ is set as 0.90, the time history map, trajectory on phase plane, and Poincaré plot of $x_{n}$ are illustrated in Figure 12. It can be seen that the vibration of $x_{n}$ under 500 N m load condition is a two-cycle motion, while those under 1000 and 2000 N m load condition are both single-cycle motion. Consequently, at the same excitation frequency, the vibration characteristics of the hypoid gear system under different loads are significantly different. Under the same load condition, the vibration features of the gear system change with the excitation frequency $ω_{h}$ . The influence of the time-varying meshing stiffness and the load on the dynamics cannot be neglected in the analysis of dynamics of hypoid gear system. In addition, the effects of the excitation frequency on the vibration features of gear system should be analyzed.

Figure 12.

The normalized displacement and normalized velocity of the pinion along X-axis, Y-axis, and Z-axis under different loads are shown in Figure 13. It can be seen that under the three load conditions, the vibration amplitudes along the offset direction are all the largest, followed by those along the axial direction of the pinion, while those along the gear axial direction are minimum. This rule also applies to the comparison of the normalized velocities of the pinion along the three axes. Since the hypoid gears are supported by the tapered roller bearings, these vibration features are consistent to the actual motion state of the gear system.

Figure 13.

In order to study the influence of the excitation frequency on the vibration characteristics of the system, the bifurcation diagram of $x_{n}$ relative to $ω_{h}$ under different load conditions are illustrated in Figure 13. For the sake of demonstrating the effects of the fitting accuracy of the meshing stiffness on the dynamics characteristics of the gear system, the bifurcation diagrams of $x_{n}$ relative to $ω_{h}$ with different meshing stiffness fitted by second-order Fourier series are illustrated in Figure 14.

Figure 14.

Bifurcation diagram of x_n relative to ω_h with different meshing stiffness fitted by eighth-order Fourier series: (a) 500 N m load condition, (b) 1000 N m load condition, and (c) 2000 N m load condition.

It can be seen from Figure 14 that as the excitation frequency changes, the local characteristics of the normalized relative displacement $x_{n}$ of the hypoid gear pair under different loads are significantly different. Firstly, the vibration amplitudes of $x_{n}$ under different loads are different, which can be clearly seen within the range of $ω_{h}$ from 0.6 to 0.8. Secondly, the vibration state of $x_{n}$ under different loads are also apparently different, which can be found within the range of $ω_{h}$ from 1.4 to 1.8. However, there is no obvious difference in the global vibration characteristics of $x_{n}$ under the three loads, which start as a single-cycle motion and undergo the multi-period, quasi-period, chaotic motion state, and finally return to a single-cycle motion.

By comparing the bifurcation diagrams under the same load condition in Figures 14 and 15, it can be found that the influence of the fitting accuracy of the time-varying stiffness on the motion of the gear system cannot be ignored. For instance, there are visible differences between the local vibration state of $x_{n}$ within the range of $ω_{h}$ from 1.4 to 1.6 under 500 N m load condition in Figures 14 and 15. Nevertheless, the local and global vibration states of $x_{n}$ are not significantly different under 1000 and 2000 N m load condition in Figures 14 and 15. Therefore, it can be concluded that the motion state of the hypoid gear system under light load condition is more sensitive to the fitting accuracy of the meshing stiffness.

Figure 15.

Bifurcation diagram of x_n relative to ω_h with different meshing stiffness fitted by second-order Fourier series: (a) 500 N m load condition, (b) 1000 N m load condition, and (c) 2000 N m load condition.

Conclusion

A detailed mathematical model has been given in this study for computing the meshing stiffness of a hypoid gear pair by measuring the contact forces and the static and dynamic transmission error. The reliability and accuracy of the calculation method are verified by using the same method to compute the stiffness of a spur gear pair. A generalized 8-dof nonlinear dynamics model of the gear drive is established considering the time-varying stiffness, backlash, dynamic transmission error, elastic supports, and neglecting the shimmy of the gears. Meshing stiffness obtained by the above analysis is fitted by an eighth-order Fourier series and imported into the dynamics model. The influence of the meshing stiffness under different loads on the dynamics response of the system is analyzed and compared. The effects of the excitation frequency on the motion characteristics of the gear pair are also analyzed. The conclusion can be summarized as follows:

a. Meshing stiffness of the hypoid gear pair changes periodically, and the period is same to the changing period of the contact forces and the transmission error, which equals to the meshing period of the gear pair. Unlike the variation of the meshing stiffness of spur gear pair, the meshing stiffness of hypoid gear pair will not change suddenly, thereby the gears can smoothly contact and the meshing impact can be reduced. The meshing stiffness increases with the increase of the load, mainly because the increase of the contact forces enhances the deformations of the tooth profile and the tooth surface, so that the average value of the meshing stiffness increases. The load has a major effect on the amplitude of meshing stiffness, and the consideration of load is essential in the dynamics analysis of hypoid gears.

b. The accuracy and goodness of the meshing stiffness curves fitted by the eighth-order Fourier series are better than those fitted by the second-order Fourier series. The fitting of the low-order function for time-varying meshing parameters may be not accurate enough, which may cause errors in the dynamic analysis results. In the future studies on the dynamics of hypoid gears, the time-varying meshing stiffness should be represented as accurately as possible.

c. Under different load conditions, the vibration amplitudes along the offset direction are all the largest, followed by that along the axial direction of the pinion, while those along the gear axial direction are minimum. This rule also applies to the comparison of the normalized velocities of the pinion along the three axes. Since the hypoid gears are supported by the tapered roller bearings, these vibration features are consistent to the actual motion state of the gear system. During the assembly process of the hypoid gear pair, it is necessary to restrain the displacement of the gears along the offset distance.

d. As the excitation frequency increases, the hypoid gear pair studied in this paper experiences singe-cycle, multi-cycle, quasi-periodic, and chaotic motion state and finally returns to a single-cycle motion. At the same excitation frequency, there are big differences in the steady-state dynamics features of the hypoid gear system corresponding to different meshing stiffness. Under the same load condition, the vibration features of the gear system change with the excitation frequency. Hence, the influence of the time-varying meshing stiffness and the load on the dynamics cannot be neglected in the analysis of dynamics of hypoid gear system. The effects of the excitation frequency on the vibration features of gear system should be analyzed, and the influence of the fitting accuracy of the meshing stiffness on the motion of the gear system also should be considered. Besides, the conclusion can be drawn from comparative analysis that the motion state of the hypoid gear system under light load condition is more sensitive to the fitting accuracy of the meshing stiffness.

e. Compared with previous and similar hypoid gear system dynamics models, the high-order time-varying meshing stiffness and transmission error under different loads are introduced into the dynamics model instead of using the mean value and the first harmonics. Therefore, the effect of load variation on the behavior of gear system dynamics can be found in detail. Furthermore, the dynamics behavior of hypoid gear system considering the variation of the meshing damping, frequency, and tooth clearance under different loads is also discussed. However, in order to simplify the calculation process, the 12 DOF dynamic system is simplified to 8 DOF in this paper, which may lead to slight deviation between the calculated results and the real state. Moreover, there are assembly errors between the meshing hypoid gears, and the influence of misalignments is not considered in this paper, which is also an important content of follow-up research.

Footnotes

Appendix

(A1)

{\begin{matrix} k_{m_500} = 1 - 0.21039 \cos (ω τ) - 0.28979 \sin (ω τ) + 0.030914 \cos (2 ω τ) + 0.005569 \sin (2 ω τ) \\ - 0.0393 \cos (3 ω τ) + 0.036809 \sin (3 ω τ) + 0.010475 \cos (4 ω τ) + 0.001201 \sin (4 ω τ) \\ + 0.012192 \cos (5 ω τ) + 0.007266 \sin (5 ω τ) + 0.00572 \cos (6 ω τ) - 0.00555 \sin (6 ω τ) \\ + 0.004254 \cos (7 ω τ) - 0.00337 \sin (7 ω τ) - 0.00278 \cos (8 ω τ) + 0.00041 \sin (8 ω τ) \\ k_{m_1000} = 1 - 0.08226 \cos (ω τ) - 0.24247 \sin (ω τ) + 0.064409 \cos (2 ω τ) - 0.04905 \sin (2 ω τ) \\ + 0.016165 \cos (3 ω τ) + 0.024721 \sin (3 ω τ) - 0.00719 \cos (4 ω τ) - 0.00176 \sin (4 ω τ) \\ + 0.007085 \cos (5 ω τ) - 0.00445 \sin (5 ω τ) + 0.0042 \cos (6 ω τ) + 0.006798 \sin (6 ω τ) \\ - 0.00574 \cos (7 ω τ) + 0.003942 \sin (7 ω τ) - 0.00293 \cos (8 ω τ) - 0.00348 \sin (8 ω τ) \\ k_{m_2000} = 1 - 0.10954 \cos (ω τ) - 0.11944 \sin (ω τ) + 0.057602 \cos (2 ω τ) - 0.03099 \sin (2 ω τ) \\ + 0.01652 \cos (3 ω τ) + 0.021336 \sin (3 ω τ) - 0.01339 \cos (4 ω τ) + 0.010087 \sin (4 ω τ) \\ - 0.00412 \cos (5 ω τ) - 0.00922 \sin (5 ω τ) + 0.005518 \cos (6 ω τ) - 0.00267 \sin (6 ω τ) \\ + 0.002282 \cos (7 ω τ) + 0.003951 \sin (7 ω τ) - 0.00342 \cos (8 ω τ) + 0.002184 \sin (8 ω τ) \end{matrix}

The dynamics model of the gear drive described as equation (20) can be represented as Matrix form as

(A2)

M {\overset{\cdot\cdot}{q}} + C {\overset{\cdot}{q}} + K {q} = {F}

where M , C , K are the mass, damping, and stiffness matrix of the gear dynamics system respectively, F is the load vector of the system, ${\overset{\cdot\cdot}{q}}$ , ${\overset{\cdot}{q}}$ , ${q}$ are the generalized acceleration, velocity, and displacement vector respectively.

(A3)

M = diag [m_{1 x}, m_{1 y}, m_{1 z}, I_{1 x}, m_{2 x}, m_{2 y}, m_{2 z}, I_{2 y}]

(A4)

C = [\begin{matrix} c_{x 1} - a_{1} a_{7} c_{m} & - a_{3} a_{7} c_{m} & a_{5} a_{7} c_{m} & - a_{3} a_{7} r_{p 1} c_{m} \\ a_{1} a_{8} c_{m} & c_{y 1} + a_{3} a_{8} c_{m} & - a_{5} a_{8} c_{m} & a_{3} a_{8} r_{p 1} c_{m} \\ - a_{1} a_{9} c_{m} & - a_{3} a_{9} c_{m} & c_{z 1} + a_{5} a_{9} c_{m} & - a_{3} a_{9} r_{p 1} c_{m} \\ - a_{1} r_{1} c_{m} & - a_{3} r_{1} c_{m} & a_{5} r_{1} c_{m} & a_{3} r_{1} r_{p 1} c_{m} \\ a_{1} a_{7} c_{m} & a_{3} a_{7} c_{m} & - a_{5} a_{7} c_{m} & a_{3} a_{7} r_{p 1} c_{m} \\ - a_{1} a_{8} c_{m} & - a_{3} a_{8} c_{m} & a_{5} a_{8} c_{m} & - a_{3} a_{8} r_{p 1} c_{m} \\ a_{1} a_{9} c_{m} & a_{3} a_{9} c_{m} & - a_{5} a_{9} c_{m} & a_{3} a_{9} r_{p 1} c_{m} \\ a_{1} a_{2} c_{m} & a_{3} r_{2} c_{m} & - a_{5} a_{2} c_{m} & a_{3} r_{2} r_{p 1} c_{m} \end{matrix} \begin{matrix} a_{2} a_{7} c_{m} & a_{4} a_{7} c_{m} & - a_{6} a_{7} c_{m} & a_{2} a_{7} r_{p 2} c_{m} \\ - a_{2} a_{8} c_{m} & - a_{4} a_{8} c_{m} & a_{6} a_{8} c_{m} & - a_{2} a_{8} r_{p 2} c_{m} \\ a_{2} a_{9} c_{m} & a_{4} a_{9} c_{m} & - a_{6} a_{9} c_{m} & a_{2} a_{9} r_{p 2} c_{m} \\ a_{2} r_{1} c_{m} & a_{4} r_{1} c_{m} & - a_{6} r_{1} c_{m} & a_{2} r_{1} r_{p 2} c_{m} \\ c_{x 2} - a_{2} a_{7} c_{m} & - a_{4} a_{7} c_{m} & a_{6} a_{7} c_{m} & - a_{2} a_{7} r_{p 2} c_{m} \\ a_{2} a_{8} c_{m} & c_{y 2} + a_{4} a_{8} c_{m} & - a_{6} a_{8} c_{m} & a_{2} a_{8} r_{p 2} c_{m} \\ - a_{2} a_{9} c_{m} & - a_{4} a_{9} c_{m} & c_{z 2} + a_{6} a_{9} c_{m} & - a_{2} a_{9} r_{p 2} c_{m} \\ - a_{2} r_{2} c_{m} & - a_{4} r_{2} c_{m} & a_{6} r_{2} c_{m} & - a_{2} r_{2} r_{p 2} c_{m} \end{matrix}]

(A5)

K = [\begin{matrix} k_{x 1} - a_{1} a_{7} c_{m} & - a_{3} a_{7} k_{m} & a_{5} a_{7} k_{m} & - a_{3} a_{7} r_{p 1} k_{m} \\ a_{1} a_{8} k_{m} & c_{y 1} + a_{3} a_{8} k_{m} & - a_{5} a_{8} k_{m} & a_{3} a_{8} r_{p 1} k_{m} \\ - a_{1} a_{9} k_{m} & - a_{3} a_{9} k_{m} & c_{z 1} + a_{5} a_{9} k_{m} & - a_{3} a_{9} r_{p 1} k_{m} \\ - a_{1} r_{1} k_{m} & - a_{3} r_{1} k_{m} & a_{5} r_{1} k_{m} & a_{3} r_{1} r_{p 1} k_{m} \\ a_{1} a_{7} k_{m} & a_{3} a_{7} k_{m} & - a_{5} a_{7} k_{m} & a_{3} a_{7} r_{p 1} k_{m} \\ - a_{1} a_{8} k_{m} & - a_{3} a_{8} k_{m} & a_{5} a_{8} k_{m} & - a_{3} a_{8} r_{p 1} k_{m} \\ a_{1} a_{9} k_{m} & a_{3} a_{9} k_{m} & - a_{5} a_{9} k_{m} & a_{3} a_{9} r_{p 1} k_{m} \\ a_{1} a_{2} k_{m} & a_{3} r_{2} k_{m} & - a_{5} a_{2} k_{m} & a_{3} r_{2} r_{p 1} k_{m} \end{matrix} \begin{matrix} a_{2} a_{7} k_{m} & a_{4} a_{7} k_{m} & - a_{6} a_{7} k_{m} & a_{2} a_{7} r_{p 2} k_{m} \\ - a_{2} a_{8} k_{m} & - a_{4} a_{8} k_{m} & a_{6} a_{8} k_{m} & - a_{2} a_{8} r_{p 2} k_{m} \\ a_{2} a_{9} k_{m} & a_{4} a_{9} k_{m} & - a_{6} a_{9} k_{m} & a_{2} a_{9} r_{p 2} k_{m} \\ a_{2} r_{1} k_{m} & a_{4} r_{1} k_{m} & - a_{6} r_{1} k_{m} & a_{2} r_{1} r_{p 2} k_{m} \\ k_{x 2} - a_{2} a_{7} k_{m} & - a_{4} a_{7} k_{m} & a_{6} a_{7} k_{m} & - a_{2} a_{7} r_{p 2} k_{m} \\ a_{2} a_{8} k_{m} & k_{y 2} + a_{4} a_{8} k_{m} & - a_{6} a_{8} k_{m} & a_{2} a_{8} r_{p 2} k_{m} \\ - a_{2} a_{9} k_{m} & - a_{4} a_{9} k_{m} & k_{z 2} + a_{6} a_{9} k_{m} & - a_{2} a_{9} r_{p 2} k_{m} \\ - a_{2} r_{2} k_{m} & - a_{4} r_{2} k_{m} & a_{6} r_{2} k_{m} & - a_{2} r_{2} r_{p 2} k_{m} \end{matrix}]

(A6)

{F} = [\begin{matrix} - a_{7} c_{m} {\overset{\cdot}{e}}_{n} (t) - a_{7} k_{m} e_{n} (t) \\ a_{8} c_{m} {\overset{\cdot}{e}}_{n} (t) + a_{8} k_{m} e_{n} (t) \\ - a_{9} c_{m} {\overset{\cdot}{e}}_{n} (t) - a_{9} k_{m} e_{n} (t) \\ T_{1} - r_{1} c_{m} {\overset{\cdot}{e}}_{n} (t) - r_{1} k_{m} e_{n} (t) \\ a_{7} c_{m} {\overset{\cdot}{e}}_{n} (t) + a_{7} k_{m} e_{n} (t) \\ - a_{8} c_{m} {\overset{\cdot}{e}}_{n} (t) - a_{8} k_{m} e_{n} (t) \\ a_{9} c_{m} {\overset{\cdot}{e}}_{n} (t) + a_{9} k_{m} e_{n} (t) \\ - T_{2} + r_{2} c_{m} {\overset{\cdot}{e}}_{n} (t) - r_{2} k_{m} e_{n} (t) \end{matrix}]

(A7)

{\overset{\cdot}{q}} = {d {q} / dt = [{\overset{\cdot}{X}}_{1}, {\overset{\cdot}{Y}}_{1}, {\overset{\cdot}{Z}}_{1}, {\overset{\cdot}{θ}}_{1 x}, {\overset{\cdot}{X}}_{2}, {\overset{\cdot}{Y}}_{2}, {\overset{\cdot}{Z}}_{2}, {\overset{\cdot}{θ}}_{2 y}]}^{T}

(A8)

{\overset{\cdot\cdot}{q}} = {d^{2} {q} / d t^{2} = [{\overset{\cdot\cdot}{X}}_{1}, {\overset{\cdot\cdot}{Y}}_{1}, {\overset{\cdot\cdot}{Z}}_{1}, {\overset{\cdot\cdot}{θ}}_{1 x}, {\overset{\cdot\cdot}{X}}_{2}, {\overset{\cdot\cdot}{Y}}_{2}, {\overset{\cdot\cdot}{Z}}_{2}, {\overset{\cdot\cdot}{θ}}_{2 y}]}^{T}

Handling Editor: Sharmili Pandian

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors would like to thank the Universities Natural Science Research Project of Anhui Province of China (Grant No. 2024AH050475), the National Natural Science Foundation of China (Grant No. 51875150), the Aeronautical Science Foundation of China (Grant No. 2018ZEP4006), and Anhui Agricultural University youth Science Fund project (Grant No. rc412011) for the support given to this research.

ORCID iD

Xiaole Wang

References

Litvin

Fuentes

Gear geometry and applied theory. 2nd ed. New York, NY: Cambridge University Press, 2004.

Litvin

Fuentes

Hayasaka

Design, manufacture, stress analysis, and experimental tests of low-noise high endurance spiral bevel gears. Mech Mach Theory 2006; 41: 83–118.

Wang

Fong

ZH.

Fourth-order kinematic synthesis for face-milling spiral bevel gears with modified radial motion (MRM) correction. J Mech Des 2006; 128: 457–467.

Simon

Machine-tool settings to reduce the sensitivity of spiral bevel gears to tooth errors and misalignments. J Mech Des 2008; 130: 082603.

Artoni

Bracci

Gabiccini

, et al. Optimization of the loaded contact pattern in hypoid gears by automatic topography modification. J Mech Des 2009; 131: 011008.

Gonzalez-Perez

Fuentes

Hayasaka

Analytical determination of basic machine-tool settings for generation of spiral bevel gears from blank data. J Mech Des 2010; 132: 101002.

Simon

Optimization of face-hobbed hypoid gears. Mech Mach Theory 2014; 77: 164–181.

Wang

Zhang

Wan

Global synthesis for face milled spiral bevel gears with zero transmission errors. J Mech Des 2016; 138: 033302.

Temirkhan

Amrin

Spitas

, et al. Convergence and accuracy problems of the conventional TCA model–critical analysis and novel solution for crowned spur gears. Proc IMechE, Part C: J Mechanical Engineering Science 2024; 238: 724–736.

10.

Litvin

Chen

, et al. Application of finite element analysis for determination of load share, real contact ratio, precision of motion, and stress analysis. J Mech Des 1996; 118: 561–567.

11.

Litvin

Fuentes

Fan

, et al. Computerized design, simulation of meshing, and contact and stress analysis of face-milled formate generated spiral bevel gears. Mech Mach Theory 2002; 37: 441–459.

12.

Vogel

Griewank

Bär

Direct gear tooth contact analysis for hypoid bevel gears. Comput Methods Appl Mech Eng 2002; 191: 3965–3982.

13.

Temirkhan

Tariq

Spitas

, et al. Parametric design of straight bevel gears based on a new tooth contact analysis model. Arch Appl Mech 2023; 93: 4181–4196.

14.

Fan

High-order tooth flank form error correction for face-milled spiral bevel and hypoid gears. J Mech Des 2008; 130: 072601.

15.

Fan

Tooth surface error correction for face-hobbed hypoid gears. J Mech Des 2010; 132: 011004.

16.

Sakaridis

Spitas

Non-linear modeling of gear drive dynamics incorporating intermittent tooth contact analysis and tooth eigenvibrations. Mech Mach Theory 2019; 136: 307–333.

17.

Cheng

Dynamics of high-speed hypoid and bevel geared rotor systems. PhD Dissertation, Ohio State University, USA, 2000.

18.

Wang

Gear mesh characteristics and dynamics of hypoid gear rotor system. PhD Dissertation, University of Alabama, USA, 2002.

19.

Peng

Coupled multi-body dynamic and vibration analysis of hypoid and bevel geared rotor system. PhD Dissertation, University of Cincinnati, USA, 2010.

20.

Spitas

Coupled multi-DOF dynamic contact analysis model for the simulation of intermittent gear tooth contacts, impacts and rattling considering backlash and variable torque. Proc IMechE, Part C: J Mechanical Engineering Science 2016; 230: 1022–1047.

21.

Hua

Lim

Peng

, et al. Dynamic analysis of spiral bevel geared rotor systems applying finite elements and enhanced lumped parameters. Int J Automot Technol 2012; 13: 97–107.

22.

Shi

Nonlinear time-varying dynamic modeling of vehicle driveline systems with emphasis on hypoid gear excitation and response. PhD Dissertation, University of Cincinnati, USA, 2016.

23.

Yang

Gao

Deng

, et al. Nonlinear dynamics of spiral bevel gears and hypoid gear drives. J Aerosp Power 2004; 19: 54–57.

24.

Wang

Fang

Xie

, et al. Vibration and stability analysis of a spinning shaft with arbitrary boundaries subjected to partial load. J Vib Control 2024; 30: 1378–1390.

25.

Wang

Zhao

Chen

, et al. Vibration and stability analysis of a spinning hollow functionally graded beam with exponentially varying cross-section. Int J Struct Stab Dyn 2024; 24: 2450061.

26.

Wang

Lim

Dynamics of a hypoid gear pair considering the effects of time-varying mesh parameters and backlash nonlinearity. J Sound Vib 2007; 308: 302–329.

27.

Wang

Lim

TC.

Effect of tooth mesh stiffness asymmetric nonlinearity for drive and coast sides on hypoid gear dynamics. J Sound Vib 2009; 319: 885–903.

28.

Zheng

Influence of asymmetric mesh stiffness on dynamics of spiral bevel gear transmission system. Math Probl Eng 2010; 2010: 124–148.

29.

Song

Zhu

Lim

, et al. Parametric analysis of gear mesh and dynamic response of loaded helical beveloid transmission with small shaft angle. J Mech Des 2012; 134: 084501.

30.

Yavuz

Saribay

Cigeroglu

Nonlinear dynamic analysis of a spiral bevel geared system. In: DiMaio

Castellini

(eds) Rotating machinery, hybrid test methods, vibro-acoustics & laser vibrometry, volume 8. Conference Proceedings of the Society for Experimental Mechanics Series. Cham: Springer, 2017, pp.31–40.

31.

Kiyohiko

The meshing test on helical gears under load transmission: 1st report, the appropriate formula for deflections of gear tooth. JSME Int J 1972; 38: 896–904.

32.

Kuang

Yang

. An estimate of mesh stiffness and load sharing ratio of a spur gear pair. In: Proceedings of the ASME journal of 12th international power transmission and gearing conference, Scottsdale, AZ, USA, 13–16 September 1992, pp.1–9. New York: ASME.

33.

Gosselin

Gagnon

Cloutier

Accurate tooth stiffness of spiral bevel gear teeth by the finite strip method. J Mech Des 1998; 120: 599–605.

34.

Mennem

RC.

Parametrically excited vibrations in spiral bevel geared systems. PhD Dissertation, Purdue University, USA, 2004.

35.

Mohammadpour

Theodossiades

Rahnejat

Multiphysics investigations on the dynamics of differential hypoid gears. J Vib Acoust 2014; 136: 041007.

36.

Karagiannis

Theodossiades

Rahnejat

On the dynamics of lubricated hypoid gears. Mech Mach Theory 2012; 48: 94–120.

37.

Fang

ZD.

Dynamic load calculation of hypoid gear. Automobile Eng 1994; 16: 92–97. (in Chinese)

38.

Yang

Shi

Zhang

Dynamic analysis of spiral bevel and hypoid gears with high-order transmission errors. J Sound Vib 2018; 417: 149–164.