In this study, we focus on the critical transmission component of spiral bevel gears, which are widely used in high-power and extreme service conditions due to their excellent torque transmission capacity, high contact ratio, smooth operation, compact structure, and low noise. However, under high-speed and heavy-load conditions, spiral bevel gears are prone to crack faults, especially under sudden load changes and stress concentration. The complex geometric tooth surfaces, spatial curvature characteristics, and special local conjugate meshing processes of spiral bevel gears pose challenges for contact analysis and stiffness calculation. This paper systematically investigates the analytical calculation method for time-varying meshing stiffness under crack fault conditions, combining geometric design principles, manufacturing methods, and meshing transmission mechanisms.
We begin by deriving the tooth surface processing equations based on the geometric design and meshing mechanisms of spiral bevel gears. The hypothetical generating gear formed by the movement of the processing tool in space is completely conjugated with the machined gear, accurately establishing the complex tooth surface in real space. The tooth surface is represented by a parabolic blade profile, and the generating surface consists of two parts: the working surface formed by the rotation of the parabolic blade and the transition surface formed by the circular arc profile. The mathematical representation of the generating surface is given by:
$$ \mathbf{r}_g^{(a)}(s_g, \theta_g) = \begin{bmatrix}
\left[ R_g \pm (s_g + s_{g0}) \sin \alpha_g \pm a_c s_g^2 \cos \alpha_g \right] \cos \theta_g \\
\left[ R_g \pm (s_g + s_{g0}) \sin \alpha_g \pm a_c s_g^2 \cos \alpha_g \right] \sin \theta_g \\
– (s_g + s_{g0}) \cos \alpha_g + a_c s_g^2 \sin \alpha_g
\end{bmatrix} $$
where $s_g$ and $\theta_g$ are surface coordinates, $a_c$ is the parabola coefficient, $\alpha_g$ is the blade angle at point M, and $R_g$ is the radius of the blade tip. The unit normal vector of the generating surface is derived as:
$$ \mathbf{n}_g^{(a)}(s_g, \theta_g) = \frac{\mathbf{N}_g^{(a)}}{|\mathbf{N}_g^{(a)}|}, \quad \mathbf{N}_g^{(a)} = \frac{\partial \mathbf{r}_g^{(a)}}{\partial s_g} \times \frac{\partial \mathbf{r}_g^{(a)}}{\partial \theta_g} $$
The tooth surface processing equation is applied to establish the real spatial curved tooth surface, as shown in the projection below. The tooth surface is divided into the working part and the transition part, and the complete tooth profile is formed by rotating the projection.

