A general contact stiffness model for elastic bodies and its application in time-varying mesh stiffness of gear drive

Abstract

As a critical element of time-varying mesh stiffness (TVMS), contact stiffness of a gear drive has been defined based on simplified Hertzian contact stiffness or semi-empirical nonlinear Hertzian contact stiffness in previous works. This study proposes a general contact stiffness model for elastic bodies through piecewise linear interpolation of contact pressure. The TVMS of a spur gear drive is determined through potential energy method and proposed contact stiffness model verified by Hertzian contact theory and finite-element method. Then, the influence of applied load on contact stiffness is studied, and the differences among proposed contact stiffness, simplified Hertzian contact stiffness, and nonlinear Hertzian contact stiffness are analyzed. Results show that contact stiffness increases with the applied load, and the TVMS based on the proposed contact stiffness model is the smallest among the three contact stiffness models. Effects of tooth width and input torque on the TVMS are further discussed. The TVMS becomes bigger with increased tooth width and input torque, but the increase rate decreases as tooth width or input torque increases. These findings indicate that reasonable matching of design parameters is beneficial for increasing load capacity and optimizing the dynamic performance of gear systems.

Keywords

Contact stiffness time-varying mesh stiffness potential energy method elastic contact gear drive

1. Introduction

As a periodic and time-varying internal excitation of a gear system, time-varying mesh stiffness (TVMS) significantly affects vibration and noise [1 –6]. The TVMS of the gear system consists of the following types of stiffness: contact, bending, axial compression, fillet foundation, and shear stiffness [7]. However, in previous studies, the contact stiffness of the gear drive was defined on the basis of simplified or semi-empirical contact stiffness formulas due to complex contact mechanical characteristics between two elastic bodies. In addition, the contact status is also the key factor affecting the gear transmission performance and dynamic characteristics [8]. To calculate the TVMS accurately and establish a more reliable dynamic model, we have to investigate the general contact stiffness model for the gear system.

Experimental [9 –11], finite element (FE) [1,12 –22], and analytical methods [5,12,15 –34] have been used in calculating the TVMS. As the experimental method is the most direct and accurate approach, techniques such as photoelasticity [9], strain gauge [10], and digital image correlation [11] have been widely applied. However, these experimental methods are costly and cannot be widely applied. In contrast to the experimental methods, the TVMS of the gear system has drawn increasing attention with the advent of efficient numerical methods. Furthermore, two-dimensional (2D) FE [1,12,13,15 –17,22] and three-dimensional (3D) FE methods [14,18 –21] have been developed to verify the TVMS. The FE method has lower cost than the experimental method and can be applied to determine the TVMS with tooth profile modification [17,20] or crack [14 –16,18,19]. Although the FE method is widely used in calculating the TVMS, the FE method is time-consuming and complex because it needs extremely high mesh refinement to ensure its computational accuracy.

The analytical method is simpler and more effective than the experimental and FE methods, and many researchers have focused on developing a highly accurate analytical method to calculate the TVMS. Yang and Lin [23] first calculated the TVMS by considering simplified Hertzian contact stiffness (SHCS), bending stiffness, and axial stiffness using potential energy method. Subsequently, the potential energy method was improved by considering shear stiffness [35] and fillet foundation stiffness [12]. However, in the aforementioned studies, the nonuniform cantilever beam model of gear tooth is based on the base circle and the calculation error appears when the tooth number is larger or smaller. To eliminate this error, Ma et al. [16] presented a TVMS model that established a nonuniform cantilever beam of gear tooth based on the root circle and demonstrated that the model is more accurate than cantilever beam model of gear tooth based on the base circle. With the development of the potential energy method, researchers have investigated the TVMS of cracked [5,15,16,18,19,25 –28,30], spalling defect [29], tooth-profile modified [17,20], and helical gears [30 –33].

From the aforementioned studies, significant efforts have been exerted to accurately determine the TVMS for the gear system. However, as an indispensable component of the TVMS, the contact stiffness has rarely been mentioned. Since Yang and Sun [36] first derived the simplified formula of Hertzian contact stiffness in 1985, this contact stiffness model has been popularly applied to the TVMS modeling of gear pairs [26,27,37 –39]. Considering the nonlinear characteristics of the contacting tooth surfaces, Cornell [40] adopted an accurate formula in which the Hertzian contact stiffness is nonlinear. This contact stiffness model has also been applied in many research works, such as [20,22,41]. However, the nonlinear Hertzian contact stiffness (NHCS) is a semi-empirical formula, and its engineering application faces considerable limitations. Recently, Sainsot and Velex [42] presented a thorough review and analysis of the classic contact deflection formulae in gear simulations and proposed a general formula for contact stiffness in gear teeth. Based on contact mechanics and gear engagement principle, the analytical formula of local contact stiffness is determined [43]. The above research work is very representative, but their work builds on the Hertz contact model and also introduces some empirical formulas, which may bring some errors.

With the widespread application of high-performance gears, neither the simplified nor the semi-empirical contact stiffness formulas can meet the calculation requirements of high-precision TVMS. The contact stiffness significantly affects the TVMS, which is the key factor that influences the load capacity of the contacting tooth surfaces and the dynamic performance of the gear pair. Therefore, a general contact stiffness model for the gear system and practical design solutions to improve the load capacity and dynamic performance of gear system have to be investigated. This study focuses on the calculation method of contact stiffness for elastic bodies and its application to the TVMS modeling of the gear system. The influence of design parameters on gear TVMS is also examined. A numerical contact stiffness model with high computational accuracy for elastic bodies is proposed by piecewise linear interpolating the contact pressure. The TVMS of spur gears is determined in accordance with the stiffness models of axial compression, bending, fillet foundation, and shear, which are derived by the potential energy method and proposed contact stiffness model. The proposed contact model is verified by Hertzian contact theory and FE method, and the differences among proposed contact stiffness, SHCS and NHCS are analyzed. Then, the TVMS of a spur gear drive with three contact stiffness models are studied. The effects of the main design parameters on the TVMS of the gear drive are subsequently discussed.

