Power-combining spiral bevel gear systems use two input pinions to drive a single output gear through two meshes. The arrangement is compact, suits intersecting shafts, and can transmit high torque at high speed, which makes it attractive for helicopter main reductions, marine propulsion, and wind-turbine power trains. The two meshes also make the dynamics less forgiving: time-varying mesh stiffness, transmission error, backlash, and damping act together, and a small change in one parameter can move the system between periodic, multi-periodic, and chaotic response. The practical question is therefore not only whether one operating point is stable, but where stable windows exist in parameter space, how wide they are, and which parameter combinations form their boundaries.
System architecture and nonlinear model
The model treats the two input pinions and the output gear as rigid bodies. Each body has three translational coordinates and one torsional coordinate. The output gear is denoted g0, the pinions g1 and g2, and the two meshes M1 and M2. The axes are orthogonal, and the two gear pairs are geometrically symmetric. Tooth friction is neglected, but backlash and combined transmission error are included at each mesh. The pinions have 23 teeth and the output gear has 78; both pinions are left-hand, the output gear is right-hand, the normal module is 4.25 mm, the normal pressure angle is 20°, and the mean spiral angle is 35°.
The mesh force along the line of action is written as
$$F_{0i}=k_i(t)x_{0i}+c_i\dot{x}_{0i},\qquad i=1,2$$
where x0i is the relative displacement at mesh i, and the mesh stiffness varies periodically:
$$k_i(t)=k_m\left[1+\varepsilon\sin(\Omega t+\psi_i)\right]$$
Here km is the mean mesh stiffness, ε is the stiffness fluctuation coefficient, Ω is the meshing frequency, and ψi is the mesh phase, set to zero in the initial phase convention. The combined transmission error is introduced as
$$e_i(t)=-\hat{E}\sin(\Omega t+\psi_i)$$
with Ê the error amplitude. The relative mesh displacement combines the translational and torsional projections of the three bodies and subtracts the transmission error. After the rigid-body rotation is eliminated through relative coordinates, the formulation has 11 generalized coordinates: nine translational coordinates and two relative mesh displacements.
| Component | Parameter | Value |
|---|---|---|
| Input pinions g1, g2 | Tooth number | 23 |
| Input torque | 400 N·m | |
| Moment of inertia | 0.018 kg·m² | |
| Support stiffness | 2.08×108 N/m | |
| Support damping | 558.12 N·s/m | |
| Output gear g0 | Tooth number | 78 |
| Output torque | 2713 N·m | |
| Moment of inertia | 0.042 kg·m² | |
| Support stiffness | 7.83×108 N/m | |
| Support damping | 651.02 N·s/m | |
| Meshes M1, M2 | Normal module | 4.25 mm |
| Normal pressure angle | 20° | |
| Mean spiral angle | 35° | |
| Mean mesh stiffness | 7.79×108 N/m | |
| Mesh damping | 729.86 N·s/m | |
| Stiffness fluctuation coefficient ε | 0.1–0.4 | |
| Mesh damping ratio ξ | 0.01–0.45 |
For numerical work, time is scaled as τ = ωnt, meshing frequency as ω = Ω/ωn, and displacement by a characteristic length. The equivalent mass of gear i is mi′ = Ji/ri², and the reference frequency is
$$\omega_n=\left[k_{m1}\left(\frac{1}{m_0′}+\frac{1}{m_1′}\right)\right]^{1/2}$$
The resulting dimensionless equations are integrated with a variable-step Runge–Kutta method. The model is deterministic and assumes rigid gear bodies, symmetric meshes, sinusoidal stiffness fluctuation, and no tooth friction; those assumptions define the scope of the parameter maps below.
Mapping solution domains in a two-parameter plane
A one-dimensional bifurcation diagram shows how a response changes along a single parameter cut, but it can miss interactions between two parameters. A two-parameter domain map treats the parameter plane itself as the object of study. The plane is divided into uniform rectangular cells, and each cell is represented by a parameter pair. For every cell, the governing equations are integrated to a steady attractor, the attractor is sampled on a Poincaré section, and the sample sequence is classified by its return period. A cell is assigned to a period-n domain when the iterates close after n returns; a sequence that does not close is classified as chaotic.
The classification uses a point-set distance on the Poincaré section. If the maximum norm between two attractor point sets is below ε1 = 1.0×10−2, the two cells are treated as belonging to the same solution domain. Cells whose neighbors carry a different label form the boundary set. These boundary cells approximate bifurcation sets: they mark parameter combinations where the attractor changes, a periodic window opens or closes, or a chaotic region begins.
In the maps discussed here, each axis is divided into 401 intervals, giving 401² = 160,801 cells per plane. The labels P1, P2, P3, P4, P8, and P16 denote period-1, period-2, period-3, period-4, period-8, and period-16 attractors. “Chaotic” denotes cells whose Poincaré iterates did not settle to a periodic closure within the iteration limit. The labels are numerical classifications over a specified parameter box, not analytical proofs of global stability.
Meshing-frequency bifurcation under error excitation
The interaction between meshing frequency and transmission-error amplitude is visible along a one-dimensional cut. At dimensionless error amplitude E = 28, three bifurcation points appear at ω = 0.917, 0.935, and 0.951. Over ω ∈ [0.9, 0.954], the response remains periodic, but the period-2 branch shows a jump-like rupture: part of the branch contracts without a gradual precursor, and the system maintains two branches until it enters the period-1 state. The event near ω = 0.935 is consistent with a sudden change that can drive the mesh into collision-induced instability. Above ω = 0.954, the response enters a sustained chaotic interval.
Increasing the error amplitude changes this path. At E = 29, the period-2 branch is further suppressed, and the response above ω = 0.951 becomes more irregular. At E = 30, the right-hand period-2 branch disappears and the system briefly returns to a single-period response. The bifurcation points therefore move along the frequency axis as the period-2 branches contract. The combined action of error and meshing frequency can suppress multi-period vibration over a limited range, but beyond a threshold it reinforces chaotic response. This non-monotonic effect is important for design: reducing transmission error is not simply a linear stabilizing action when the operating point lies near a branch boundary.
Error–stiffness and error–frequency domains
In the combined-error and stiffness-fluctuation plane, E ∈ [1.0, 2.6] and ε ∈ [0.01, 0.1], the map contains period-1, period-4, and chaotic domains. Period-1 regions appear as windows embedded in the chaotic background, with relatively clear boundaries. The low-error portion, approximately E < 1.22, is dominated by period-4 cells. Within this box, changing E produces more frequent state switching than changing ε; the error axis therefore carries the stronger sensitivity, even though the numerical span of ε is narrower. Boundary cells again identify critical parameter values where adjacent domains meet, and these cells trace the bifurcation set across the plane.
The error–frequency plane is more fragmented. Period-1, period-2, period-3, and chaotic cells coexist. Period-3 cells are scattered near the chaotic boundary, where unstable or weakly stable attractors make the local dynamics sensitive to small parameter changes. Period-1 and period-2 domains occupy most of the plane, while chaotic bands separate them. Period-1 regions appear on both sides of the chaotic zone. As E decreases from 2.5 to 1.2, the transition from period-1 to period-2 includes period-doubling behavior. Near ω ∈ [0.76, 0.8], period-3 and chaotic cells form a mixed band that can be treated as an unstable local region. Some boundary cells also indicate a direct transition from period-1 to chaos without an intervening period-doubling cascade. Maximum Lyapunov exponents along three fixed frequency cuts are consistent with the cell labels, supporting the numerical classification.
Damping ratio, high-period windows, and global occupancy
Damping ratio is a strong global control parameter. In the damping-ratio and error-amplitude plane, reducing ξ allows chaotic domains to dominate. Period-doubling windows appear in the sequence P2, P4, P8, and P16 as damping decreases. The high-period windows are narrow: only a small number of P16 cells lie near the front boundary of the chaotic region. Their limited area means that a coarse grid can miss them, while a finer grid sharply increases computational cost. The practical implication is that low damping raises the risk of unstable vibration, and the shape of the high-period windows is highly sensitive to the excitation parameters.
A direct check at E = 0.35 illustrates the damping effect. For ξ = 0.1, the Poincaré section is chaotic; for ξ = 0.3, it is period-4; for ξ = 0.4, it is period-2. The sequence is consistent with the domain map: lower damping moves the response toward chaos, while higher damping stabilizes lower-period motion.
Domain occupancy provides a simple global comparison. If ni is the number of cells classified as state i and N is the total number of cells, the occupancy share is S = ni/N. For period-4 cells, the damping–backlash plane gives 16.67%, while the damping–error plane gives 20.21%.
| Parameter plane | Period-4 occupancy |
|---|---|
| Damping ratio–backlash, ξ–b | 16.67% |
| Damping ratio–error amplitude, ξ–E | 20.21% |
These percentages are grid occupancy shares, not probabilities of failure. They depend on the chosen parameter box and cell resolution. Even so, the comparison indicates that the period-4 state occupies more of the sampled damping–error plane than of the damping–backlash plane, so the error amplitude is the stronger global influence on that state within these windows.
Design implications and limits
For power-combining spiral bevel gear systems, the useful output of a domain map is not a single stable point but a picture of where periodic windows are wide, where boundaries are close, and where chaotic bands interrupt a design path. A transmission-error limit chosen only from a one-dimensional sweep can land near a branch contraction or a fragmented period-3 band. A two-parameter map makes those risks visible: boundary cells mark critical combinations, period windows show where multi-period response can be avoided, and occupancy shares compare how strongly different excitations reshape the global response.
The maps also show why damping and transmission error should be treated together. Lower damping expands chaotic domains and narrows high-period windows; increased error can suppress a period-2 branch at one frequency and intensify irregular response at another. Stiffness fluctuation adds a third axis of sensitivity. For preliminary design, the practical strategy is to keep the operating point away from fragmented boundary bands, preserve margin in damping, and specify error limits with reference to the local window structure rather than to a single nominal bifurcation point.
The conclusions are bounded by the model. Rigid bodies, frictionless teeth, symmetric meshes, and sinusoidal stiffness fluctuation exclude effects such as tooth-surface friction, manufacturing mistuning, thermal distortion, and load sharing that varies with housing compliance. The domain classification is also grid-based: narrow period-16 windows and thin chaotic bands can be under-resolved, and a finer mesh may move a boundary by several cells. Experimental confirmation is therefore necessary before applying a specific map to a hardware design. Within those limits, the parameter-plane approach provides a compact way to connect multi-source excitation to the global periodic and chaotic structure of power-combining spiral bevel gears.