Abstract
As one of the most critical components in transmission systems, spur gears are widely applied in aerospace, marine, precision machine tools, and transportation equipment. Under cyclic meshing loads, fatigue cracks often initiate at the tooth root fillet where stress concentration is severe. Once a crack appears, it propagates along a certain path until tooth fracture, causing unexpected failures and economic losses. This paper systematically investigates the crack propagation mechanism, the influence of gear blank structural parameters on crack paths, the time-varying meshing stiffness (TVMS) of cracked spur gears, the dynamic response of a cracked gear pair, and the extraction of fault features from vibration signals. The extended finite element method (XFEM) is improved by introducing enriched shape functions and a hole-enrichment function to handle the discontinuity and singularity at the crack tip. The improved XFEM model is validated by comparing stress intensity factors (SIF) with literature results. The effects of crack nucleation position, web outer diameter, and the relative position between the web hole and the crack nucleation point on the crack propagation path are studied. The J-integral method is used to calculate the SIF. A modified energy method based on a cantilever beam model is adopted to compute the TVMS of healthy and cracked spur gears under different crack depths and crack angles. A 6-degree-of-freedom (6-DOF) lumped parameter dynamic model is established, and the Runge-Kutta method is employed to solve the vibration differential equations. Gear crack experiments with artificial cracks of different depths are carried out using wire electrical discharge machining. The vibration signals are analyzed in time and frequency domains, and the envelope spectrum is obtained by Hilbert transform. Furthermore, a fault feature extraction method combining adaptive variational mode decomposition (VMD) and K-singular value decomposition (KSVD) is proposed to extract the weak fault features from noisy signals. The results show that the improved XFEM can effectively simulate spur gear root crack propagation; the crack path tends to propagate toward the gear rim when the nucleation point is closer to the root circle, when the web outer diameter increases, or when the web hole is closer to the crack. The SIF increases with crack growth and is dominated by mode I. The TVMS decreases with increasing crack depth and crack angle, and the crack depth has a more significant effect than the crack angle. The dynamic response reveals that the impact amplitude increases with crack depth, and the envelope spectrum shows a clear growth of the shaft rotational frequency amplitude and sidebands. The proposed VMD-KSVD method can automatically determine the number of intrinsic mode functions and the penalty factors, and the fault characteristic frequencies and amplitudes of each IMF increase with crack severity. These findings provide valuable insights for spur gear crack propagation prediction, dynamic modeling, and fault diagnosis.
Keywords: spur gears; crack propagation; extended finite element method; time-varying meshing stiffness; vibration characteristics; fault feature extraction; variational mode decomposition; K-SVD
1. Introduction
Spur gears are fundamental power transmission elements in rotating machinery. Their reliability directly affects the safety of the entire system. Among gear failure modes, tooth root fatigue crack and subsequent tooth fracture are the most common and dangerous. Cracks usually initiate at the tooth root fillet due to high stress concentration caused by cyclic bending and shear stresses. Under repeated loading, the crack propagates through the gear tooth or rim, eventually leading to tooth breakage or rim fracture. Understanding the crack propagation path and the dynamic behavior of cracked spur gears is essential for gear design, condition monitoring, and life prediction.
The crack propagation path is affected by many factors, including the initial crack length, initial crack angle, gear blank geometry, web structure, web holes, load conditions, and material properties. Traditional finite element method (FEM) requires remeshing at each crack growth step, which is computationally expensive and less efficient. The extended finite element method (XFEM) overcomes this drawback by allowing the crack to propagate independently of the mesh. However, the standard XFEM may suffer from issues near the crack tip due to the presence of blended (mixing) elements. Therefore, a modified XFEM is needed to improve the enrichment functions, especially for spur gears with web holes.
For the vibration analysis of cracked spur gears, the time-varying meshing stiffness is the most important internal excitation. Numerous studies have proposed analytical methods to calculate TVMS, including the potential energy method, the Weber energy method, and the Ishikawa method. The potential energy method models a gear tooth as a cantilever beam and considers Hertzian, bending, shear, axial compressive, and foundation deformations. When a crack exists, the bending and shear stiffness change significantly. Thus, accurate calculation of TVMS for different crack depths and crack angles is necessary for dynamic analysis.
The dynamic response of a gear system with a cracked tooth contains characteristic sidebands around the meshing frequency and its harmonics. However, measured vibration signals are often contaminated by noise and other interfering sources. Therefore, advanced signal processing techniques are required to extract the weak fault features. Variational mode decomposition (VMD) is a robust non-recursive decomposition method that can separate signal components with narrow bandwidths. But the number of modes K and the penalty factor α must be chosen properly. The K-singular value decomposition (KSVD) algorithm is a dictionary learning method that can denoise the signal and preserve its transient characteristics. Combining VMD and KSVD with an adaptive scale-space approach can automatically determine the parameters and improve the fault feature extraction efficiency.
This paper aims to provide a comprehensive study on spur gear root crack propagation and the vibration characteristics of cracked spur gears. The main contributions are:
- Deriving the stress and displacement fields at the crack tip using the Westergaard stress function and the complex variable method.
- Improving the XFEM with two enriched shape functions and a hole-enrichment function to accurately simulate the crack propagation path in spur gears with webs and holes.
- Simulating the influence of crack nucleation position, web outer diameter, and web hole position on the crack path and SIF.
- Calculating the TVMS of healthy and cracked spur gears using an improved energy method.
- Establishing a 6-DOF dynamic model and solving it using the Runge-Kutta method to obtain the vibration response.
- Validating the dynamic characteristics through experiments and extracting fault features using an adaptive VMD-KSVD method.
2. Crack Propagation Fundamentals
2.1 Crack Modes
Depending on the loading direction relative to the crack surface, cracks are classified into three modes: opening mode (mode I), sliding mode (mode II), and tearing mode (mode III). In spur gear tooth root cracks, the crack is subjected to both normal stress and shear stress during meshing. Therefore, the crack is typically a mixed-mode crack composed of mode I and mode II. The stress field near the crack tip can be derived from the Westergaard stress function using the theory of complex functions. The stress components in the local crack-tip coordinate system are expressed as:
$$
\begin{cases}
\sigma_x = \dfrac{K_I}{\sqrt{2\pi r}} \cos\dfrac{\theta}{2}\left(1 – \sin\dfrac{\theta}{2}\sin\dfrac{3\theta}{2}\right) – \dfrac{K_{II}}{\sqrt{2\pi r}} \sin\dfrac{\theta}{2}\left(2 + \cos\dfrac{\theta}{2}\cos\dfrac{3\theta}{2}\right) \\[6pt]
\sigma_y = \dfrac{K_I}{\sqrt{2\pi r}} \cos\dfrac{\theta}{2}\left(1 + \sin\dfrac{\theta}{2}\sin\dfrac{3\theta}{2}\right) + \dfrac{K_{II}}{\sqrt{2\pi r}} \sin\dfrac{\theta}{2}\cos\dfrac{\theta}{2}\cos\dfrac{3\theta}{2} \\[6pt]
\tau_{xy} = \dfrac{K_I}{\sqrt{2\pi r}} \cos\dfrac{\theta}{2}\sin\dfrac{\theta}{2}\cos\dfrac{3\theta}{2} + \dfrac{K_{II}}{\sqrt{2\pi r}} \cos\dfrac{\theta}{2}\left(1 – \sin\dfrac{\theta}{2}\sin\dfrac{3\theta}{2}\right)
\end{cases}
$$
where \(K_I\) and \(K_{II}\) are the mode I and mode II stress intensity factors (SIF), \(r\) is the distance from the crack tip, and \(\theta\) is the angle measured from the crack plane. The displacement components are given by:
$$
\begin{cases}
u_x = \dfrac{(1+\nu)}{E}\sqrt{\dfrac{r}{2\pi}}\left[K_I\cos\dfrac{\theta}{2}\left(\dfrac{3-\nu}{1+\nu} – \cos\theta\right) + K_{II}\sin\dfrac{\theta}{2}\left(\dfrac{5+\nu}{1+\nu} + \cos\theta\right)\right]\\[8pt]
u_y = \dfrac{(1+\nu)}{E}\sqrt{\dfrac{r}{2\pi}}\left[K_I\sin\dfrac{\theta}{2}\left(\dfrac{3-\nu}{1+\nu} – \cos\theta\right) – K_{II}\cos\dfrac{\theta}{2}\left(\dfrac{1-3\nu}{1+\nu} + \cos\theta\right)\right]
\end{cases}
$$
Here, \(E\) is the elastic modulus and \(\nu\) is Poisson’s ratio.
2.2 Stress Intensity Factor Calculation
The SIF characterizes the singular stress field at the crack tip. In this study, the J-integral method is adopted because of its path independence, which avoids the need for extremely fine meshes at the crack tip. The J-integral is defined as:
$$
J = \int_L \left( w \, dy – \left( T_x \dfrac{\partial u_x}{\partial x} + T_y \dfrac{\partial u_y}{\partial x} \right) ds \right)
$$
where \(L\) is a closed contour around the crack tip, \(w\) is the strain energy density, \(T_x\) and \(T_y\) are traction components on the contour, and \(u_x\), \(u_y\) are displacement components. For mixed-mode cracks, the J-integral relates to the SIFs as:
$$
J = \frac{K_I^2 + K_{II}^2}{E’}
$$
where \(E’ = E\) for plane stress and \(E’ = E/(1-\nu^2)\) for plane strain.
2.3 Crack Propagation Criterion
The maximum circumferential tensile stress criterion is widely used for predicting the direction of crack propagation. The circumferential stress at the crack tip is:
$$
\sigma_\theta = \frac{1}{\sqrt{2\pi r}}\cos\frac{\theta}{2}\left[\frac{K_I}{2}(1+\cos\theta) – \frac{3K_{II}}{2}\sin\theta\right]
$$
Define the circumferential stress intensity factor as:
$$
K_\theta = \lim_{r\to0}\sqrt{2\pi r}\,\sigma_\theta = \frac{1}{2}\cos\frac{\theta}{2}\left[K_I(1+\cos\theta) – 3K_{II}\sin\theta\right]
$$
The crack propagation direction \(\theta_0\) is found by maximizing \(K_\theta\):
$$
\frac{\partial K_\theta}{\partial \theta}\Big|_{\theta=\theta_0} = 0, \quad \frac{\partial^2 K_\theta}{\partial \theta^2}\Big|_{\theta=\theta_0} < 0
$$
Solving the above equation yields the crack deflection angle:
$$
\theta_0 = 2\arctan\left[\frac{1}{4}\left(\frac{K_I}{K_{II}} \pm \sqrt{\left(\frac{K_I}{K_{II}}\right)^2 + 8}\right)\right]
$$
When the maximum circumferential stress intensity factor reaches the critical value \(K_{\theta C}\), crack instability occurs.
3. Improved Extended Finite Element Method for Spur Gear Crack Propagation
3.1 Basics of XFEM
The extended finite element method is based on the partition of unity concept. The displacement field is enriched with additional shape functions to represent the discontinuity of the crack surface and the singularity of the crack tip. The standard XFEM displacement approximation is:
$$
u(\mathbf{x}) = \sum_{i \in S} N_i(\mathbf{x})\,u_i + \sum_{i \in S_h} N_i(\mathbf{x})\,H(f(\mathbf{x}))\,a_i + \sum_{i \in S_c} N_i(\mathbf{x})\,\phi(\mathbf{x})\,b_i
$$
where \(N_i(\mathbf{x})\) are standard finite element shape functions, \(u_i\) are nodal displacement vectors, \(H\) is the Heaviside step function, \(f(\mathbf{x})\) is the level set function, \(\phi(\mathbf{x})\) are crack-tip enrichment functions, \(a_i\) and \(b_i\) are additional degrees of freedom, \(S_h\) is the set of nodes enriched by the Heaviside function, and \(S_c\) is the set of nodes enriched by the crack-tip functions. The crack-tip enrichment functions are:
$$
\phi(\mathbf{x}) = \left[\sqrt{r}\sin\frac{\theta}{2},\; \sqrt{r}\sin\frac{\theta}{2}\sin\theta,\; \sqrt{r}\cos\frac{\theta}{2},\; \sqrt{r}\cos\frac{\theta}{2}\sin\theta\right]
$$
3.2 Modified Enrichment Functions
In the standard XFEM, when only a part of an element is enriched, the partition of unity property is not satisfied in the blended element. This can lead to a loss of convergence. To overcome this, the enriched functions are shifted and a ramp function is introduced. The modified Heaviside enrichment is:
$$
P_1(\mathbf{x}) = \sum_{i} N_i(\mathbf{x})\left[H(f(\mathbf{x})) – H(f(\mathbf{x}_i))\right]
$$
and the modified crack-tip enrichment is:
$$
Q_1(\mathbf{x}) = \sum_{i} N_i(\mathbf{x})\,\phi(\mathbf{x})\,R(\mathbf{x})
$$
where \(R(\mathbf{x})\) is a gradual (ramp) function that is zero at the boundary of the enriched region and unity inside the fully enriched element. Moreover, for spur gears with web holes, an additional enrichment term is introduced to represent the hole discontinuity:
$$
F(\mathbf{x}) = \begin{cases} 1, & \mathbf{x} \in \Omega^+ \\ 0, & \mathbf{x} \in \Omega^- \end{cases}
$$
The final improved displacement approximation becomes:
$$
u^h(\mathbf{x}) = \sum_{i \in S} N_i(\mathbf{x})\,u_i + \sum_{i \in S_h} N_i(\mathbf{x})\left[H(f(\mathbf{x})) – H(f(\mathbf{x}_i))\right]a_i + \sum_{i \in S_c} N_i(\mathbf{x})\,\phi(\mathbf{x})\,R(\mathbf{x})\,b_i + \sum_{i \in S_h} N_i(\mathbf{x})\,F(\mathbf{x})\,c_i
$$
where \(S_h\) is the set of nodes enriched by the hole function and \(c_i\) are the corresponding additional degrees of freedom.
3.3 Finite Element Model of Spur Gear
To simulate the crack propagation path in spur gears, a two-tooth segment of the driving gear is modeled using ANSYS Workbench. A fine mesh is generated in the tooth root and the web hole areas. The total model consists of 6,576 elements and 8,758 nodes. One side of the model is fully constrained in both \(x\) and \(y\) directions, while the rotational degree of freedom is also constrained. A torque \(T = 141.9\,\mathrm{N\cdot m}\) and an angular velocity \(\Omega = 1600\,\mathrm{rad/s}\) are applied. The meshing force is decomposed into tangential and radial components:
$$
F_t = \frac{2000\,T}{d}, \quad F_r = F_t \tan 20^\circ
$$
with \(d = (z + 2h_a^*) m\). The following table summarizes the gear pair parameters used in the simulation.
| Parameter | Value |
|---|---|
| Number of teeth \(z\) | 28 |
| Module \(m\) (mm) | 4 |
| Face width \(B\) (mm) | 20 |
| Web thickness \(B_w\) (mm) | 10 |
| Web inner radius \(r_{wi}\) (mm) | 36 |
| Web outer radius \(r_{wo}\) (mm) | 40–48 |
| Web hole diameter \(\Phi_s\) (mm) | 12 |
| Hole-crack angle \(\theta_P\) (deg) | -30 to 30 |
| Elastic modulus \(E\) (GPa) | 210 |
| Poisson’s ratio \(\nu\) | 0.3 |
| Density (g/cm³) | 7.83 |
3.4 Validation of the Improved XFEM
To verify the accuracy of the improved XFEM, the crack propagation path is compared with the finite element result reported in the literature. The paths are almost identical with a slight deviation in the later growth stage, confirming the effectiveness of the method. In addition, the stress intensity factors computed by the improved XFEM are compared with literature values for three different crack nucleation positions. The maximum relative error is 4.73%, which is within the acceptable engineering range. The comparison is shown in Table 1.
| Nucleation position | SIF | Literature (MPa·mm1/2) | Improved XFEM (MPa·mm1/2) | Error (%) |
|---|---|---|---|---|
| 1 | \(K_I\) | 848.91 | 868.60 | 2.32 |
| \(K_{II}\) | 146.36 | 152.56 | 4.23 | |
| 3 | \(K_I\) | 938.25 | 961.33 | 2.46 |
| \(K_{II}\) | 223.39 | 233.13 | 4.36 | |
| 5 | \(K_I\) | 1105.34 | 1139.61 | 3.10 |
| \(K_{II}\) | 298.74 | 312.87 | 4.73 |
4. Simulation Results and Discussion of Crack Propagation Path
4.1 Effect of Crack Nucleation Position
Five different crack nucleation points are selected on the tooth root fillet, labeled from 1 to 5 (from the tooth side to the rim side). The simulation results show that when the crack originates near the tooth root circle, the propagation path tends to bend toward the gear rim. Conversely, when the crack starts closer to the tooth flank, the path tends to remain inside the tooth and eventually leads to tooth fracture. The crack propagation path is initially directed toward the rim during single-tooth meshing, then gradually deflects in a parabolic curve during double-tooth meshing. The equivalent stress distribution indicates that the crack growth rate is higher in the early stage and decreases in the later stage.
4.2 Effect of Web Outer Diameter
Four web outer radii (40, 42, 46, 48 mm) are simulated. The results demonstrate that the crack path is deflected toward the gear rim as the web outer radius increases. For \(r_{wo} = 46\) mm and 48 mm, the crack paths for nucleation positions 4 and 5 obviously propagate into the rim, indicating a higher risk of rim fracture. Thus, a larger web outer diameter makes the rim more compliant and attracts the crack toward the rim.
4.3 Effect of Web Hole Position
The relative position between the web hole and the crack nucleation point is characterized by the angle \(\theta_P\), which ranges from -30° to 30°. The simulations show that the web hole exerts an attractive effect on the crack path. When the hole is located at \(\theta_P = 0°\) (directly below the crack), the influence on the crack deflection is most significant. For \(\theta_P = \pm 30°\), the influence is smaller but still noticeable. To quantify the deflection, a disturbance coefficient \(C\) is introduced:
$$
C = \frac{Y’}{X’}
$$
where \(X’\) and \(Y’\) are the horizontal and vertical coordinates of the point where the crack reaches the web hole boundary. Figure 3-16 in the original paper suggests that the disturbance coefficient increases as the angle changes from -30° to 30°, especially for nucleation positions 4 and 5.
4.4 Stress Intensity Factor Analysis
The SIFs are computed using the J-integral method. The results show that \(K_I\) is more than two orders of magnitude larger than \(K_{II}\), indicating that mode I dominates the crack growth in spur gear root cracks. \(K_I\) increases in an exponential-like manner with crack depth, while \(K_{II}\) remains relatively stable for small cracks. The variation of \(K_I\) with crack depth for different nucleation positions is shown in the following table:
| Depth (mm) | Pos 1 | Pos 2 | Pos 3 | Pos 4 | Pos 5 |
|---|---|---|---|---|---|
| 1 | 868.6 | 905.3 | 961.3 | 1042.7 | 1139.6 |
| 2 | 1024.5 | 1088.2 | 1156.4 | 1210.5 | 1285.3 |
| 3 | 1190.8 | 1262.1 | 1340.9 | 1398.2 | 1467.5 |
| 4 | 1356.2 | 1441.5 | 1523.7 | 1590.4 | 1658.9 |
It is observed that at positions 4 and 5, the SIF changes more obviously with crack depth, which aligns with the crack path deflection toward the rim. The web outer diameter and the web hole position also affect the SIF values. Larger web outer radii lead to higher SIF, especially for deep cracks.
5. Time-Varying Meshing Stiffness of Spur Gears
5.1 Potential Energy Method for Healthy Gears
The potential energy method treats a gear tooth as a cantilever beam with varying cross-section. The total meshing stiffness consists of Hertzian stiffness, bending stiffness, shear stiffness, axial compressive stiffness, and gear foundation stiffness. The formulas are given as follows:
Hertzian stiffness:
$$
K_h = \frac{\pi E L}{4(1-\nu^2)}
$$
Bending stiffness reciprocal:
$$
\frac{1}{K_b} = \int_{-\alpha_1}^{\alpha_2} \frac{12\left[1 + \cos\alpha_1(\alpha_2-\alpha)\sin\alpha – \cos\alpha\right]^2(\alpha_2-\alpha)\cos\alpha}{E L\left[\sin\alpha + (\alpha_2-\alpha)\cos\alpha\right]^3} d\alpha
$$
Shear stiffness reciprocal:
$$
\frac{1}{K_s} = \int_{-\alpha_1}^{\alpha_2} \frac{1.2(1+\nu)(\alpha_2-\alpha)\cos\alpha\cos^2\alpha_1}{E L\left[\sin\alpha + (\alpha_2-\alpha)\cos\alpha\right]} d\alpha
$$
Axial compressive stiffness reciprocal:
$$
\frac{1}{K_a} = \int_{-\alpha_1}^{\alpha_2} \frac{(\alpha_2-\alpha)\cos\alpha\sin^2\alpha_1}{2 E L\left[\sin\alpha + (\alpha_2-\alpha)\cos\alpha\right]} d\alpha
$$
Foundation stiffness reciprocal:
$$
\frac{1}{K_f} = \frac{\cos^2\alpha_m}{E L}\left\{L^*\left(\frac{u_f}{S_f}\right)^2 + M^*\left(\frac{u_f}{S_f}\right) + P^*\left(1 + Q^*\tan^2\alpha_m\right)\right\}
$$
where \(L^*\), \(M^*\), \(P^*\), \(Q^*\) are polynomial functions of \(u_f\) and \(\theta_f\) as defined in the original paper. The total mesh stiffness for one pair of teeth is:
$$
\frac{1}{K_e} = \frac{1}{K_h} + \frac{1}{K_{b1}} + \frac{1}{K_{s1}} + \frac{1}{K_{a1}} + \frac{1}{K_{f1}} + \frac{1}{K_{b2}} + \frac{1}{K_{s2}} + \frac{1}{K_{a2}} + \frac{1}{K_{f2}}
$$
5.2 Cracked Tooth Stiffness
When a crack appears at the tooth root, the cross-sectional area and moment of inertia are reduced. Depending on the crack depth \(q_i\) and the crack angle \(\theta_i\), two cases are considered.
Case 1 (shallow crack): \(u_c < u_0\) or \(u_c \ge u_0\) but \(\alpha_1 \le \alpha_c\). The bending stiffness reciprocal is:
$$
\frac{1}{K_b’} = \int_{-\alpha_1}^{\alpha_2} \frac{12\left[1+\cos\alpha_1(\alpha_2-\alpha)\sin\alpha-\cos\alpha\right]^3(\alpha_2-\alpha)\cos\alpha}{E L\left[\sin\alpha_2 – \frac{q_i}{R_{b1}}\sin\theta_i + \sin\alpha + (\alpha_2-\alpha)\cos\alpha\right]^3} d\alpha
$$
and the shear stiffness reciprocal is:
$$
\frac{1}{K_s’} = \int_{-\alpha_1}^{\alpha_2} \frac{2.4(1+\nu)(\alpha_2-\alpha)\cos\alpha\cos^2\alpha_1}{E L\left[\sin\alpha_2 – \frac{q_i}{R_{b1}}\sin\theta_i + \sin\alpha + (\alpha_2-\alpha)\cos\alpha\right]} d\alpha
$$
Case 2 (deep crack): \(u_c \ge u_0\) and \(\alpha_1 > \alpha_c\). The integrals are split into two parts, as shown in the original equations (4.24) and (4.25). In this paper, the calculations are performed for the gear pair whose parameters are listed in Table 4.5. The crack depths and angles used in the simulation are summarized in the following tables.
| Case | Healthy | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 |
|---|---|---|---|---|---|---|---|---|---|---|
| Depth (mm) | 0 | 0.6 | 1.0 | 1.4 | 2.0 | 2.5 | 3.0 | 3.6 | 4.0 | – |
| Case | Healthy | 1 | 2 | 3 | 4 | 5 | 6 | 7 |
|---|---|---|---|---|---|---|---|---|
| Angle (deg) | – | 30 | 45 | 60 | 75 | 90 | 105 | – |
| Parameter | Driving gear | Driven gear |
|---|---|---|
| Number of teeth \(z\) | 25 | 53 |
| Module \(m\) (mm) | 3 | 3 |
| Face width \(B\) (mm) | 40 | 40 |
| Pressure angle \(\alpha\) (deg) | 20 | 20 |
| Poisson’s ratio \(\nu\) | 0.3 | 0.3 |
| Elastic modulus \(E\) (GPa) | 210 | 210 |
| Base circle radius \(R_b\) (mm) | 35.2 | 74.7 |
| Mass \(m\) (kg) | 1.45 | 2.08 |
| Moment of inertia \(I\) (kg·m²) | 1.33×10⁻⁴ | 2.42×10⁻⁴ |
| Support stiffness \(k_x\), \(k_y\) (N/m) | 1.8×10⁸ | 1.8×10⁸ |
| Support damping \(c_x\), \(c_y\) (N·s/m) | 6.56×10⁴ | 6.56×10⁴ |
| Mesh damping \(c_m\) (N·s/m) | 67 | 67 |
| Load torque \(T\) (Nm) | 40 | 40 |
| Speed (rpm) | 800 | – |
5.3 Results of TVMS
The calculated TVMS curves are shown in the original paper. The key findings are listed below:
- The inverse sum of bending and shear stiffness increases in an exponential-like manner with crack depth.
- The single-tooth mesh stiffness decreases with increasing crack depth. The reduction is small for shallow cracks but becomes significant for deep cracks, especially in the early and late meshing regions.
- The combined mesh stiffness also decreases with crack depth, with a more pronounced reduction during double-tooth engagement.
- Increasing the crack angle also reduces the TVMS, but the effect is much smaller than that of the crack depth.
The following table gives the computed combined mesh stiffness (in N/m) at the pitch point for different crack depths.
| Crack depth (mm) | Healthy | 1 | 2 | 3 | 4 |
|---|---|---|---|---|---|
| TVMS (10⁸ N/m) | 2.83 | 2.72 | 2.58 | 2.41 | 2.20 |
6. Dynamic Analysis of Cracked Spur Gear System
6.1 Six-Degree-of-Freedom Dynamic Model
A 6-DOF lumped parameter model is established to investigate the vibration characteristics of the spur gear transmission system. The model includes two rotational degrees of freedom (torsional motions \(\theta_1\), \(\theta_2\)) and four translational degrees of freedom (\(x_1, y_1\) for the driving gear and \(x_2, y_2\) for the driven gear). The equations of motion are given by:
$$
\begin{cases}
m_1 \ddot{x}_1 + c_{x1}\dot{x}_1 + k_{x1}x_1 – F_f = 0\\[4pt]
m_2 \ddot{x}_2 + c_{x2}\dot{x}_2 + k_{x2}x_2 – F_f = 0\\[4pt]
m_1 \ddot{y}_1 + c_{y1}\dot{y}_1 + k_{y1}y_1 + (k_m y + c_m \dot{y}) = 0\\[4pt]
m_2 \ddot{y}_2 + c_{y2}\dot{y}_2 + k_{y2}y_2 – (k_m y + c_m \dot{y}) = 0\\[4pt]
I_1 \ddot{\theta}_1 – R_{b1}(k_m y + c_m \dot{y}) – M_1 = T_1\\[4pt]
I_2 \ddot{\theta}_2 + R_{b2}(k_m y + c_m \dot{y}) – M_2 = -T_2
\end{cases}
$$
where \(y = y_1 – y_2 + R_{b1}\theta_1 – R_{b2}\theta_2\) is the relative displacement along the line of action. The frictional force \(F_f\) is:
$$
F_f = \mu F_n
$$
with the normal meshing force:
$$
F_n = k_m\left[(y_1 – y_2) – (R_{b1}\theta_1 – R_{b2}\theta_2)\right] + c_m\left[(\dot{y}_1 – \dot{y}_2) – (R_{b1}\dot{\theta}_1 – R_{b2}\dot{\theta}_2)\right]
$$
The system equations are solved using the fourth-order Runge-Kutta method. The gear mesh frequency is calculated as:
$$
f_m = f_{01} z_1 = f_{02} z_2
$$
For the experiment, the driving speed is 800 rpm, the first-stage teeth numbers are 23 and 39, so the intermediate shaft speed is \(800 \times 23/39 = 471.79\) rpm, giving \(f_0 = 7.86\) Hz. The mesh frequency of the second-stage gear pair is \(f_m = 7.86 \times 39 = 306.54\) Hz (approximately 307 Hz).
6.2 Experimental Setup
A standard two-stage spur gearbox is used in the experiment. The gearbox is driven by a three-phase motor and loaded by a magnetic powder brake. Vibration acceleration signals are measured using a PCB 353B18 accelerometer mounted on the bearing housing near the intermediate shaft. A LMS data acquisition system is used to record the signals. The sampling frequency is 10240 Hz and the sampling duration is 30 s. Cracks are machined by wire electrical discharge machining (WEDM) with different depths: 1 mm, 2 mm, 3 mm, and 4 mm, all with a crack angle of 75°. The test gear is the driven gear of the second stage with 53 teeth.
6.3 Vibration response analysis
The time-domain acceleration response of the driving gear in the \(y\)-direction is obtained from the dynamic model. For the healthy gear, the response is relatively steady with an amplitude range of about -5 to 5 m/s². For a cracked gear with a 2 mm crack, the response shows periodic impacts with an amplitude range of about -80 to 80 m/s². The frequency spectrum of the cracked gear exhibits sidebands around the mesh frequency and its harmonics due to the periodic modulation caused by the cracked tooth. The envelope spectrum obtained after band-pass filtering and Hilbert transform reveals a clear peak at the shaft rotational frequency \(f_0 = 7.86\) Hz and its multiples. The amplitudes of the rotational frequency for different crack depths are summarized in Table below:
| Crack depth (mm) | Rotational frequency amplitude | Envelope spectrum entropy |
|---|---|---|
| 1 | 64.65 | 7.96 |
| 2 | 75.47 | 7.93 |
| 3 | 96.71 | 7.80 |
| 4 | 114.14 | 7.77 |
It is observed that as the crack depth increases, the amplitude at the rotational frequency increases, while the information entropy of the envelope spectrum decreases. This indicates that a smaller entropy corresponds to a more severe fault. The sidebands around the mesh frequency also become more abundant and stronger for deeper cracks, which is a typical sign of localized gear tooth damage.
7. Fault Feature Extraction using Adaptive VMD-KSVD
7.1 Variational Mode Decomposition
VMD decomposes a signal into \(K\) intrinsic mode functions (IMFs) with limited bandwidth. The constrained variational problem is:
$$
\min_{\{u_k\},\{\omega_k\}} \left\{ \sum_k \left\| \partial_t \left[ \left(\delta(t) + \frac{j}{\pi t}\right) * u_k(t) \right] e^{-j\omega_k t} \right\|_2^2 \right\}
$$
$$
\text{s.t.} \quad \sum_k u_k(t) = f(t)
$$
Here, \(u_k\) and \(\omega_k\) are the \(k\)-th component and its center frequency. The augmented Lagrangian is iteratively solved via the alternating direction method of multipliers (ADMM). The frequency-domain update equation is:
$$
\hat{u}_k^{n+1}(\omega) = \frac{\hat{f}(\omega) – \sum_{i\neq k}\hat{u}_i(\omega) + \hat{\lambda}(\omega)/2}{1 + 2\alpha(\omega – \omega_k)^2}
$$
where \(\alpha\) is the penalty factor. The center frequency is updated using:
$$
\omega_k^{n+1} = \frac{\int_0^\infty \omega |\hat{u}_k(\omega)|^2 d\omega}{\int_0^\infty |\hat{u}_k(\omega)|^2 d\omega}
$$
7.2 K-Singular Value Decomposition
The KSVD algorithm is a dictionary learning method that alternates between sparse coding and dictionary updating. Using orthogonal matching pursuit (OMP) for sparse coding and SVD for updating the dictionary atoms, KSVD can effectively denoise the signal while preserving transient features. The flow of KSVD is shown in the original paper. In this study, KSVD is first applied to the raw gear vibration signal to reduce noise and enhance the fault-related transients.
7.3 Adaptive Parameter Selection for VMD
To avoid manually setting the number of modes \(K\) and the penalty factor \(\alpha\), an improved scale-space representation is employed. The spectrum of the denoised signal is convolved with a Gaussian kernel at different scales:
$$
L(m,t) = \sum_{n=-M}^{M} f(m-n) g(n;t)
$$
$$
g(n;t) = \frac{1}{\sqrt{2\pi t}} e^{-n^2/(2t)}
$$
The local minima of \(L(m,t)\) are tracked across scales, and the scale-space curves \(S_i = \{s | R_{i-1} \to R_i\}\) are generated. An Otsu threshold selects the number of significant minima, which determines the number of IMFs \(K\). The initial center frequencies and the weighted penalty factors are computed as:
$$
\hat{\omega}_k^{init} = \frac{\int_0^\infty \omega |\hat{f}_{B_k}(\omega)|^2 d\omega}{\int_0^\infty |\hat{f}_{B_k}(\omega)|^2 d\omega}
$$
$$
\alpha_k = \alpha \left(\frac{f_s – B_k}{f_s}\right)
$$
where \(B_k\) is the bandwidth of the \(k\)-th segment and \(f_s\) is the sampling frequency. This adaptive VMD is combined with the KSVD-preprocessed signal to extract the fault features.
7.4 Results of Feature Extraction
The proposed VMD-KSVD method is applied to the measured signals of healthy and cracked gears. The scale-space analysis indicates that the signal can be divided into three significant frequency bands for all crack depths. The corresponding center frequencies are close to the mesh frequency and its second and third harmonics. The adaptive parameters for different crack depths are listed in the following table:
| Crack depth (mm) | IMF1 | IMF2 | IMF3 | |||
|---|---|---|---|---|---|---|
| \(\hat{\omega}_1\) (Hz) | \(\alpha_1\) | \(\hat{\omega}_2\) (Hz) | \(\alpha_2\) | \(\hat{\omega}_3\) (Hz) | \(\alpha_3\) | |
| 1 | 307.11 | 1921.87 | 614.22 | 1920.98 | 921.33 | 1960.94 |
| 2 | 306.72 | 1921.56 | 613.44 | 1922.12 | 920.16 | 1961.20 |
| 3 | 307.34 | 1920.96 | 614.68 | 1921.22 | 922.02 | 1960.90 |
| 4 | 306.41 | 1922.20 | 612.82 | 1921.45 | 919.23 | 1960.55 |
The decomposed IMFs in both time and frequency domains are analyzed. For the healthy gear, the IMF amplitudes are small, and the maximum amplitude in the frequency domain is about 0.4. In contrast, for the cracked gear, the IMF amplitudes are much larger and increase with crack depth. The maximum frequency-domain amplitudes for the first IMF are approximately 2.2, 3.3, 3.5, and 3.1 for crack depths of 1, 2, 3, and 4 mm, respectively. The third and fourth IMFs contain higher frequency components with smaller amplitudes. The sidebands around the mesh frequency become more visible as the crack propagates. This demonstrates that the adaptive VMD-KSVD method can effectively extract the fault features related to spur gear root cracks.
| Condition | Healthy | 1 mm | 2 mm | 3 mm | 4 mm |
|---|---|---|---|---|---|
| IMF1 amplitude | 0.4 | 2.2 | 3.3 | 3.5 | 3.1 |
8. Discussion
The results presented in this paper show that the improved XFEM is a powerful tool for simulating spur gear crack propagation. The crack path can be accurately predicted without remeshing, even when the gear has a complex blank with holes. The SIF analysis reveals that mode I dominates the crack propagation, which is consistent with the fact that the tooth root is primarily subjected to bending stress. The presence of the web and its outer diameter significantly alters the crack path. A larger web outer diameter reduces the rim stiffness and makes the crack more prone to propagate into the rim. The web hole also attracts the crack, causing an additional deflection. The disturbance coefficient \(C\) introduced in this paper provides a quantitative measure of the hole effect.
The TVMS calculation demonstrates that the stiffness reduction depends strongly on crack depth. This is because the cross-section at the crack location becomes smaller, reducing the bending and shear stiffness. The mesh stiffness is an input to the dynamic model, and its reduction directly leads to larger vibration amplitudes and impact features. The experimental results confirm that the envelope spectrum amplitude at the shaft rotational frequency increases with crack depth, and the information entropy decreases. These metrics can be used as indicators of spur gear crack severity.
The proposed adaptive VMD-KSVD algorithm overcomes the limitation of traditional VMD by automatically determining the decomposition parameters. The scale-space approach uses the spectral characteristics of the signal itself to estimate the number of modes and the penalty factors. Combined with the KSVD denoising, the fault features are extracted more reliably. The increase in spectral amplitude of the extracted IMFs with crack depth verifies the sensitivity of the proposed method to early crack faults.
9. Conclusion
This paper comprehensively investigates the crack propagation path and vibration characteristics of spur gears with root cracks. The following conclusions are drawn:
- The stress field near the tooth root crack tip in spur gears has a singular behavior, and the complex variable method effectively describes the stress and displacement distributions. The J-integral is suitable for computing the SIF because of its path independence.
- The improved XFEM with shifted Heaviside and ramp crack-tip enrichment functions, along with the hole-enrichment function, accurately simulates crack propagation in spur gears with webs and holes. The maximum error in SIF is below 5% compared to literature values.
- The crack nucleation position, web outer diameter, and web hole position significantly influence the crack propagation path. When the nucleation point is near the root circle, the crack tends to grow toward the rim. Increasing the web outer diameter also promotes rim fracture. The web hole acts as a crack attractor, and the disturbance coefficient reaches its maximum for holes near the root fillet.
- Mode I dominates the SIF for spur gear root cracks; \(K_I\) increases exponentially with crack depth, whereas \(K_{II}\) remains relatively small.
- The TVMS of a cracked spur gear decreases with increasing crack depth and crack angle. The effect of crack depth is more pronounced than that of the crack angle, especially for deep cracks.
- The 6-DOF dynamic model predicts that cracks induce periodic impacts in the time-domain response. The envelope spectrum amplitude at the shaft rotational frequency increases as the crack depth increases from 1 mm to 4 mm, while the envelope spectrum entropy decreases. These features can be used to assess the crack severity.
- The adaptive VMD-KSVD method automatically determines the number of IMFs and their corresponding center frequencies and penalty factors. It successfully extracts fault features from noisy gear vibration signals. The spectral amplitudes of the IMFs increase with crack depth, indicating that the method is effective for early spur gear root crack detection.
Future work will focus on applying the proposed methodology to other gear types, such as helical gears and planetary gears, and on developing more advanced diagnostic techniques for gear tooth spalling and pitting faults.