2. General contact stiffness model

2.1. Contact model

As described in Figure 1, according to elastic theory, when an elastic body is subjected to a load Q at the point (τ, 0), the normal displacement at point (x, 0) is as follows [44]:

y (x) = \frac{2 (1 - υ^{2}) Q}{π E} \ln | \frac{d_{0}}{x - τ} | - \frac{(1 + υ) Q}{π E},

(1)

Figure 1.

Normal displacement under load Q.

where E is Young’s modulus, υ is Poisson’s ratio, and d₀ is a constant depth, where the normal displacement of each point is zero.

To eliminate the influence of d₀, a certain point (x_r, 0) is selected as the basic point. Then, the relative normal displacement between points (x, 0) and (x_r, 0) is given as follows:

d (x) = y (x) - y (x_{r}) = \frac{2 (1 - υ^{2}) Q}{π E} \ln | \frac{x_{r} - τ}{x - τ} | .

(2)

For the contact of two elastic bodies, the contact pressure is continuously distributed. Figure 2 describes the contact pressure distribution in the contact area [τ_a, τ_b], where the contact pressure is continuously distributed, τ_a and τ_b are the edges of the contact area. Therefore, for the contact of two elastic bodies, the influence of distributed loads should be considered.

Figure 2.

Pressure distribution between τ = τ_a and τ = τ_b.

Equation (2) can be rewritten as follows:

d (x) = \frac{2 (1 - υ^{2})}{π E} \int_{τ_{a}}^{τ_{b}} \ln | \frac{x_{r} - τ}{x - τ} | P (τ) d τ .

(3)

As shown in Figure 3, when two elastic bodies contact at point O, the initial gap between points a₁(x, y₁) and a₂(x, y₂) is h(x), and the contact zone is Ω_c under load Q.

Figure 3.

Contact model of two elastic bodies.

The contact equation is expressed as follows:

