Spiral bevel gears transmit motion and torque between intersecting shafts, usually at 90°, and are chosen where high load capacity, a large transmission ratio, reliable operation, and smooth meshing are required. Their tooth contact is oblique and progresses along the face width, so the mesh process combines time-varying stiffness, backlash, static transmission error, and load fluctuation. These features make the gear pair a nonlinear dynamical system: small changes in operating parameters can move the response from a periodic orbit to quasi-periodic or chaotic motion, with consequences for vibration, impact, and tooth contact.
Spiral bevel gears used for right-angle power transmission.
Much of the established nonlinear gear literature concerns spur gears and planetary sets. Spiral bevel gears are geometrically more complex, and their dynamic assessment has relied mainly on phase-space reconstruction, Poincaré maps, and Lyapunov exponents. Recurrence plots and recurrence quantification analysis (RQA) offer a complementary route: they preserve the visual structure of repeated states and return scalar measures of diagonal and vertical line organization. The analysis below uses that route to examine three control parameters—dimensionless meshing frequency \(\omega_h\), stiffness coefficient \(\varepsilon\), and damping coefficient \(\xi_m\)—and to identify operating regions where the response is simple and stable.
Model Structure and Nonlinear Sources
The gear pair is represented as two rigid bodies, each with three translations and one rotation, giving eight degrees of freedom. The normal relative displacement at the mesh point is
\[
\lambda_n=(X_1-X_2)c_1-(Y_1-Y_2)c_2-(Z_1-Z_2+r_1\theta_1-r_2\theta_2)c_3-e_n(t),
\]
where \(r_1\) and \(r_2\) are mesh-point radii, \(e_n(t)\) is the normal static transmission error, and \(c_1\)–\(c_3\) project the coordinate motion onto the mesh normal. The transmission error is represented by a harmonic series,
\[
e_n(t)=\sum_{\ell=1}^{N_e} A_{e\ell}\cos(\ell\Omega_h t+\Phi_{e\ell}),
\]
with meshing frequency \(\Omega_h\), harmonic amplitude \(A_{e\ell}\), and phase \(\Phi_{e\ell}\). The normal mesh force is
\[
F_n=k_h(t)f(\lambda_n)+c_h\dot{\lambda}_n,
\]
and its coordinate components are
\[
\begin{aligned}
F_x&=-F_n(\sin\alpha_n\sin\beta_m+\cos\alpha_n\sin\beta_m\sin\delta_1),\\
F_y&=F_n(\sin\alpha_n\sin\delta_1-\cos\alpha_n\sin\beta_m\cos\delta_1),\\
F_z&=F_n\cos\alpha_n\sin\beta_m.
\end{aligned}
\]
The projection coefficients are \(c_1=\cos\delta_1\sin\alpha_n\), \(c_2=\cos\delta_1\cos\alpha_n\sin\beta_m\), \(c_3=\cos\alpha_n\cos\beta_m\), \(c_4=\sin\alpha_n\cos\beta_m+\cos\alpha_n\sin\beta_m\sin\delta_1\), \(c_5=\sin\alpha_n\sin\delta_1-\cos\alpha_n\sin\beta_m\cos\delta_1\), and \(c_6=\cos\alpha_n\sin\beta_m\). The time-varying mesh stiffness and the backlash function are
\[
k_h(t)=k_m+\sum_{\ell=1}^{N_k}A_{k\ell}\cos(\ell\Omega_h t+\Phi_{k\ell}),
\]
\[
f(\lambda_n)=
\begin{cases}
\lambda_n-b_m, & \lambda_n \gt b_m,\\
\lambda_n, & |\lambda_n|\le b_m,\\
\lambda_n+b_m, & \lambda_n \lt -b_m,
\end{cases}
\]
where \(k_m\) is the mean stiffness and \(b_m\) is half the average normal backlash. After nondimensionalization, the response is governed by the dimensionless time \(\tau=\omega_n t\), frequency ratio \(\omega_h=\Omega_h/\omega_n\), stiffness ratio \(\tilde{k}_h=k_h/k_m\), and clearance \(b=b_m/b_c\). Only the fundamental harmonic is retained in the numerical case: \(\tilde{k}_h=1+0.2\cos(\omega_h\tau)\) and the dimensionless error excitation is \(0.2\omega_h^2\cos(\omega_h\tau)\). The baseline parameter set is \(\beta_m=35^\circ\), \(\alpha_n=20^\circ\), \(\xi_{ij}=0.01\), \(\xi_{h1}=\xi_{h2}=0.0125\), \(\xi_m=0.05\), \(f_{pv}=0\), \(f_{pm}=0.5\), \(k_{ij}=1.25\), and \(k_{h1}=k_{h2}=0.5\). Time series are obtained by numerical integration of the dimensionless equations.
Recurrence Diagnostics
Phase-space reconstruction converts a scalar response into a trajectory that preserves the essential dynamics. For delay \(\tau_d\) and embedding dimension \(m\), the reconstructed vector is
\[
\mathbf{X}_t=[x(t),x(t+\tau_d),\ldots,x(t+(m-1)\tau_d)].
\]
The delay is selected from the autocorrelation function \(C(\tau_d)\), using the first drop to \(C(0)/e\), which opens the attractor without making successive coordinates independent. The embedding dimension is selected by the false-nearest-neighbor method: when \(m\) is too small, points that are distant in the full state appear close after projection; the fraction of such false neighbors falls sharply once \(m\) is sufficient. These settings provide the trajectory used for recurrence analysis.
A recurrence plot marks pairs of states that fall within a threshold \(\epsilon_r\):
\[
R_{ij}=\Theta\left(\epsilon_r-\|\mathbf{X}_i-\mathbf{X}_j\|\right),\qquad i,j=1,\ldots,N,
\]
where \(\Theta(s)=1\) for \(s\ge 0\) and \(\Theta(s)=0\) for \(s \lt 0\). A marked pair appears as a black point at \((i,j)\), and the plot is symmetric about the main diagonal. The threshold is adjusted so that the recurrence rate is approximately \(0.02\), within the \(0.01\)–\(0.05\) range commonly used for comparable plots. This keeps point density from dominating visual differences between operating conditions.
Three RQA measures are used here. The recurrence rate is
\[
RR=\frac{1}{N^2}\sum_{i,j=1}^{N}R_{ij}.
\]
The recurrence entropy is the Shannon entropy of the diagonal-line length distribution,
\[
ENTR=-\sum_{l=l_{\min}}^{N}p(l)\ln p(l),
\]
where \(p(l)\) is the probability density of diagonal lines of length \(l\). Higher \(ENTR\) indicates a broader range of recurrence times and greater dynamical complexity. The trapping time is the mean length of vertical lines,
\[
TT=\frac{\sum_{v=v_{\min}}^{N}vP(v)}{\sum_{v=v_{\min}}^{N}P(v)},
\]
where \(P(v)\) is the number of vertical lines of length \(v\). \(TT\) measures how long the trajectory remains in a recurrent state. In plot structure, diagonal lines parallel to the main diagonal indicate deterministic periodic recurrence; segmented diagonals indicate period-2 or multi-period motion; vertical lines indicate laminar states; blocks indicate nonstationarity; and complex patterned recurrence indicates chaotic folding.
Effect of Meshing Frequency
Table 1 gives the recurrence measures for eight dimensionless meshing frequencies at \(\varepsilon=0.2\) and \(\xi_m=0.05\). The recurrence rate remains close to \(0.02\) throughout, so the differences are carried by line organization rather than by point density.
| \(\omega_h\) | RR | ENTR | TT |
|---|---|---|---|
| 0.30 | 0.01992 | 5.6903 | 23.7937 |
| 0.55 | 0.02009 | 4.0923 | 13.0949 |
| 0.70 | 0.02015 | 5.9160 | 10.3184 |
| 0.90 | 0.01904 | 4.8286 | 15.4563 |
| 1.30 | 0.02085 | 5.0137 | 5.7481 |
| 1.65 | 0.02066 | 5.9390 | 27.6131 |
| 1.85 | 0.02071 | 4.6106 | 4.0128 |
| 2.00 | 0.02031 | 4.5170 | 3.6408 |
At \(\omega_h=0.30\), \(0.55\), \(0.70\), \(0.90\), and \(1.30\), the recurrence plots are dominated by diagonal lines parallel to the main diagonal. The spacing changes with frequency, so the motions are periodic but have different periods. At \(0.70\), segmented short and long diagonal lines indicate period-2 motion; the attractor is a narrow closed loop with an amplitude range of \(1.92\)–\(1.96\). At \(0.90\), the attractor is distorted but remains periodic, with an amplitude range of \(1.8\)–\(2.0\). At \(1.30\), the attractor indicates multi-period motion and the amplitude range contracts to \(1.85\)–\(2.0\).
The response changes sharply at \(\omega_h=1.65\). The attractor folds and crosses, spans approximately \(-2\) to \(5\), and the recurrence plot develops a complex regular pattern rather than a set of simple parallel lines. This case has the largest \(ENTR\) (\(5.9390\)) and the largest \(TT\) (\(27.6131\)) in Table 1, consistent with chaotic motion. At \(1.85\) and \(2.00\), long diagonal lines with different spacing return, indicating quasi-periodic motion; \(ENTR\) falls to \(4.6106\) and \(4.5170\), and \(TT\) falls to \(4.0128\) and \(3.6408\).
\(TT\) alone is not a chaos indicator. At \(\omega_h=0.30\), \(TT=23.7937\) is large even though the motion is periodic, because the trajectory remains in a recurrent state for a long time. The diagnostic value comes from reading \(TT\), \(ENTR\), and the recurrence pattern together. For design, the interval around \(\omega_h=1.65\) is the least favorable in this dimensionless map and should be avoided when the actual speed can be scaled to that range.
Effect of Mesh Stiffness Coefficient
Table 2 varies the stiffness coefficient \(\varepsilon\) at \(\omega_h=1.65\) for two damping levels, \(\xi_m=0.05\) and \(\xi_m=0.07\). The response is non-monotonic: periodic, chaotic, and multi-period windows alternate as stiffness rises.
| \(\xi_m\) | \(\varepsilon\) | RR | ENTR | TT |
|---|---|---|---|---|
| 0.05 | 0.1 | 0.0203 | 4.9762 | 4.4200 |
| 0.05 | 0.3 | 0.0206 | 4.9004 | 17.8529 |
| 0.05 | 0.5 | 0.0197 | 5.1669 | 44.6562 |
| 0.05 | 0.6 | 0.0196 | 5.3563 | 19.1473 |
| 0.05 | 0.8 | 0.0203 | 4.0895 | 4.3997 |
| 0.05 | 1.0 | 0.0192 | 3.9290 | 4.1642 |
| 0.07 | 0.1 | 0.0204 | 4.0680 | 4.4279 |
| 0.07 | 0.3 | 0.0209 | 4.8913 | 18.1799 |
| 0.07 | 0.5 | 0.0210 | 5.0877 | 18.2166 |
| 0.07 | 0.6 | 0.0204 | 5.2858 | 38.3728 |
| 0.07 | 0.8 | 0.0208 | 3.9940 | 4.5159 |
| 0.07 | 1.0 | 0.0201 | 4.4031 | 4.3746 |
At \(\xi_m=0.05\), \(\varepsilon=0.1\) and \(0.3\) produce single closed curves. At \(\varepsilon=0.5\), the attractor folds and entangles, and the recurrence plot contains short diagonal segments with varying spacing—a chaotic signature. At \(\varepsilon=0.6\), vertical lines appear in the recurrence plot, indicating laminar states and multi-period motion. At \(\varepsilon=0.8\) and \(1.0\), the attractor returns to a single curve, and the amplitude range is smallest at \(\varepsilon=1.0\). The quantitative peak is split: \(TT\) reaches \(44.6562\) at \(\varepsilon=0.5\), while \(ENTR\) reaches \(5.3563\) at \(\varepsilon=0.6\).
At \(\xi_m=0.07\), \(\varepsilon=0.3\) produces a strongly nonstationary recurrence plot: the upper-right region is dense with black points while the lower-left is largely white, a drift pattern rather than a stationary attractor. At \(\varepsilon=0.5\), vertical lines indicate multi-period motion. At \(\varepsilon=0.6\), scattered points and blocks coexist, and both \(TT\) (\(38.3728\)) and \(ENTR\) (\(5.2858\)) reach their maxima for this damping level. Stable periodic motion returns at \(\varepsilon=0.1\), \(0.8\), and \(1.0\).
For both damping levels, the high-stiffness end \(\varepsilon=0.8\)–\(1.0\) is consistently simple: \(ENTR\) falls to about \(3.9\)–\(4.4\) and \(TT\) to about \(4.2\)–\(4.5\). The route is not monotonic—the system can pass through chaotic and multi-period windows before returning to a periodic state—but the high-stiffness region is the more stable one in this model.
Effect of Damping Coefficient
Table 3 varies the damping coefficient \(\xi_m\) at \(\omega_h=1.65\) for \(\varepsilon=0.2\) and \(\varepsilon=0.5\). Increasing damping from a low value does not immediately simplify the motion; complexity first rises to a mid-range peak and then falls at the highest damping level.
| \(\varepsilon\) | \(\xi_m\) | RR | ENTR | TT |
|---|---|---|---|---|
| 0.2 | 0.01 | 0.0202 | 4.7725 | 17.5770 |
| 0.2 | 0.05 | 0.0200 | 4.9910 | 29.3635 |
| 0.2 | 0.07 | 0.0196 | 5.7224 | 18.6100 |
| 0.2 | 0.10 | 0.0211 | 4.0549 | 4.5733 |
| 0.5 | 0.01 | 0.0192 | 4.6865 | 19.8113 |
| 0.5 | 0.05 | 0.0197 | 5.1669 | 44.6562 |
| 0.5 | 0.07 | 0.0190 | 5.2409 | 38.8932 |
| 0.5 | 0.10 | 0.0204 | 5.0867 | 17.7056 |
At \(\varepsilon=0.2\), \(\xi_m=0.05\) and \(0.07\) produce more complex attractors and recurrence plots, with vertical black lines appearing in the plots; the \(0.07\) case has the larger black area and the highest \(ENTR\) (\(5.7224\)). At \(\xi_m=0.10\), the attractor becomes a single closed curve, the amplitude range shrinks from about \(-3\) to \(5\) down to \(1\)–\(3\), and both measures fall to \(ENTR=4.0549\) and \(TT=4.5733\).
At \(\varepsilon=0.5\), the same non-monotonic pattern appears. At \(\xi_m=0.05\), the attractor is twisted and folded, and the recurrence plot is complex; \(TT\) reaches \(44.6562\), the largest value in Table 3. At \(\xi_m=0.07\), the recurrence plot shows a regular pattern with vertical lines, consistent with multi-period motion; \(ENTR\) reaches \(5.2409\). At \(\xi_m=0.01\) and \(0.10\), long diagonal lines dominate, indicating a return to periodic motion; at \(0.10\), \(TT\) is \(17.7056\) and \(ENTR\) is \(5.0867\). The high-damping case is therefore simpler than the mid-range peak, but it does not erase all small-scale recurrence structure at \(\varepsilon=0.5\).
The practical reading is that a mid-range damping coefficient can enlarge the chaotic window, while a sufficiently high damping coefficient restores a simple attractor. The two RQA measures also peak at different damping values—\(TT\) near \(0.05\) and \(ENTR\) near \(0.07\)—so neither scalar alone captures the transition.
Operating Guidelines and Scope
Within this dimensionless model, stable operation of spiral bevel gears is favored by three conditions: a stiffness coefficient near the high end of the tested range, a damping coefficient near the high end, and a meshing frequency away from the complex region around \(\omega_h=1.65\). At high stiffness and high damping, the reconstructed attractor is close to a single closed curve, the recurrence plot is dominated by long lines parallel to the main diagonal, and both \(ENTR\) and \(TT\) are small. Under the opposite conditions, the attractor folds and crosses, the recurrence plot becomes complex, and the response acquires chaotic characteristics.
These are trends from a lumped-parameter, fundamental-harmonic model, not absolute vibration levels for a specific gearbox. The gears are treated as rigid bodies; only the first harmonic of mesh stiffness and transmission error is retained; friction and thermal effects are outside the formulation. The results therefore define a dimensionless operating map that must be scaled to the tooth count, speed range, and load spectrum of a given spiral bevel gear pair before it is used as a forbidden-speed or damping-selection chart. Within that scope, recurrence plots and RQA provide a compact way to distinguish periodic, multi-period, quasi-periodic, and chaotic response regimes and to connect them to measurable design parameters.