Next, we employ the tooth contact analysis method to simulate the meshing process of the spiral bevel gear pair. The contact trajectory, contact ellipse, and transmission error curve are analyzed in detail. The meshing process requires continuous tangency between the pinion and gear tooth surfaces, represented by the position and normal vectors in the fixed coordinate system $S_h$. The equations for the pinion and gear tooth surfaces in $S_h$ are:
$$ \mathbf{r}_h^{(1)}(s_p, \theta_p, \psi_1, \phi_1) = \mathbf{M}_{hb1} \mathbf{M}_{b11}(\phi_1) \mathbf{r}_1(s_p, \theta_p, \psi_1) $$
$$ \mathbf{r}_h^{(2)}(s_g, \theta_g, \psi_2, \phi_2) = \mathbf{M}_{hb2} \mathbf{M}_{b22}(\phi_2) \mathbf{r}_2(s_g, \theta_g, \psi_2) $$
The tangency conditions are given by:
$$ \mathbf{r}_h^{(1)} – \mathbf{r}_h^{(2)} = 0, \quad \mathbf{n}_h^{(1)} – \mathbf{n}_h^{(2)} = 0 $$
with additional constraints from the meshing equations $f_{1p}(s_p, \theta_p, \psi_1) = 0$ and $f_{2g}(s_g, \theta_g, \psi_2) = 0$. The transmission error function is calculated as:
$$ \Delta \phi_2(\phi_1) = \phi_2(\phi_1) – \frac{N_1}{N_2} \phi_1 $$
The contact analysis results in a continuous contact area characterized by elongated contact ellipses, whose major and minor axes lengths and directions are key parameters determined by the principal curvatures and directions of the contact surfaces.
Based on the meshing analysis, we develop an analytical model for the time-varying meshing stiffness of spiral bevel gears using the slicing method and potential energy method. The gear tooth is divided into multiple slices along the tooth width, and the stiffness components for each slice are calculated. The total meshing stiffness is the sum of the stiffnesses of all engaged tooth pairs at each meshing instant. The potential energy includes contributions from the tooth root transition curve and the contact force acting point to the base circle.
The force analysis model for a healthy spiral bevel gear slice is shown below. The meshing force $\mathbf{F}_m$ is decomposed into tangential force $\mathbf{F}_t$, radial force $\mathbf{F}_r$, and axial force $\mathbf{F}_a$:
$$ \begin{cases}
F_t = F_m \cos(\alpha) \cos(\beta) \\
F_r = F_m \left( \sin(\alpha) \cos(\delta_G) – \cos(\alpha) \sin(\beta) \sin(\delta_G) \right) \\
F_a = F_m \left( \sin(\alpha) \sin(\delta_G) + \cos(\alpha) \sin(\beta) \cos(\delta_G) \right)
\end{cases} $$
where $\alpha$ is the pressure angle, $\beta$ is the spiral angle, and $\delta_G$ is the pitch cone angle. The relationship between the meshing force and the load torque $T$ is:
$$ F_m = \frac{T}{r_p \cos(\alpha) \cos(\beta)} $$
where $r_p$ is the radial distance from the meshing point to the rotation axis. The stiffness components for the $k$-th slice include axial compressive stiffness $k_{ka}$, shear stiffness $k_{ks}$, bending stiffness $k_{kb}$, Hertzian contact stiffness $k_{kh}$, and gear foundation stiffness $k_{kf}$. The expressions for these stiffnesses are:
$$ \frac{1}{k_{ka}} = \int_{-\alpha_{k1}}^{\alpha_{k2}} \frac{-\sin^2 \alpha_{k1} R_{kb}}{E A_{xk}} d\alpha_k + \frac{\sin^2 \alpha_{k1} d_{k1}}{E A_{dk1}} $$
$$ \frac{1}{k_{ks}} = \int_{-\alpha_{k1}}^{\alpha_{k2}} \frac{-1.2 \cos^2 \alpha_{k1} R_{kb}}{G A_{xk}} d\alpha_k + \frac{1.2 \cos^2 \alpha_{k1} d_{k1}}{G A_{dk1}} $$
$$ \frac{1}{k_{kb}} = \int_{-\alpha_{k1}}^{\alpha_{k2}} \frac{-\left[ (d_k – x_k) \cos \alpha_{k1} – h_k \sin \alpha_{k1} \right]^2 R_{kb}}{E I_{xk}} d\alpha_k + \int_0^{d_{k1}} \frac{\left[ (d_k + x_k) \cos \alpha_{k1} – h_k \sin \alpha_{k1} \right]^2}{E I_{dk1}} dx_k $$
$$ k_{kh} = \frac{\pi E \Delta w}{4(1 – \nu^2)}, \quad k_{kf} = \frac{F}{\delta_{kf}} $$
where $E$ is Young’s modulus, $G$ is the shear modulus, $\nu$ is Poisson’s ratio, $A_{xk}$ and $I_{xk}$ are the contact area and moment of inertia for the slice, and $\Delta w$ is the slice width. The meshing stiffness for the $k$-th slice is:
$$ k_{kt} = \frac{1}{\frac{1}{k_{kh}} + \frac{1}{k_{kb1(2)}} + \frac{1}{k_{ks1(2)}} + \frac{1}{k_{ka1(2)}} + \frac{1}{k_{kf1(2)}}} $$
The total meshing stiffness for a single tooth pair is the sum of the stiffnesses of all slices, and the comprehensive meshing stiffness accounts for multiple tooth pairs engaged simultaneously due to the contact ratio of 2-3:
$$ k_{\text{single}} = \sum_{k=1}^n k_{kt}, \quad k_{\text{total}} = \sum_{j=1}^m k_{\text{single}} $$
For cracked spiral bevel gears, the crack affects the shear and bending stiffnesses. The crack propagation is modeled along the tooth width or tooth tip direction, with parameters such as crack length $L$, initial crack depth $q_1$, termination crack depth $q_3$, crack angle $v$, and whether the crack exceeds the tooth centerline. The stiffness calculation for cracked slices is modified based on the crack geometry. For example, when the crack does not exceed the centerline, the crack depths for the $k$-th slice are:
$$ q_{k1} = \begin{cases}
\left( q_1 – (q_1 – q_3) \frac{k-1}{k} \right) \frac{N_c – k}{N_c}, & k \leq N_c \\
0, & k > N_c
\end{cases}, \quad q_{k2} = 0 $$
where $N_c$ is the number of cracked slices. The modified shear and bending stiffnesses are calculated by integrating over the cracked region, considering the reduced cross-sectional area and moment of inertia.
We analyze the time-varying meshing stiffness under various crack conditions, such as different crack lengths, depths, and angles. The results show that the stiffness decreases significantly as the crack propagates, and the reduction is more pronounced when the crack exceeds the tooth centerline. The table below summarizes the geometric design and processing parameters of the spiral bevel gear used in the analysis.
| Parameter | Pinion | Gear |
|---|---|---|
| Number of teeth | 32 | 53 |
| Module (mm) | 3 | 3 |
| Face width (mm) | 14 | 20 |
| Pitch cone angle (°) | 41.56 | 48.43 |
| Face cone angle (°) | 43.07 | 49.67 |
To validate the analytical model, we compare the calculated time-varying meshing stiffness with finite element analysis results. The finite element model simulates the meshing process, and the stiffness values are extracted from the contact forces and deformations. The comparison shows good agreement between the analytical and finite element results, confirming the accuracy of the proposed method.
Furthermore, we conduct experimental validation by measuring the vibration acceleration signals of a spiral bevel gear pair with a preset crack fault. The experimental setup includes an AC asynchronous motor, dynamic torque sensor, data acquisition card, and accelerometers. The spiral bevel gear is tested at speeds of 2000 rpm and 3000 rpm, and the acceleration signals are compared with those calculated using the analytical model with the same crack parameters. The time-domain and frequency-domain analyses reveal similar characteristics, such as periodic impacts and sideband modulation, in both the calculated and measured signals, verifying the effectiveness of the stiffness calculation method.
In conclusion, this study provides a comprehensive analytical framework for calculating the time-varying meshing stiffness of spiral bevel gears under crack faults. The key contributions include the derivation of tooth surface processing equations, detailed meshing contact analysis, development of stiffness models for healthy and cracked gears, and validation through finite element simulation and experiments. The results demonstrate the significant impact of crack parameters on meshing stiffness and offer insights for gear design optimization, fault diagnosis, and dynamic simulation. Future work will explore the application of this method to other gear types and more complex fault scenarios.