y_{1} (x) + y_{2} (x) + h (x) - δ {\begin{matrix} = 0 & x \in Ω_{c} \\ > 0 & x \notin Ω_{c} \end{matrix},

(4)

where δ is the elastic approach. Taking the displacement relative to point (x_r, 0), we obtain

[y_{1} (x) - y_{1} (x_{r})] + [y_{2} (x) - y_{2} (x_{r})] + [h (x) - h (x_{r})] {\begin{matrix} = 0 & x \in Ω_{c} \\ > 0 & x \notin Ω_{c} \end{matrix} .

(5)

Introducing a variable W(x) = h(x_r) – h(x), we obtain

\frac{2}{π E^{*}} \int_{Ω_{c}} \ln | \frac{x_{r} - τ}{x - τ} | P (τ) d τ = W (τ),

(6)

where $\frac{1}{E^{*}} = ((1 - υ_{1}^{2} / E_{1}) + (1 - υ_{2}^{2} / E_{2}))$ . E_i and μ_i (i = 1, 2) are Young’s modulus and Poisson’s ratio of elastic body i, respectively.

The contact pressure P needs to satisfy the following equation:

\int_{Ω_{c}} P (τ) d τ = Q .

(7)

Then, the contact stiffness is written as follows [45]:

k_{c} = \frac{Q}{δ} \cdot L,

(8)

where L is the contact length.

2.2. Numerical solution

As depicted in Figure 4, divided the contact zone into M interval, and there are M + 1 nodes, where every node is τ_j (j = 1, 2,…, M + 1). Based on linear interpolation method, the contact pressure P(τ) at interval [τ_j, τ_j+1] is written as follows:

P (τ) = \frac{τ - τ_{j}}{Δ τ} (P_{j + 1} - P_{j}) + P_{j} .

(9)

Figure 4.

Interpolation of contact pressure.

Substituting equation (9) into equations (6) and (7), the following equations are obtained:

\frac{2}{π E^{*}} \sum_{j = 1}^{M} \int_{τ_{j}}^{τ_{j + 1}} \ln | \frac{x_{r} - τ}{x - τ} | P (τ) d τ = W (x),

(10)

\frac{Δ τ}{2} \sum_{j = 1}^{M} (P_{j} + P_{j + 1}) = Q .

(11)

Amount of elastic deformation of any node due to contact pressure P_j can be expressed as follows:

\sum_{j = 1}^{M + 1} F_{ij} P_{j} = W_{i},

(12)

where F_ij is the influence coefficient. To obtain the influence coefficient, we need to integrate the following equation:

S_{j} (x) = \int_{τ_{j}}^{τ_{j + 1}} \ln | x - τ | P (τ) d τ .

(13)

Substituting equations (9) and (13) into equation (10), we obtain

W_{i} = \frac{2}{π E^{*}} \sum_{j = 1}^{M} [S_{j} (x_{r}) - S_{j} (x_{i})] = \frac{2}{π E^{*}} \sum_{j = 1}^{M + 1} [T_{j} (x_{r}) - T_{j} (x_{i})] P_{j} .

(14)

According to equations (12) and (14), the influence coefficients F_ij are written as follows:

F_{ij} = \frac{2}{π E^{*}} (T_{j} (x_{r}) - T_{j} (x_{i})) .

(15)

When the influence coefficients F_ij are determined, equations (11) and (12) show that the contact pressure can be addressed. The solution flow chart is presented in Figure 5, where a sufficiently large contact zone containing the actual contact zone is given and divided into M equal parts (M + 1 nodes). By calculating the influence coefficients F_ij and solving equations (11) and (12), we can obtain the contact pressure P_j (j = 1, 2,3,…, M + 1). It should be noted that the nodes whose contact pressure is less than 0 are removed from the contact zone and the contact pressure at these nodes is set to 0. Then, the process is repeated by updating the elastic approach δ until each P_j is no less than 0 and the relative difference between the total load and the given load is less than ε₀ (an extremely small value and can be set to 10^–6). Then, δ is obtained. The contact stiffness is determined by equation (8).

Figure 5.

Solution flow chart of elastic bodies contact.

3. TVMS modeling of spur gears

In the gear system, the deformations of axial compression, bending, fillet foundation, shear, and contact coexist. For a gear tooth under applied load Q, the deformation potential energy U_i (i = a, b, f, s) can be described as follows:

U_{i} = \frac{Q^{2}}{2 k_{i}},

(16)

where k refers to the stiffness. Subscript i (i = a, b, f, s) denotes the axial compression, bending, fillet foundation, and shear.

The contact deformation potential energy U_c between two elastic bodies is expressed as follows:

U_{c} = \frac{Q^{2}}{2 k_{c}},

(17)

where k_c is the contact stiffness.

For a contacting tooth pair, the total potential energy U_t can be obtained by Chen and Shao [46].

U_{t} = \sum U_{i} = U_{a 1} + U_{b 1} + U_{f 1} + U_{s 1} + U_{a 2} + U_{b 2} + U_{f 2} + U_{s 2} + U_{c},

(18)

where subscripts 1 and 2 represent pinion and gear, respectively.

Substituting equations (16) and (17) into equation (18), we obtain

U_{t} = \frac{Q^{2}}{2 k_{a 1}} + \frac{Q^{2}}{2 k_{b 1}} + \frac{Q^{2}}{2 k_{f 1}} + \frac{Q^{2}}{2 k_{s 1}} + \frac{Q^{2}}{2 k_{a 2}} + \frac{Q^{2}}{2 k_{b 2}} + \frac{Q^{2}}{2 k_{f 2}} + \frac{Q^{2}}{2 k_{s 2}} + \frac{Q^{2}}{2 k_{c}} .

(19)

Solving equations (16)–(19) simultaneously, we can express the meshing stiffness of contacting tooth pair i (i = 1, 2) as

k_{i} = \frac{1}{\frac{1}{k_{a 1}} + \frac{1}{k_{b 1}} + \frac{1}{k_{f 1}} + \frac{1}{k_{s 1}} + \frac{1}{k_{a 2}} + \frac{1}{k_{b 2}} + \frac{1}{k_{f 2}} + \frac{1}{k_{s 2}} + \frac{1}{k_{c}}} .

(20)

The TVMS of spur gears is given as follows:

k = {\begin{matrix} k_{1} & single - tooth contact \\ k_{1} + k_{2} & double - teeth contact \end{matrix} .

(21)

The cantilever beam model of the gear tooth is established in Figure 6. Point B represents the intersection of root circle and transition curve, and point C is the intersection of the transition and involute curves. The expressions of U_a, U_b, and U_s are written as follows [16]:

U_{a} = \int_{x_{B}}^{x_{C}} \frac{{(Q \sin β)}^{2}}{2 E A_{x 1}} d x_{1} + \int_{x_{C}}^{x_{β}} \frac{{(Q \sin β)}^{2}}{2 E A_{x 2}} d x_{2},

(22)

U_{b} = \int_{x_{B}}^{x_{C}} \frac{Q^{2} {[\cos β (x_{β} - x_{1}) - y_{β} \sin β]}^{2}}{2 E I_{x 1}} d x_{1} + \int_{x_{C}}^{x_{β}} \frac{Q^{2} {[\cos β (x_{β} - x_{2}) - y_{β} \sin β]}^{2}}{2 E I_{x 2}} d x_{2},

(23)

U_{s} = \int_{x_{B}}^{x_{C}} \frac{1.2 {(Q \cos β)}^{2}}{2 G A_{x 1}} d x_{1} + \int_{x_{C}}^{x_{β}} \frac{1.2 {(Q \cos β)}^{2}}{2 G A_{x 2}} d x_{2},

(24)

where x_i (i = B, C) is the X-coordinate of point i. y_β represents the Y-coordinate of arbitrary points on the involute curve. x₁ and x₂ are the X-coordinates of arbitrary points on the transition and involute curves, respectively. I_i, A_i (i = x₁, x₂) are the moment of inertia and cross-sectional area at coordinate x = i, respectively. G is the shear modulus.

Figure 6.

Cantilever beam model of gear tooth.

The transition curve coordinates (x₁, y₁) and involute tooth profile coordinates (x₂, y₂) are determined by Zhou and Wang [47]:

{\begin{matrix} x_{1} = (r - x_{c} - R_{o} \cos γ) \cos φ_{1} + (x_{c} \tan γ + R_{o} \sin γ) \sin φ_{1} \\ y_{1} = (r - x_{c} - R_{o} \cos γ) \sin φ_{1} - (x_{c} \tan γ + R_{o} \sin γ) \cos φ_{1} \end{matrix},

(25)

{\begin{matrix} x_{2} = [r - \frac{1}{2} (r φ_{2} - y_{0}) \sin (2 α)] \cos φ_{2} + (r φ_{2} - y_{0}) \cos^{2} α \sin φ_{2} \\ y_{2} = [r - \frac{1}{2} (r φ_{2} - y_{0}) \sin (2 α)] \sin φ_{2} - (r φ_{2} - y_{0}) \cos^{2} α \cos φ_{2} \end{matrix},

(26)

where the relevant parameters are listed in the work of Zhou and Wang [47].

According to equations (16) and (22)–(24), the following equations can be obtained by adopting the roll angle of the hob φ₁ and φ_2.

\frac{1}{k_{a}} = \int_{φ_{B}}^{φ_{C}} \frac{\sin^{2} β}{E A_{x 1}} \frac{d x_{1}}{d φ_{1}} d φ_{1} + \int_{φ_{C}}^{φ_{β}} \frac{\sin^{2} β}{E A_{x 2}} \frac{d x_{2}}{d φ_{2}} d φ_{2},

(27)

\frac{1}{k_{b}} = \int_{φ_{B}}^{φ_{C}} \frac{{[\cos β (x_{β} - x_{1}) - y_{β} \sin β]}^{2}}{E I_{x 1}} \frac{d x_{1}}{d φ_{1}} d φ_{1} + \int_{φ_{C}}^{φ_{β}} \frac{{[\cos β (x_{β} - x_{2}) - y_{β} \sin β]}^{2}}{E I_{x 2}} \frac{d x_{2}}{d φ_{2}} d φ_{2},

(28)

\frac{1}{k_{s}} = \int_{φ_{B}}^{φ_{C}} \frac{1.2 \cos^{2} β}{G A_{x 1}} \frac{d x_{1}}{d φ_{1}} d φ_{1} + \int_{φ_{C}}^{φ_{β}} \frac{1.2 \cos^{2} β}{G A_{x 2}} \frac{d x_{2}}{d φ_{2}} d φ_{2} .

(29)

where φ_B and φ_C are the hob roll angle of the points B and C, respectively.

The fillet foundation stiffness of the gear [12] is written as follows:

\frac{1}{k_{f}} = \frac{\cos^{2} β}{EB} {L^{*} {(\frac{u_{f}}{S_{f}})}^{2} + M^{*} (\frac{u_{f}}{S_{f}}) + P^{*} (1 + Q^{*} \tan^{2} β)},

(30)

where the parameters of u_f, S_f, L^*, M^*, and Q^* are presented by Sainsot et al. [12].

4. Results and discussion

4.1. Validation of proposed contact stiffness model

Hertzian contact theory is the most classical method used to solve the contact problem of elastic bodies. As shown in Figure 3, two elastic cylinders come into contact at point O under load Q. The half contact width a_H and maximum contact pressure P_max of the Hertzian contact model are given as follows:

a_{H} = \sqrt{\frac{4 Q (\frac{1 - μ_{1}^{2}}{E_{1}} + \frac{1 - μ_{2}^{2}}{E_{2}})}{π L (\frac{1}{R_{1}} + \frac{1}{R_{2}})}},

(31)

P_{max} = \sqrt{\frac{Q}{π L} \frac{(\frac{1}{R_{1}} + \frac{1}{R_{2}})}{(\frac{1 - μ_{1}^{2}}{E_{1}} + \frac{1 - μ_{2}^{2}}{E_{2}})}},

(32)

where R_i (i = 1, 2) is the radius of elastic body i, respectively.

The geometric and material parameters of two elastic cylinders are presented in Table 1. As shown in Figure 7(a), two elastic cylinders come into contact at the origin point O. According to the solving method for elastic body contact proposed in Section 2.2, taking 1.2 times the Hertzian contact width as the preset contact zone. Along the OX direction, discretized the preset contact zone into M (M = 69) points. The contact pressure distribution of the proposed model and Hertzian contact model are illustrated in Figure 7(b). Obviously the contact pressure distribution is almost the same in both models. In order to further study the effectiveness of the proposed model, the number of discrete points M into which the preset contact zone is divided is increased, and the results of maximum contact pressure and contact width are shown in Table 2. According to Hertzian contact theory, the maximum contact pressure and contact width are calculated as 837.82 MPa and 0.7598 mm, respectively. According to the proposed model, except for the maximum contact pressure is calculated as 837.88 MPa for M = 69, the relative differences with Hertzian contact is 0.0072%. The maximum contact pressures are all calculated as 837.82 MPa for M = 99, 199, and 299, and the relative differences with Hertzian contact are both almost zero. Furthermore, the contact width is calculated as 0.7665 mm for M = 69, 0.7737 mm for M = 99, 0.7698 mm for M = 199, and 0.7624 mm for M = 299. The differences relative to Hertzian contact are 0.882%, 1.829%, 1.316%, and 0.342%, respectively. The maximum contact pressure of the proposed model is almost the same as that of the Hertzian contact model, whereas a small deviation is observed between the results of the proposed contact width and the Hertzian contact width due to the limited discrete point. Nonetheless, the maximum relative difference is less than 0.4% for M = 299. Further extrapolation shows that the difference between the contact width of the proposed model and that of the Hertzian contact model will decrease continuously as M increases.

Table 1.

Parameters of contact of two elastic cylinders.

Parameters	Values	Parameters	Values
Cylindrical O₁ radius	R ₁ = 35 mm	Contact length	L = 30 mm
Cylindrical O₂ radius	R ₂ = 88 mm	Young’s modulus	E ₁ = E₂ = 201 GPa
Load	Q = 15,000 N	Poisson’s ratio	μ ₁ = μ₂ = 0.3

Figure 7.

Contact of two elastic cylinder: (a) initial contact position and (b) contact pressure.

Table 2.

Results of Hertzian contact and proposed method.

Item		Maximum contact pressure (MPa)		Contact width (mm)		Contact stiffness (N/mm)
Hertzian contact		837.82	Difference (%)	0.7598	Difference (%)	–
Proposed method	M = 69	837.88	0.0072	0.7665	0.882	1.0456 × 10⁶
	M = 99	837.82	0	0.7737	1.829	1.0660 × 10⁶
	M = 199	837.82	0	0.7698	1.316	1.0774 × 10⁶
	M = 299	837.82	0	0.7624	0.342	1.0888 × 10⁶

The FE method is further used to verify the proposed contact model. See Figure 8(a), a 2D FE model with plane strain elements (CPE4R) is presented in ABAQUS. The contact pressure distribution of the FE and proposed models is plotted in Figure 8(b). The number of mesh grids at the actual contact zone is approximately 50 in this FE model. The contact pressure distribution of the proposed model and FE contact model are almost the same. Owing to the limited number of mesh grids, the maximum contact pressure is calculated as 837.12 MPa in the FE model, and the relative difference between the proposed contact model and FE model is 0.0836%. Moreover, the contact width is calculated as 0.7879 mm in the FE model, and the relative difference is −3.2%.

Figure 8.

FE contact model of two elastic cylinders: (a) mesh grid model and (b) contact pressure between FE and proposed models.

The proposed contact model is verified through the preceding analysis, and the contact stiffness is listed in Table 2. Along with an increase in M, the contact stiffness increases slightly and is calculated as 1.0888 × 10⁶ N/mm for M = 299. When M increases from 69 to 299, the contact stiffness increases by only 4.1%. Moreover, as illustrated in Figure 9, the sum of the maximum deformation of the two elastic bodies in the FE model is 1.335 × 10^–2 mm, and the contact stiffness is calculated as 1.1236 × 10⁶ N/mm. The relative difference between the proposed and FE model is −3.1%, which further demonstrates the correctness of the proposed model.

Figure 9.

Results of contact displacement by FE model.

The contact stiffness of the proposed model and FE model across different applied load are illustrated in Figure 10. In both models, the contact stiffness increase with an increasing applied load. However, with the increase of load, the increasing rate of the contact stiffness decreases gradually. Further analysis of Figure 10 shows that even if the contact stiffness increases with the increase of the load, the increase value of the contact stiffness is not large. During the change of the load from 3000 to 240,000 N, the contact stiffness only increases from 9.4159 × 10⁵ N/mm to 1.5662 × 10⁶ N/mm. This means that the load is increased by a factor of 80, while the value of the contact stiffness is only increased by a factor of 1.6634. Moreover, the minimum relative difference of proposed contact stiffness and the FE contact stiffness model is −1.64% and the maximum relative difference is −3.81%. Therefore, the validity of the proposed contact stiffness model is further verified.

Figure 10.

Contact stiffness across different applied load.

4.2. Comparison of different contact stiffness models

In the previous studies, contact stiffness is defined by SHCS [36] or NHCS [40]. The expression of SHCS is

k_{h} = \frac{π E^{*} L}{2} .

(33)

The NHCS is written as follows:

k_{nh} = \frac{{E_{eq}}^{0.9} L^{0.8} Q^{0.1}}{1.275},

(34)

where 1/E_eq = (1/E₁ + 1/E₂)/2.

In accordance with equations (33)–(34) and the parameters in Table 1, the SHCS and NHCS are obtained as 5.2043 × 10⁶ and 1.8479 × 10⁶ N/mm, respectively. This condition indicates that the contact stiffness derived from the proposed model is approximately 0.21 times that of SHCS and 0.59 times that of NHCS. Proposed contact stiffness, SHCS, NHCS, and FE contact stiffness across different applied load are plotted in Figure 11. Since the influence of load is not considered in SHCS, SHCS does not change with the change of load. Proposed contact stiffness, NHCS, and FE contact stiffness increase with the increase of load, and they almost have the same variation trend. While NHCS is greater than proposed contact stiffness at the same load. As the load increased from 3000 to 240,000 N, SHCS increases by a minimum of 232.29% and a maximum of 452.72% over proposed contact stiffness, and NHCS increases by a minimum of 55.68% and a maximum of 71.41% over proposed contact stiffness. Therefore, in engineering application, the use of SHCS and NHCS will make the value of contact stiffness larger. The minimum relative difference between SHCS and FE model is 224.84% and the maximum relative difference between SHCS and FE model is 443.68%. The minimum relative difference of NHCS and the FE model is 52.19% and the maximum relative difference is 64.34%. According to the analysis in Section 4.1, the minimum relative difference of the proposed contact stiffness and the FE model is −1.64%, and the maximum relative difference is −3.81%. Therefore, it can be concluded that the stiffness model proposed in this work is more accurate.

Figure 11.

Three typical contact stiffness across different applied load.

4.3. TVMS of spur gear drive

Since the contact stiffness is related to the applied load, an iterative method is needed to calculate the TVMS of spur gear pair. The calculation flow chart of TVMS of spur gear pair is shown in Figure 12. A meshing process of a gear pair is divided into N_d nodes. The axial compression stiffness (k_a₁, k_a₂), bending stiffness (k_b₁, k_b₂), fillet foundation (k_f₁, k_f₂), and shear stiffness (k_s₁, k_s₂) can be determined by the gear design parameters because they are independent of the load. Then, the initial TVMS is obtained by adopting SHCS, and the initial load sharing factor (l_sf₁) of one meshing gear pair is written as follows:

l_{sf 1} = \frac{k_{1}}{k} = {\begin{matrix} 1 & single - tooth contact \\ \frac{k_{1}}{k_{1} + k_{2}} & double - teeth contact \end{matrix}

(35)

Figure 12.

Solution flow chart of TVMS of spur gear drive.

One meshing cycle is averagely divided into N_d parts, the initial load sharing factor (l_q, _sf1) of the gear pair at q position can be obtained, and the contact load will be determined. According to the proposed contact stiffness model and the energy method stiffness model, TVMS of the gear pair at q position can be obtained, and the load sharing factor l’_{q, sf1} can be obtained according to equation (35). TVMS of the gear pair at q position can be output if l_q, _sf1 and l’_{q, sf1} satisfy the following equation:

| (l_{q, sf 1} - l_{q, sf 1}^{'}) / l_{q, sf 1} | < ε_{1}

(36)

where ε₁ is a small value that can be set to 10^–6.

If l_q, _sf1 and l’_{q, sf1} do not satisfy equation (36), set l_q, _sf1 = l’_{q, sf1}, and repeat the above process until l_q, _sf1 and l’_{q, sf1} satisfy equation (36). Then, the TVMS at the next meshing position is calculated according to the above process until q = N_d.

The design parameters of a spur gear drive are presented in Table 3, and the mesh stiffness of single-tooth pair and the TVMS values of the gear pair with SHCS, NHCS, and proposed contact stiffness are depicted in Figure 13. Since TVMS of gear pair obtained by NHCS is also related to load, and the load of gear pair in the meshing process is time-varying. Therefore, the determination of TVMS of gear pair obtained by NHCS should be carried out in accordance with the calculation flow chart of TVMS obtained by the proposed contact stiffness in Figure 12, just need to change the k_c by k_nh.

Table 3.

Design parameters of a spur gear drive.

Parameters	Values	Parameters	Values
Tooth number	N ₁ = 31, N₂ = 48	Hub bore radius	r _o1 = 25 mm, r_o2 = 25 mm
Module	m = 4 mm	Input torque	T_in = 300 N m
Pressure angle	α = 20°	Young’s modulus	E ₁ = 206 GPa, E₂ = 206 GPa
Tooth width	B_w = 30 mm	Poisson’s ratio	μ ₁ = 0.3, μ₂ = 0.3

Figure 13.

Mesh stiffness: (a) mesh stiffness of single tooth pair and (b) TVMS of gear drive.

As shown in Figure 13(a), the mesh stiffness of the single-tooth pair obtained by SHCS is the largest while that of the single tooth pair obtained by proposed contact stiffness is the smallest. Proposed contact stiffness and NHCS are affected by the load, and the load of spur gears mutates at the inflection point of the single-tooth and double-tooth contact. Therefore, the mesh stiffness obtained by the proposed contact stiffness and NHCS show a small mutation at these two meshing positions. The variation trends of the TVMS obtained by the three methods are the same (Figure 13(b)) because only the calculation method of the contact stiffness is changed. Based on the parameters given in Table 3, the TVMS calculated by the proposed contact stiffness model is the smallest and that by SHCS is the largest. At double-tooth contact, the difference between the maximum TVMS calculated by SHCS and proposed contact stiffness is 9.8382 × 10⁷ N/m, the relative difference is 18.69%, the difference between the maximum TVMS calculated by NHCS and proposed contact stiffness is 2.2318 × 10⁷ N/m, and the relative difference is 4.24%. At single-tooth contact, the difference between the maximum TVMS calculated by SHCS and proposed contact stiffness is 5.8072 × 10⁷ N/m, and the relative difference is 19.94%; the difference between the maximum TVMS calculated by NHCS and by the proposed contact stiffness is 1.5108 × 10⁷ N/m, and the relative difference is 5.19%. These results mean that the differences of TVMS by SHCS or NHCS are relatively large.

The TVMS of the gear drive is generally influenced by the contact stiffness, which is significantly determined by tooth width and input torque. To aid gear designers in maximizing load capacity and optimizing dynamic performance, we discuss the effects of these parameters on the TVMS of the spur gear drive in the following section.

4.4. Effects of tooth width

The mesh stiffness of the single tooth pair and TVMS of the gear drive across different tooth widths are depicted in Figure 14. Along with an increase in tooth width, the mesh stiffness of the single-tooth pair increases. As illustrated in Figure 14(a), the maximum mesh stiffness of the single tooth pair is calculated as 1.9828 × 10⁸ N/m for B_w = 20 mm, 2.9120 × 10⁸ N/m for B_w = 30 mm, and 4.7524 × 10⁸ N/m for B_w = 50 mm. Even though the influence of the tooth width on axial compression stiffness, bending stiffness, fillet foundation stiffness, and shear mesh stiffness is linear, the mesh stiffness of the single-tooth pair non-linearly increases with the tooth width because the contact stiffness is related to the load. The TVMS of gear drive becomes larger with increasing tooth width is illustrated in Figure 14(b). The maximum mesh stiffness of the gear drive at the double-tooth contact is calculated as 3.5769 × 10⁸ N/m for B_w = 20 mm, 5.2632 × 10⁸ N/m for B_w = 30 mm, and 8.6054 × 10⁸ N/m for B_w = 50 mm. This means that although TVMS increases with the increase of load, the rate of increase is less and less.

Figure 14.

Mesh stiffness across different tooth widths: (a) mesh stiffness of single tooth pair and (b) TVMS of gear drive.

4.5. Effects of input torque

Figure 15 illustrates the mesh stiffness of the single tooth pair and the TVMS of the gear drive across different input torques. The mesh stiffness of the single tooth pair increases with an increasing input torque. As shown in Figure 15(a), the maximum mesh stiffness of the single tooth pair is calculated as 2.8570 × 10⁸ N/m for T_in = 200 N m, 2.9120 × 10⁸ N/m for T_in = 300 N m, and 2.9979 × 10⁸ N/m for T_in = 500 N m. Since the axial compression stiffness, bending stiffness,0 fillet foundation stiffness, and shear mesh stiffness are independent of the load and the contact stiffness depends on the load, TVMS vary with the load. The TVMS of gear drive becomes larger with increasing input torque is illustrated in Figure 15(b). The maximum mesh stiffness of the gear drive at the double-tooth contact is calculated as 5.1772 × 10⁸ N/m for T_in = 200 N m, 5.2632 × 10⁸ N/m for T_in = 300 N m, and 5.3994 × 10⁸ N/m for T_in = 500 N m. These results indicate that the TVMS of the gear drive increases by 1.66% as the input torque increases from 200 to 300 N m, and 2.59% as the input torque increases from 300 to 500 N m. It indicates that the increase rate of TVMS with the increase of torque becomes smaller and smaller.

Figure 15.

Mesh stiffness across different input torques: (a) mesh stiffness of single tooth pair and (b) TVMS of gear drive.

5. Conclusion

In this study, a numerical contact stiffness model for elastic bodies is developed and applied to the TVMS modeling of the spur gear drive. According to the theory of elasticity, the contact stiffness model of elastic bodies is obtained by a numerical model that employs a piecewise linear interpolation method. The TVMS of the spur gears is determined in accordance with the stiffness models of axial compression, bending, fillet foundation, and shear which are derived by the potential energy method and proposed contact stiffness model. The proposed contact model is verified by Hertzian contact theory and FE method, and the differences among proposed contact stiffness, SHCS, and NHCS are analyzed. Then, the TVMS values of a spur gear drive with three contact stiffness models are further investigated. The effects of tooth width and input torque on the TVMS of the gear drive are discussed. The main conclusions of this study are as follows:

The proposed contact stiffness model of elastic bodies is verified by the Hertzian contact and FE models. Even though the contact stiffness increases with the load, the increasing rate decreases gradually.

Compared with the SHCS model and NHCS model, the proposed contact stiffness is the minimum, which means that the TVMS of the gear pair by the Hertzian simplified contact stiffness or NHCS are larger.

The TVMS of the gear pair increases with the increased tooth width and input torque, but the increase rate decreases as the tooth width or input torque increases. These findings indicate that reasonable matching of design parameters is beneficial for increasing gear load capacity and optimizing dynamic performance.

The proposed contact stiffness model for elastic bodies can be well extended to cylindrical roller bearing, cam, and friction wheel, which have line contact.

Footnotes

Appendix 1 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: This study was funded by the National Natural Science Foundation of China (grant nos: 52075153 and 52005051), the Key Research and Development Program of Hunan Province (grant no: 2020WK2032), and the Natural Science Foundation of Hunan Province (grant no: 2021JJ40580).

ORCID iD

Changjiang Zhou

References

Ambarisha

Parker

RG.

Nonlinear dynamics of planetary gears using analytical and finite element models. J Sound Vib 2007; 302: 577–595.

Zhou

Wang

, et al. Prediction and validation of dynamic characteristics of a valve train system with flexible components and gyroscopic effect. Mech Mach Theory 2021; 157: 104222.

Feng

, et al. Experimental and theoretical investigation on rotodynamic characterization of hybrid porous tilting pad bearings. Mech Syst Signal Pr 2022; 164: 108245.

Zhou

Zhao

Tang

, et al. A comparative study on asperity peak modeling methods. Chin J Mech Eng 2021; 34: 61.

Chen

Zhai

Shao

, et al. Analytical model for mesh stiffness calculation of spur gear pair with non-uniformly distributed tooth root crack. Eng Fail Anal 2016; 66: 502–514.

Wang

Tang

Zhou

, et al. Wear life prediction method of crowned double helical gear drive in point contact mixed elastohydrodynamic lubrication. Wear 2021; 484–485: 204041.

Chen

Zhou

Zhai

, et al. Improved analytical calculation model of spur gear mesh excitations with tooth profile deviations. Mech Mach Theory 2020; 149: 103838.

Chen

Ning

Wang

, et al. An improved dynamic model of spur gear transmission considering coupling effect between gear neighboring teeth. Nonlinear Dynam 2021; 106: 339–357.

Pandya

Parey

Experimental investigation of spur gear tooth mesh stiffness in the presence of crack using photoelasticity technique. Eng Fail Anal 2013; 34: 488–500.

10.

Raghuwanshi

Parey

Experimental measurement of gear mesh stiffness of cracked spur gear by strain gauge technique. Measurement 2016; 86: 266–275.

11.

Raghuwanshi

Parey

Experimental measurement of spur gear mesh stiffness using digital image correlation technique. Measurement 2017; 111: 93–104.

12.

Sainsot

Velex

Duverger

Contribution of gear body to tooth deflections-a new bidimensional analytical formula. Trans Am Soc Mech Eng J Mech Des 2004; 126: 748–752.

13.

Wang

Howard

Finite element analysis of high contact ratio spur gears in mesh. J Tribol 2005; 127: 469–483.

14.

Jia

Howard

Comparison of localised spalling and crack damage from dynamic modelling of spur gear vibrations. Mech Syst Signal Pr 2006; 20: 332–349.

15.

Pang

Feng

, et al. Improved time-varying mesh stiffness model of cracked spur gears. Eng Fail Anal 2015; 55: 271–287.

16.

Song

Pang

, et al. Time-varying mesh stiffness calculation of cracked spur gears. Eng Fail Anal 2014; 44: 179–194.

17.

Zeng

Feng

, et al. An improved analytical method for mesh stiffness calculation of spur gears with tip relief. Mech Mach Theory 2016; 98: 64–80.

18.

Chen

Zhang

Zhai

, et al. Improved analytical methods for calculation of gear tooth fillet-foundation stiffness with tooth root crack. Eng Fail Anal 2017; 82: 72–81.

19.

Feng

, et al. Meshing characteristics of spur gear pair under different crack types. Eng Fail Anal 2017; 80: 123–140.

20.

Sun

Huangfu

, et al. A revised time-varying mesh stiffness model of spur gear pairs with tooth modifications. Mech Mach Theory 2018; 129: 261–278.

21.

Wang

Chen

Zhao

, et al. An analytical-finite-element method for calculating mesh stiffness of spur gear pairs with complicated foundation and crack. Eng Fail Anal 2018; 94: 339–353.

22.

Xie

Hua

Lan

, et al. Improved analytical models for mesh stiffness and load sharing ratio of spur gears considering structure coupling effect. Mech Syst Signal Pr 2018; 111: 331–347.

23.

Yang

Lin

Hertzian damping, tooth friction and bending elasticity in gear impact dynamics. J Mech Transm Autom Des 1987; 109: 189–196.

24.

Zuo

Parey

Simulation of spur gear dynamics and estimation of fault growth. J Sound Vib 2008; 317: 608–624.

25.

Chaari

Fakhfakh

Haddar

Analytical modelling of spur gear tooth crack and influence on gearmesh stiffness. Eur J Mech A: Solid 2009; 28: 461–468.

26.

Chen

Shao

Dynamic simulation of spur gear with tooth root crack propagating along tooth width and crack depth. Eng Fail Anal 2011; 18: 2149–2164.

27.

Mohammed

Rantatalo

Aidanpää

J-O.

Improving mesh stiffness calculation of cracked gears for the purpose of vibration-based fault analysis. Eng Fail Anal 2013; 34: 235–251.

28.

Wan

Cao

, et al. An improved time-varying mesh stiffness algorithm and dynamic modeling of gear-rotor system with tooth root crack. Eng Fail Anal 2014; 42: 157–177.

29.

Saxena

Parey

Chouksey

Time varying mesh stiffness calculation of spur gear pair considering sliding friction and spalling defects. Eng Fail Anal 2016; 70: 200–211.

30.

Wan

Cao

, et al. Mesh stiffness calculation using an accumulated integral potential energy method and dynamic analysis of helical gears. Mech Mach Theory 2015; 92: 447–463.

31.

Feng

, et al. An improved analytical method for calculating time-varying mesh stiffness of helical gears. Meccanica 2018; 53: 1131–1145.

32.

Huangfu

Chen

, et al. Deformation and meshing stiffness analysis of cracked helical gear pairs. Eng Fail Anal 2019; 95: 30–46.

33.

Mechefske

A new model for the single mesh stiffness calculation of helical gears using the slicing principle. Iran J Sci Technol Trans Mech Eng 2019; 43: 503–515.

34.

Zhou

Wang

, et al. Nonlinear tribo-dynamic model and experimental verification of a spur gear drive under loss-of-lubrication condition. Mech Syst Signal Pr 2021; 153: 107509.

35.

Tian

Dynamic simulation for system response of gearbox including localized gear faults. Master’s Thesis, University of Alberta, Edmonton, AB, Canada, 2004.

36.

Yang

Sun

A rotary model for spur gear dynamics. J Mech Transm Autom Des 1985; 107: 529–535.

37.

Mohammed

Rantatalo

Aidanpää

J-O.

Dynamic modelling of a one-stage spur gear system and vibration-based tooth crack detection analysis. Mech Syst Signal Pr 2015; 54–55: 293–305.

38.

Pang

Feng

, et al. Evaluation of optimum profile modification curves of profile shifted spur gears based on vibration responses. Mech Syst Signal Pr 2016; 70–71: 1131–1149.

39.

Wang

Zhou

Wang

, et al. A novel contact model for rough surfaces using piecewise linear interpolation and its application in gear wear. Wear 2021; 476: 203685.

40.

Cornell

Compliance and stress sensitivity of spur gear teeth. J Mech Des 1981; 103: 447–459.

41.

Xie

Hua

Han

, et al. Analytical formulas for gear body-induced tooth deflections of spur gears considering structure coupling effect. Int J Mech Sci 2018; 148: 174–190.

42.

Sainsot

Velex

On contact deflection and stiffness in spur and helical gears. Mech Mach Theory 2020; 154: 104049.

43.

Dai

Long

Chen

, et al. An improved analytical model for gear mesh stiffness calculation. Mech Mach Theory 2021; 159: 104262.

44.

Webster

Sayles

RS.

A numerical model for the elastic frictionless contact of real rough surfaces. J Tribol 1986; 108: 314–320.

45.

Cao

Chen

Jiang

Nonlinear dynamics of a spur gear pair with force-dependent mesh stiffness. Nonlinear Dynam 2020; 99: 1227–1241.

46.

Chen

Shao

Mesh stiffness calculation of a spur gear pair with tooth profile modification and tooth root crack. Mech Mach Theory 2013; 62: 63–74.

47.

Zhou

Wang

An adhesive wear prediction method for double helical gears based on enhanced coordinate transformation and generalized sliding distance model. Mech Mach Theory 2018; 128: 58–83.