# Spiral Bevel Gears: nonlinear dynamic model based on accurate static stiffness evaluation

Moslem Molaie, Farhad Samani, Antonio Zippo, Francesco Pellicano

<figure><img src="image_1.png" /></figure>

## To cite this version:

Moslem Molaie, Farhad Samani, Antonio Zippo, Francesco Pellicano. Spiral Bevel Gears: nonlinear dynamic model based on accurate static stiffness evaluation. Journal of Sound and Vibration, 2023, 544, pp.117395. 10.1016/j.jsv.2022.117395. hal-03889150

HAL Id: hal-03889150

https://hal.science/hal-03889150

Submitted on 7 Dec 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Spiral Bevel Gear: nonlinear dynamic model based on accurate static stiffness evaluation

## Moslem Molaie

Department of Engineering "Enzo Ferrari"
University of Modena and Reggio Emilia, Modena, Italy
Email: moslem_molaie@unimore.it
ORCID: 0000-0002-2316-0803

## Farhad S. Samani

Associate Professor
Department of Mechanical Engineering
Shahid Bahonar University of Kerman, Kerman, Iran
Email: farhad.samani@uk.ac.ir
ORCID: 0000-0001-7887-1154

## Antonio Zippo

Department of Engineering "Enzo Ferrari"
Centre InterMech MoRe
University of Modena and Reggio Emilia, Modena, Italy
Email: antonio.zippo@unimore.it
ORCID: 0000-0001-6206-2619

## Francesco Pellicano

Corresponding author, Full Professor
Department of Engineering "Enzo Ferrari"
Centre InterMech MoRe
University of Modena and Reggio Emilia, Modena, Italy
Email: francesco.pellicano@unimore.it
ORCID: 0000-0003-2465-6584

**Abstract** In the present paper non-linear dynamics of a spiral bevel gear pair with backlash are investigated in order to clarify the internal excitations of major importance from the vibration point of view: manufacturing errors in the teeth profile, teeth spacing errors, and elastic deformation of the teeth. In some conditions, like in the case of backside contact, the destructive effect of internal excitations can be intensified leading to complex dynamics; for such reasons here backside contacts and reverse rotation are investigated in detail using a nonlinear time-varying model. The effect of damping is investigated as well. A one-DOF model is developed in order to study the dynamic behavior; the resulting a nonlinear differential equation with time-varying mesh stiffness is solved via numerical integration based on an adaptive step-size implicit Runge-Kutta scheme. The dynamic response of the system is analyzed through time histories, phase portraits, bifurcation diagrams, and Poincaré maps. Results show that for small backlash values, the possibility of backside contact increases. Meanwhile, by increasing the backlash value, the amplitude vibration of the gear rotation rises as well. By comparing the dynamic response of the system with different damping ratios, the results show that higher damping effectively reduces gear vibration resonance, although the probability of unsteady response still exists.

**Keywords:** nonlinear vibration, spiral bevel gear, gear mesh stiffness, forward and reverse motions, backside contact.

## List of Symbols

$a_j, b_j$ Fourier coefficients

$b$ Half of the gear backlash

$c$ Damping coefficient between the mesh gear teeth of the pairs

$C_{eq}$ Equivalent damping coefficient

$E$ Module of elasticity

$e(t)$ Time-varying circumferential no-load transmission error
$F_n$ The normal dynamic load for the driven gear

$F_z$ The Z-component of the normal dynamic load for the driven gear

$I_p^x, I_g^y$ Rotary inertia of pinion and gear

$m_{eq}$ Equivalent mass

$N_1$ Teeth number of the pinion

$n$ Gear ratio of the gear pair

$N_p$ Number of samples for mesh stiffness computation

$k_0$ Average value of torsional mesh stiffness of the gear pair

$K_{eq}$ Equivalent mesh stiffness of the gear pair

$K_m$ The torsional mesh stiffness of the gear pair

$r_p, r_g$ Base radii of the pinion and the gear

$r_{pitch-p}$ pitch radius of the pinion at the mid-section

$r_{pitch-g}$ pitch radius of the gear at the mid-section

$S$ Number of harmonics

$T_{eq}$ Equivalent applied torque on the driven gear

$T_l$ Constant driver torque

$T_m$ Constant breaking torque

$w$ Face width

$\alpha$ Normal pressure angle

$\beta$ The spiral angle

$\gamma_s$ Input shaft speed

$\zeta$ Damping ratio

$\theta_1$ Driver angular displacement

$\theta_2$ Driven angular displacement

$\theta_b$ Angular backlash

$\lambda$ Linear dynamic transmission error along the line of action

$\lambda_\theta$ Angular dynamic transmission error

$\nu$ Poisson ratio

$\tau$ Non dimensional parameter for time

$\omega_m$ Fundamental mesh frequency

$\omega_n$ Natural mesh frequency
# 1. Introduction

Bevel gears embrace a wide range of applications in different engineering fields: in aerospace, bevel gears are present in helicopter gearboxes, in terrestrial vehicles, they are often used in the final stage of powertrain; in heavy industry whenever it is required to transmit a high load between non-parallel shafts. There are many types of bevel gears, the most significant of which are the spiral bevel gears (SBGs). The main advantage of SBGs is the high-power transmission density compared with other types of bevel gear pairs, for instance, straight or Zerol bevel gears.

The vibration and dynamic behavior of gear pairs have been the subject of intense research in the past [1-5] since their influence is pervasive on noise and stress distribution, and consequently, on the lifetime, the efficiency and the capability of transmitting power. Yavuz et al. [6] proposed a nonlinear dynamic model for the SBG pair of a gear train that took into account the stiffness of the shafts and bearings. The Tredgold approach was used to optimize straight bevel gear models by Motahar et al. [3]. Yassine et al. [7] studied the dynamic responses of a two-stage straight bevel system with eccentricity, profile error, and a fractured tooth. The Newmark approach was used to determine the dynamic response. Masoumi et al. [8] investigated complex dynamics and chaos in high speed balanced planetary gears, proving that the nonlinearity and the stiffness variation can induce vibration in inherently balanced systems.

As the main source of vibration in most of gear trains is the static transmission error (STE), this topic attracted the attention of an increasing number of scientists and engineers, with a consequent large production of scientific and technical papers. The total deflection of gears during meshing is given by the summation of different deformations: teeth bending, shear and surface indentation. Kickbush et al. [9] proposed two finite element models (two-dimensional and three-dimensional) to approximate the mesh stiffness (MS). Tang et al. [10] investigated how two distinct STEs affected SBG's dynamic responsiveness; a predesigned parabolic function and a sine function were the two STEs evaluated. Some other research investigated how faults impact gear pairs' dynamic behavior [11, 12]. Peng et al. [13] proposed a new method for estimating the loaded transmission error (LTE) taking into account the bearing supports' influence. Using a finite element analysis technique, Wang et al. [14] examined the time-varying mesh stiffness of a gear pair with cracked teeth. In Ref. [15] a new approach for creating very accurate nonlinear finite element models was proposed and the effect of misalignments was investigated in the case of helical gears.

The main goal of the present study is to investigate the static transmission error of a SBG and its nonlinear dynamic behavior by considering an example of SBG which is frequently used in transmission power system of helicopters. The accurate static analysis of LTE is essential to calculate the system's mesh stiffness and plays a significant role in understanding the gear stress status under working conditions. As the example of SBG is used for a helicopter gearbox, the method to manufacture the gear pair, is face milling process which provides the possibility of grinding. The static analysis is validated
with both experimental and finite element method (FEM) results by extracting the root stress and tooth force. Due to the mismatched teeth surfaces between gear and pinion, loaded tooth contact analysis (LTCA) and unloaded tooth contact analysis (UTCA) must be conducted. The main result of the static simulations consists of the evaluation of the mesh stiffness, through LTCA, for both forward and reverse motion. The dynamic model is used for investigating the nonlinear dynamic scenario, the analysis is carried out through a 1-DOF model that includes: the effect of backlash, i.e. non-smooth equation; variable stiffness, i.e. time-varying coefficients equation and consequent parametric excitation; effect of stiffness variation and phasing in the case of backside contact (reverse motion). The scenario is studied through the development of bifurcation diagrams, amplitude-frequency diagrams, Poincaré maps, and phase portraits. Besides, the effect of some important parameters, including the damping ratio and the backlash, on the frequency-response curves is investigated.

## 2. Loaded and unloaded tooth contact analysis

Bevel gears are classified into different categories based on the employed manufacturing process. The type of tooth trace determines one of the classifications: straight bevel gears, hypoid bevel gears, Zerol bevel gears and spiral bevel gears. Spiral bevel gears, when manufactured in a metal-cutting process, can either be produced in single indexing (face-milled process) or continuous indexing (face-hobbed process) operations which dictate the shape of the tooth trace [16-19], see Fig. 1.

<figure><img src="image_3.png" /><figcaption>Fig. 1. Single indexing (a) and continuous indexing (b) methods.</figcaption></figure>

These are two ubiquitous processes which are used to manufacture spiral bevel: face milling and face hobbing methods. In the face milling approach, one side of a tooth is generated by the cutter. The cutter is then withdrawn from the gear blank; the workpiece is indexed by one tooth at time; the cutter returns to this new position, and the process is repeated [20]. The main advantage of face milling process is that it can be used for cutting and grinding. There are two types of the face milling processes, Formate and Generated [18]. The ‘Formate’ face-milling process is a method with no generating motion (faster) and is used in automotive ring gears. While the ‘Generated’ face-milling process is a method with
generating motion and is used for automotive pinions as well as for both the pinion and ring gear in aerospace applications.

In contrast with the face milling method, the face hobbing process is more efficient as it requires fewer machine setups. However, the noteworthy negative point of this process is that there is no possibility of grinding the gear pair generated by this method. Due to this point, face milling is mainly used in aerospace applications, where the gears need to be ground, whereas face hobbing is a predominant method in automotive applications where machining time is more of a concern than surface finish [17]. To sum up, in face milling process, both gear and pinion design data must be given as well as the kinematics of gear cutting machine and its settings, all these parameters are necessary to create an SBG-pairs.

Due to the significant development in gear software packages during the last decade, there is the opportunity to simulate a gear pair based on analytical or FEM methods before creating gears. Moreover, the identification of defects within a transmission system, reducing equipment failures and unscheduled downtime is increasingly demanded. The major accomplishment of simulating models in software before manufacturing is to understand its efficiency and behavior under real situation and the ability to optimize and improve the gearbox design and efficiency.

The vibration analysis has proven to be one of the most effective tools for identifying mechanical faults within gear systems and choose the best model for their system. One of initial steps of conducting vibration analysis is LTCA to calculate the mesh stiffness of the system. To determine the mesh stiffness of a gear pair, an LTCA must be performed. The LTCA could be carried out using a variety of approaches: FEM, experiments, and analytical solutions (with certain simplifications). By comparing these methods, one can conclude that: i) the experimental method requires a well-equipped laboratory, it is cost and time demanding, moreover, measurements on rotating parts are limited to non-contact surfaces; ii) analytical solutions have been generally developed for spur or helical gear pairs, however, they are not perfectly suitable or precise for gear pairs with complicated geometries, such as a spiral gear pair; iii) nonlinear FEM (including contact analysis) has become increasingly popular in recent years [3, 21, 22].

SBGs are manufactured with an intentional tooth profile error, i.e. mismatch [23-25] between two teeth surfaces of gear and pinion. This mismatch is done so that the gear mesh can perform in an acceptable way even in the presence of misalignment, manufacturing errors, or high torque level that cause the gears to deflect to locations different than what was defined in their design. Although the conjugate motion provides an ideal transmission, i.e. constant transmission ratio, this is not practical in SBG; consequently, SBGs generally operate with a certain amount of geometric transmission error [24, 26]. Therefore, to obtain the correct mesh stiffness, both LTCA and UTCA must be carried out.

The process, considered in this paper, is the face-milled method. Such kind of SBGs are mainly used in aerospace power transmissions (i.e., helicopter main/tail rotor transmissions) to transmit power
between horizontal gas turbine engines and the vertical rotor shaft [27]. The gear and pinion considered in the present work were experimentally analyzed through the NASA Lewis Spiral Bevel Gear Test Facility [28]. An image of the spiral bevel gear mesh is shown in Fig. 2, the design data for the pinion and gear and the machine settings are listed in Appendix A.

a)

<figure><img src="image_4.png" /><figcaption>Fig. 2. Spiral bevel gear and pinion, a) the test specimen of the NASA Spiral Bevel Gear Rig [28], b) finite element model (Ansol Calyx)</figcaption></figure>

b)

<figure><img src="image_5.png" /></figure>

The results, which are available for the pinion root stresses, are interesting from validation point of view. The test facility used in their testing is the Spiral Bevel Gear Test Facility at NASA Glenn Research Center [29]. The pinion of the spiral bevel gear mesh was instrumented with strain gauges. A picture of the pinion used in the experiments is shown in Fig. 3. Three successive teeth of the pinion were instrumented with the five strain gauges. A total of three strain gauges were installed on the middle tooth (see Fig. 4). The strain gauges used in their study had an active gauge length of 0.38 mm (0.015 in) [29]. The contact performance of the left-hand flank pairs is analyzed, which corresponds to the concave side of the pinion tooth and the convex side of the gear tooth. The kinematics of gear cutting are incorporated in the computer programs, ANSOL-Calyx (Calyx-HypoidFaceMilled product from ANSOL company - Version 2021), to generate the spiral bevel gear model (see Appendix A).

a)

<figure><img src="image_6.png" /><figcaption>Fig. 3. Strain gauges on pinion for experiments, a) the location of strain gauges on the pinion tooth root, b) a close-up of strain gauge location [27, 30].</figcaption></figure>

b)

<figure><img src="image_7.png" /></figure>
<figure><img src="image_8.png" /><figcaption>Fig. 4. Strain gauge location on the successive tooth along the pinion fillet</figcaption></figure>

Fig. 5 shows the comparison between experimental results [29] and the present FEM analysis: the red line represents the root stress for the pinion tooth derived from experimental data (strain) and black one represents root stress extracted from simulation. The comparison shows an excellent agreement.

<figure><img src="image_9.png" /><figcaption>Fig. 5. Maximum principal stress on the root of the tooth respect to pinion for a) the heel root, b) the middle root, and c) the toe root - torque value is 7840 lb.in on the gear, — FEM (Calyx) simulation, — Experimental results [29].</figcaption></figure>

A second comparison is conducted in order to have a further proof of the accuracy of the present model. In this case the reference consists of a numerical model developed by Bibel et al. [31]: they carried out a nonlinear FEM analysis by MSC Marc in order to simulate the gear pair described in Appendix A. To simulate the model, three pinion teeth and four gear teeth were considered. The pinion was fixed (Fig. 6-a) and the gear was allowed to rotate about its axis only [26, 31].
<figure><img src="image_10.png" /><figcaption>Fig. 6. FEM model a) NASA - MSC Marc [31], b) Calyx.</figcaption></figure>

The first parameter reported by Bibel et al. [31] is the tooth force on the pinion tooth. Fig. 7 represents the results from Ref. [26] (blue and red spots), beside to the those obtained by the present simulation (Calyx).

<figure><img src="image_11.png" /><figcaption>Fig. 7. Tooth force of pinion: torque value is 9508 lb.in on the gear, ‘—’ present FEM simulation, ● nonlinear FEM (MSC Marc) [31].</figcaption></figure>

Moreover, Fig. 8 shows the maximum principal stress on the root of the tooth with respect to the heel, middle, and toe roots of the pinion. The solid line shows the simulations and the filled circles are FEM results reported by Bibel et al. [31]; at that time, using FEM to simulate an SBG model was difficult, therefore it is conducted by some simplifications. For instance, to connect flexible members together, gap elements were required. This leads to many difficulties as the gap elements needed to be oriented normal to the surfaces as closely as possible [26]. Also, due to some limitation on the computer facilities, just a few teeth of pinion and gear were modeled. These disparities lead to some differences between the present results and FEM results of Ref. [31], however the comparison is quite satisfactory.
<figure><img src="image_12.png" /><figcaption>Fig. 8. Maximum principal stress on the root of the pinion tooth: heel root (a), middle root (b), and toe root (c). Torque is 9508 lb.in on the gear. ‘—’ present FEM simulation, ● FEM (MSC Marc) [31])</figcaption></figure>

## 2.1 FEM grid sensitivity analysis

In this section the finite element grid refinement is analyzed in order to assess the quality of the simulations and to determine the appropriate finite element size for static simulation. The following points are considered: performing selected analyses on a variety of mesh sizes, noting areas with significant deformations or stresses; improving the mesh in these areas; collecting data from each mesh analysis (Fig. 9), including the result, the number of divisions in the model, and the time required for computing.

From a mesh sensitivity standpoint, the primary issue is determining if a mesh is sufficiently refined, since reducing element size results in enormous computational times. Additionally, a balance between computing time and accuracy should be achieved. From the chart in Fig. 9 and Table 1 one can check that after the 3rd step, any considerable improvement in accuracy would require a tremendous increase in computing time.

<table><caption>Table 1. Grid sensitivity analysis</caption><thead><tr><th>Case</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th></tr></thead><tbody><tr><td>Average deflection of gear [μrad]</td><td>75.86</td><td>74.77</td><td>73.60</td><td>73.49</td><td>73.52</td><td>73.49</td></tr><tr><td>Number of sections along the tooth width</td><td>3</td><td>8</td><td>13</td><td>14</td><td>15</td><td>16</td></tr><tr><td>Number of elements on the tooth section</td><td>36</td><td>74</td><td>176</td><td>176</td><td>176</td><td>176</td></tr><tr><td>Computation time [minute]</td><td>260</td><td>286</td><td>398</td><td>420</td><td>437</td><td>454</td></tr></tbody></table>
<figure><img src="image_13.png" /><figcaption>Fig. 9. Grid sensitivity analysis - torque value is 4754 lb.in on the gear.</figcaption></figure>

## 2.2 FEM results of loaded and unloaded tooth contact analysis

A valuable parameter is the maximum principal stress on the root of the pinion. Fig. 10 shows this parameter for both directions of rotation, *forward* (solid line) and *reverse* (dash line) motions. The results are brought for three different teeth of pinion; for the forward motion, tooth no. 1, which is close to lose the contact (red line), tooth no. 2, which is passing the contact (black line), and finally the tooth no. 3 (blue line), which will come to contact after tooth no. 2; for the reverse motion, the sequence is the reverse. The torque applied to the gear is 4754 lb.in; the consequent maximum value of the principal stress in the forward motion is equal to 79950 psi while this value for the reverse motion is lower, 77206 psi.

<figure><img src="image_14.png" /><figcaption>Fig. 10. Maximum principal stress on the root of pinion - torque on the gear 4754 lb.in; forward motion (solid line), reverse motion (dash line); tooth no. 1 (red line), tooth no. 2 (black line), tooth no. 3 (blue line)</figcaption></figure>
Another important parameter is the force on the tooth surface; Fig. 11 shows such parameter for the pinion, reverse motion (Fig. 11-a) and forward motion (Fig. 11-b). From these results one we find that the maximum number of teeth in contact is 2, i.e. the contact ratio is between 1 and 2.

<figure><img src="image_15.png" /><figcaption>Fig. 11. Force on the tooth surface of the pinion - torque value 4754 lb.in on the gear; a) forward motion, b) reverse motion; tooth No. 1 (red line), tooth No. 2 (black line), and tooth No. 3 (blue line)</figcaption></figure>

The contact pressure on the tooth surface could be precious results to explore the critical zone on the tooth contact surface. As the gears are spiral bevel gear, the stress distribution on the gear is not uniform along the tooth width, Fig. 12. The contact zone during reverse motion is close to the heel edge while during forward motion it comes to the toe edge. The critical stress zones for both motions are approximately at the middle of tooth surface.

<figure><img src="image_16.png" /><figcaption>Fig. 12. Contact pressure on the tooth surface - torque value is 4754 lb.in on the gear;</figcaption></figure>
As mentioned before, here the main purpose of static simulation is obtaining the mesh stiffness (MS) needed to investigate the dynamic behavior of the gear pair. To calculate the MS, it is required to know the tooth deflection by applying a certain amount of torque. Therefore, it is required to conduct the LTCA and UTCA due to the existence of some manufacturing errors such as the mismatched between two surfaces of pinion and gear teeth. Fig. 13 shows the rotational transmission error when the torque applied on the gear; besides, Fig. 14 presents the UTCA results while a very low torque (0.01lb.in) is applied on the gear, i.e., the geometric transmission error of the gear pair.

Eventually, Fig. 15 shows the difference between the transmission error obtained from LTCA, normal torque, and the one given by UTCA, negligible torque; this part of the transmission error is due to the elastic deformation only, therefore, it is the deformation to be used for calculating the meshing stiffness.

<figure><img src="image_17.png" /><figcaption>Fig. 13. Transmission error (gear rotation) after LTCA - torque 4754 lb·in; a) forward motion, b) reverse motion; — Fitted curve, ● FEM results.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 14. Rigid body transmission error (gear rotation) after UTCA - torque value is 0.01 lb·in;
a) forward motion, b) reverse motion; — Fitted curve, ● FEM results.
a) forward motion
b) reverse motion</figcaption></figure>
<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 15. Elastic transmission error after LTCA and UTCA; a) forward motion, b) reverse motion; — Fitted curve, ● FEM results.</figcaption></figure>

## 3 Dynamic physical model

Consider the spiral bevel gear pair of Fig. 16, the translational degrees of freedom for both, driver and driven gears are constrained in all directions as well as the rotations; the gears can only rotate around their axes.

<figure><img src="image_21.png" /><figcaption>Fig. 16. The dynamic model of a gear system with rotational degrees of freedom.</figcaption></figure>

The dynamic equations of motion of this system (Fig. 16) are given by [2, 3, 32]:

$$I_g^y \ddot{\theta}_g = -r_g F_z - T_l \quad (1)$$

$$I_p^x \ddot{\phi}_p = r_p F_z + T_m \quad (2)$$

Due to mounting and manufacturing error (Fig. 17) or teeth profile modifications, a local gap between mating teeth can appear, it is called geometric transmission error, $e(t)$; it can be calculated through the UTCA. The linear dynamic transmission error (DTE) along the line of action is defined as:

$$\lambda = (r_{pitch-p} \phi_p - r_{pitch-g} \theta_g)a + e(t) \quad (3)$$

where:
$$a = \cos \alpha \cos \beta$$
$$e(t) = \sum_{j=1}^{s} e_j \cos(j\omega_m t) \qquad (4)$$

$\beta$ is the spiral angle and $\alpha$ is the normal pressure angle,

<figure><img src="image_22.png" /><figcaption>Fig. 17. Geometric transmission error: lack of material along the line of action</figcaption></figure>

The dynamic load of pinion and its component along the line of action $F_n$ and along the z-axis $F_z$ can be calculated as:

$$F_n = K_m(t)f(\lambda - e) + C_m(\dot{\lambda} - \dot{e}) \qquad (5)$$

$$F_z = -F_n a \qquad (6)$$

where:

$$f(\lambda - e) = \begin{cases} \lambda + b - e, & \lambda - e > b \\ 0, & -b \le \lambda - e \le b \\ \lambda - b - e, & \lambda - e < -b \end{cases} \qquad (7)$$

$f(\lambda - e)$ is the function of the linear displacement (7), Fig. 18. This function with the multiplication of stiffness returns the restoring force function [33].

Whenever $\lambda - e$ is between $-b$ and $+b$, the contact loss happens (single-sided impact) [34, 35]. In the case of forward motion, for $\lambda - e > b$, the mesh is expected to be in the forward contact (desired situation), while if $\lambda - e < -b$, undesired backside contact happens (double-sided impact); see Ref. [35, 36].

<figure><img src="image_23.png" /><figcaption>Fig. 18. Equivalent gear model and backlash function.</figcaption></figure>

The meshing stiffness of the gear pair is a time-varying periodic function, with the fundamental mesh frequency $\omega_m = \frac{2\pi}{60} N_1 \gamma_s$; the equivalent meshing stiffness Fourier series is given by:
$$ \left\{ \begin{array}{l} \omega_m = \frac{2\pi}{60} N_1 \gamma_s \\ s = (N_p - 1)/2 \end{array} \right. \quad \Rightarrow \quad K_m(t) = k_0 + \sum_{j=1}^{S} a_j \cos(j\omega_m t) + \sum_{j=1}^{S} b_j \sin(j\omega_m t) \quad (8) $$

From the pure elastic deformation obtained from the LTCA and UTCA, see Fig. 15, the meshing stiffness for both directions of motions is extracted (see Fig. 19); moreover, Fig. 20 represents the geometric transmission error.

<figure><img src="image_24.png" /><figcaption>Fig. 19. Mesh stiffness diagram:</figcaption></figure>

a) forward motion, b) reverse motion; — Fitted curve by Fourier series, ● FEM results.

<figure><img src="image_25.png" /><figcaption>Fig. 20. Geometric transmission error on the line of action:</figcaption></figure>

a) forward motion, b) reverse motion; — Fitted curve by Fourier series, ● FEM results.

<figure><img src="image_26.png" /><figcaption>Fig. 21. Dynamic model of a reduction gear system with one DOF</figcaption></figure>

Therefore, by combining Eqs. (1) and (2), Eq. (9) is written in terms of DTE:
$$m_{eq}\ddot{\lambda} + C_m(\dot{\lambda} - \dot{e}) + K_m(t)f(\lambda - e) = T_{eq} \quad (9)$$

where,

$$m_{eq} = \left( \frac{(r_p)^2}{l_p^x} + \frac{(r_g)^2}{l_g^y} \right)^{-1}, \quad T_{eq} = \left( \frac{T_m}{r_p} \right) \quad (10)$$

In order to normalize the governing equation, new parameters are introduced as follows:

$$\zeta = \frac{C_m}{2m_{eq}\omega_n}, \quad \bar{T}_m = \frac{1}{b m_{eq} \omega_n^2} \left( \frac{T_m}{r_p} \right) \\ \tau = \omega_n t \quad \lambda' = \frac{d\lambda}{d\tau} \quad \bar{\lambda} = \frac{\lambda}{b} \quad \omega_n = \sqrt{\frac{k_0}{m_{eq}}} \quad \bar{e} = \frac{e}{b} \quad (11)$$

Consequently, Eq. (9) can be rewritten as follows:

$$\bar{\lambda}'' + \bar{K}_m(\tau)f(\bar{\lambda} - \bar{e}) + 2\zeta(\bar{\lambda}' - \bar{e}') = \bar{T}_m \quad (12)$$

$$f(\bar{\lambda} - \bar{e}) = \begin{cases} \bar{\lambda} - \bar{e} + 1, & \bar{\lambda} - \bar{e} > 1 \\ 0, & -1 \le \bar{\lambda} - \bar{e} \le 1 \\ \bar{\lambda} - \bar{e} - 1, & \bar{\lambda} - \bar{e} < -1 \end{cases} \quad (13)$$

$$\bar{K}_m(\tau) = 1 + \sum_{j=1}^{s} \frac{a_j}{m_{eq}\omega_n^2} \cos(j\omega_m \tau) + \sum_{j=1}^{s} \frac{b_j}{m_{eq}\omega_n^2} \sin(j\omega_m \tau) \quad (14)$$

Eq. (12) is a nonlinear differential equation with time-varying parameters. This equation is solved numerically based on an implicit fifth-order Runge-Kutta scheme (RADAU) coded in FORTRAN language, the algorithm is extremely stable and accurate, see e.g. Refs. [3, 21, 22, 33, 36].

## 3.1 Numerical results

As explained before, Eq.(12) describes a non-linear differential equation with time-varying coefficient, which parametrically excites system in addition to the external constant torque (power transmitted). The time responses calculated by direct simulation are used for building amplitude-frequency diagrams and bifurcation diagrams of Poincaré maps, obtained by varying the excitation frequency: the pinion rotation speed. For each simulation the transient is separated from the “steady-state” for each frequency, i.e. the first 1000 periods are deemed transient and beyond that, the influence of frequency variation is assumed to be expired, therefore, the first 1000 periods of simulation are never recorded and analyzed. The Poincaré map diagrams and phase portraits are derived from at the last 200 periods of time responses. Each Poincaré map is obtained by sampling the time histories with the same period of the meshing frequency.

For all simulations the torque applied to the pinion is 4754 lb·in (537.12 N.m) which correspond to the 50% of the nominal torque. Fig. 22 shows the bifurcation diagram of the Poincaré maps, for each frequency a one-dimensional Poincaré map section is represented, i.e. for each frequency 200 points are plotted. By means of this diagram the dynamic behavior of the system is represented, in Fig. 22 the gear pairs rotate in the forward direction. It should be noted that black dots indicate the responses during
the backward simulation, i.e. the rotation speed starts from highest to the lowest value (from high to low excitation frequency); red dots indicate the forward simulation (from low to high speed).

Let us consider the forward simulation, the system is excited from low frequency and the frequency is increasing (i.e., sweeping up the frequency), the response will grow until the point that any further increase in excitation frequency will cause a spontaneous jump in the amplitude of DTE (see Table 2 – Case D; starting point with $\frac{\omega_m}{\omega_n} = 0.3266$), after the jump a further increasing of the frequency causes a reduction of amplitude see Fig. 22b; this peak is a resonance caused by the third super-harmonic of the meshing stiffness. From the bifurcation diagram, Fig. 22a, one can see that up to $\frac{\omega_m}{\omega_n} = 0.46$ the system response is periodic (single line of the diagram). In Table 2 the list of different jumps is reported. The resonance close to $\frac{\omega_m}{\omega_n} = 0.5$ is due to the second super-harmonic of the meshing stiffness and the main resonance close to $\frac{\omega_m}{\omega_n} = 1$ is due to the fundamental harmonic; all resonances have a strong sub-harmonic character. Close to $\frac{\omega_m}{\omega_n} = 2$ an instability due to the principal parametric resonance is present, the parametric resonance branch shows a softening character. Such softening behaviors are due to the loosening of contact that take place when the inertia forces exceed the static contact force induced by the applied torque.

The sweep-down simulation completes the analysis allowing to follow all stable branches, in particular the softening branches are almost completely followed. It is interesting to note that a small portion of a hardening branch is caught, this happens for very high vibration amplitudes when the backside contact takes place; the hardening branch starts from the ending point of case A with $\frac{\omega_m}{\omega_n} = 0.5072$ (see Table 2; case A) to starting point of case B with $\frac{\omega_m}{\omega_n} = 0.4716$ (see Table 2; case B).

From a theoretical point of view, different solutions of the governing equation could exist; however, some of them are practically not feasible and cannot be physically achieved. The transition appearing in proximity of the jumps are not recorded. Besides, the system experiences unsteady phenomena, which can be chaotic or quasiperiodic, when it approaches super-harmonic resonances ($\frac{\omega_m}{\omega_n} \approx 0.46, 0.51$, respectively for backward and forward simulations).
<figure><img src="image_27.png" /><figcaption>Fig. 22. Dynamic response for forward motion, torque 4754 lb·in, a) Bifurcation diagram, b) amplitude-frequency diagram. Forward simulation (red dot), Backward simulation (black dot).</figcaption></figure>
The backlash function is the main source of nonlinearity; the system could experience three different states: drive side contact, tooth separation, and backside contact (coast-side contact). As it was explained, when the amplitude of $\lambda$ (DTE) is between -1 and 1, the tooth separation takes place (see cases A, B, C and D in Table 2). Teeth separation occurs in all the three frequency ranges: fundamental, super-harmonic, and sub-harmonic, for the backward and forward simulations. As explained, as soon as the nondimensional amplitude of the response gets lower than -1, the undesired backside contact occurs (see cases A' in Table 2) near to the super-harmonic frequency region for the backward simulation, i.e., reducing speed.

<table><caption>Table 2. The properties of starting and ending points of separation and backside contact for the forward motion</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">Backward simulation</th><th colspan="2">Forward simulation</th></tr><tr><th>starting point</th><th>ending point</th><th>starting point</th><th>ending point</th></tr></thead><tbody><tr><td rowspan="12">Teeth separation</td><td rowspan="3">A</td><td>ω<sub>m</sub>/ω<sub>n</sub></td><td>2.1043</td><td>0.5072</td><td>1.8911</td><td>2.1062</td></tr><tr><td>min (λ)</td><td>0.9967</td><td>-0.9887</td><td>0.9017</td><td>0.9996</td></tr><tr><td>max (λ)</td><td>1.5879</td><td>2.0266</td><td>1.6379</td><td>1.5871</td></tr><tr><td rowspan="3">B</td><td>ω<sub>m</sub>/ω<sub>n</sub></td><td>0.4716</td><td>0.3516</td><td>0.8490</td><td>1.1516</td></tr><tr><td>min (λ)</td><td>0.8851</td><td>0.3153</td><td>0.7031</td><td>0.9700</td></tr><tr><td>max (λ)</td><td>1.6877</td><td>1.8971</td><td>1.6216</td><td>1.5149</td></tr><tr><td rowspan="3">C</td><td>ω<sub>m</sub>/ω<sub>n</sub></td><td>3.362</td><td>0.2949</td><td>0.4851</td><td>0.5120</td></tr><tr><td>min (λ)</td><td>0.9994</td><td>0.7605</td><td>0.9384</td><td>0.9945</td></tr><tr><td>max (λ)</td><td>1.5728</td><td>1.7164</td><td>1.6438</td><td>1.5834</td></tr><tr><td rowspan="3">D</td><td>ω<sub>m</sub>/ω<sub>n</sub></td><td></td><td></td><td>0.3266</td><td>0.3362</td></tr><tr><td>min (λ)</td><td>No separation</td><td></td><td>0.9487</td><td>0.9994</td></tr><tr><td>max (λ)</td><td></td><td></td><td>1.6081</td><td>1.5728</td></tr><tr><td rowspan="3">Backside contact</td><td rowspan="3">A'</td><td>ω<sub>m</sub>/ω<sub>n</sub></td><td>0.5062</td><td>0.4726</td><td></td><td></td></tr><tr><td>min (λ)</td><td>-1.3628</td><td>-1.0129</td><td></td><td></td></tr><tr><td>max (λ)</td><td>2.4443</td><td>2.3311</td><td></td><td></td></tr></tbody></table>

A further analysis is now carried out in the case of reverse direction, the system is investigated by means of Bifurcation and RMS diagrams (see Fig. 23). Compared to forward motion, here, just tooth separation occurred (see Table 3). In Fig. 23-a, one can see two diverse trends: periodic (single line), 2T-subharmonic (2 lines) during both forward and backward simulations. No hardening branches (no backside contact) nor complex dynamics are found in this case.
<figure><img src="image_28.png" /><figcaption>b) RMS Diagram</figcaption></figure>

<figure><img src="image_29.png" /></figure>

Fig. 23. Dynamic response for reverse motion, torque 4754 Ib·in, a) Bifurcation diagram, b) amplitude-frequency diagram. Forward simulation (red dot), Backward simulation (black dot).

During reverse motion, teeth separation is the only phenomena that occurs in all the three frequency ranges: fundamental, super-harmonic, and sub-harmonic, for the backward and forward simulations (see Table 3).

<table><caption>Table 3. The properties of starting and ending points of separation and backside contact for the reverse motion</caption><thead><tr><th rowspan="2"></th><th rowspan="2"></th><th colspan="2">Backward simulation</th><th colspan="2">Forward simulation</th></tr><tr><th>starting point</th><th>ending point</th><th>starting point</th><th>ending point</th></tr></thead><tbody><tr><td rowspan="9">Teeth separation</td><td rowspan="3">A</td><td>ω<sub>m</sub>/ω<sub>n</sub></td><td>2.0418</td><td>1.3024</td><td>1.9573</td><td>2.0438</td></tr><tr><td>min (λ)</td><td>0.9993</td><td>-0.9943</td><td>0.8827</td><td>0.9998</td></tr><tr><td>max (λ)</td><td>2.1285</td><td>2.9623</td><td>2.2207</td><td>2.1279</td></tr><tr><td rowspan="3">B</td><td>ω<sub>m</sub>/ω<sub>n</sub></td><td>1.0795</td><td>0.6551</td><td>0.9172</td><td>1.0795</td></tr><tr><td>min (λ)</td><td>0.9940</td><td>-0.9880</td><td>0.6418</td><td>0.9940</td></tr><tr><td>max (λ)</td><td>2.0507</td><td>2.8879</td><td>2.2602</td><td>2.0507</td></tr><tr><td rowspan="3">C</td><td>ω<sub>m</sub>/ω<sub>n</sub></td><td>0.5072</td><td>0.4438</td><td>0.4927</td><td>0.5072</td></tr><tr><td>min (λ)</td><td>0.9987</td><td>0.5340</td><td>0.8859</td><td>0.9987</td></tr><tr><td>max (λ)</td><td>2.0575</td><td>2.3478</td><td>2.1534</td><td>2.0575</td></tr></tbody></table>

The nonlinear dynamics of the most interesting regimes are shown in Fig. 24; time histories, phase portraits, and Poincaré maps are illustrated for four different frequency ranges: Case I) $\frac{\omega_m}{\omega_n} = 0.459$, the response is non periodic and the Poincaré map presents a complex shape; Case II) $\frac{\omega_m}{\omega_n} = 1.684$ the Poincaré maps show two points, the response is 2T subharmonic due to the nonlinearity induced by the teeth separation; Case III) $\frac{\omega_m}{\omega_n} = 1.098$, forward, shows a ¼ subharmonic response; Case IV) $\frac{\omega_m}{\omega_n} = 0.512$, forward, shows another unsteady response having probably a chaotic character as the Poincaré maps appears fractal with a coherent structure.
<figure><img src="image_30.png" /></figure>

<table><caption>Fig. 24. Nonlinear vibration behavior in different frequency ratios for backward simulation and forward motion</caption><thead><tr><th>State: <br>(ω<sub>m</sub>/ω<sub>n</sub>)</th><th>Phase portrait<br>Horizontal axis: λ(τ)<br>Vertical axis: λ̇(τ)</th><th>Poincaré map<br>Horizontal axis: λ(τ<sub>0</sub> + iT)<br>Vertical axis: λ̇(τ<sub>0</sub> + iT)</th><th>Time Response<br>Horizontal axis: τ = ω<sub>n</sub>t<br>Vertical axis: λ(τ)</th></tr></thead><tbody><tr><td>I:<br>(0.459)<br>Backward<br>simulation</td><td></td><td></td><td></td></tr><tr><td>II:<br>(1.684)<br>Backward<br>simulation</td><td></td><td></td><td></td></tr><tr><td>III:<br>(1.098)<br>Forward<br>simulation</td><td></td><td></td><td></td></tr><tr><td>IV:<br>(0.512)<br>Forward<br>simulation</td><td></td><td></td><td></td></tr></tbody></table>
## 3.2. Backlash effects on the dynamic behavior

Backlash is required to avoid binding of the mating gears, which causes heat production, noise, anomalous wear, overload, and/or drive failure. In addition to preventing binding, dimensional limitations required for cost-effective production induce backlash in gear systems. Backlash is created during the gear manufacturing process by making each gear tooth one-half as thin as the backlash size required for the application. When two of such gears run together (i.e., mate), their allowances add up to the whole backlash [37]. As backlash increases due to the tooth wear, here the effect of backlash is investigated by decreasing the thickness of teeth. A static FEM simulation is carried out to calculate the mesh stiffness for each case with different backlash. Fig. 25 represents the mesh stiffness in the form of Fourier transform coefficients for forward motion. One can see that, by increasing the backlash, the mean value of mesh stiffness decreased. The values of the first harmonic for all cases are about 5 time higher than the second harmonic's values and values are insignificant from harmonic number 6th.

<figure><img src="image_31.png" /><figcaption>Fig. 25. mesh stiffness for cases with different backlash in form of Fourier transform coefficients: forward motion, ● b = 0.01 in, ● b = 0.02 in, ● b = 0.03 in, ● b = 0.04 in</figcaption></figure>

To illustrate the effect of backlash on the dynamic response of the system, the governing equation is solved by considering diverse values for the backlash (see Fig. 26). By comparing these results, it is evident that backlash affects the dynamic responses of the system from two different points of view: maximum amplitude of the response and behavior of the system when the frequency is close to the primary resonance. Increasing the backlash leads to a decrease in the probability of backside contact.
<figure><img src="image_32.png" /><figcaption>Fig. 26. Effect of backlash on dynamic response during the forward motion
• b = 0.01 in, • b = 0.02 in, • b = 0.03 in, • b = 0.04 in</figcaption></figure>

3.3 Damping ratio effects on the dynamic behavior

In this section, the role of damping ratio on the nonlinear dynamics is investigated. Different values
of the damping ratio are considered (ζ = 0.01, 0.03, and 0.05 [35]), all the other parameters are kept
fixed, the frequency-response curves are shown in Fig. 27. One can find that with the increase of the
damping ratio the amplitude of vibration decreases accordingly, and the double-sided impact would
vanish gradually (see Fig. 27 and Fig. 28); moreover, the increase of the damping ratio could control
the behavior of impact. It can also be observed that the chaotic and subharmonic resonances may exist
depending upon the damping.

<figure><img src="image_33.png" /><figcaption>Fig. 27. Effect of damping ratio on Dynamic response for the forward motion
• ζ = 0.01, • ζ = 0.03, • ζ = 0.05</figcaption></figure>
<figure><img src="image_34.png" /><figcaption>Fig. 28. Bifurcation diagram to show effect of damping ratio on Dynamic response during the forward motion; ● ζ = 0.01, ● ζ = 0.03, ● ζ = 0.05</figcaption></figure>

**Conclusion**

This paper investigates the nonlinear dynamic behavior of a spiral bevel gear pair for the two-
directions of rotation. The calculation of an accurate mesh stiffness is the first step of the study. The
statical results are carried out by executing the finite element method for the full spiral bevel gear 3D
model. Preliminary static results, i.e., root stress, force on the tooth surface, and pressure distribution
on the tooth surface are illustrated. Due to the difference of teeth trace, two different mesh stiffness are
extracted and applied into the dynamic model. The dynamic model, which is nonlinear and time-
dependent parametrically excited system, is analyzed by using an implicit Runge–Kutta algorithm. The
dynamic responses of the system are presented in two main forms: amplitude-frequency and bifurcation
diagrams, for both forward and reverse motions.

The results show that the system experienced three types of contact, i.e., drive side contact, separation tooth, and backside contact. One of the most important parameters having a remarkable effect on the dynamic response of the system is the backlash, which leads to nonlinearity in the dynamic model of the system. Four different values for transverse backlash are considered in this study: big backlash (0.04 in), where backside contact is not observed; small backlash (0.01 in, 0.02 in, and 0.03 in), where unwanted backside contact occurs. The results show that the vibration amplitude increases by increasing the backlash. Eventually, tooth separation occurs for all backlash values for some specific frequency
ranges. These results are important from the practical point of view to understand how to avoid tooth separation and backside contacts.

The effect of damping ratio on the gear dynamic characteristics is also investigated. The results show that the damping induced by the lubrication effectively reduces gear vibration in resonance and super-harmonic resonance conditions. A lubricant with a higher viscosity gives a larger damping ratio and it diminished the vibration more significantly.

### CRediT authorship contribution statement

**M. Molaie:** Conceptualization, Methodology, Software, Validation, Formal analysis, Investigation, Resources, Data Curation, Writing - Original Draft, Writing - Review & Editing, Visualization. **F. S. Samani:** Methodology, Validation, Formal analysis, Investigation, Writing - Original Draft, Writing - Review & Editing, Supervision. **A. Zippo:** Methodology, Investigation, Resources, Data Curation, Writing - Original Draft, Writing - Review & Editing, Supervision, Project administration, Funding acquisition. **F. Pellicano:** Conceptualization, Methodology, Formal analysis, Investigation, Resources, Data Curation, Writing - Original Draft, Writing - Review & Editing, Supervision, Project administration, Funding acquisition.

### Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Appendix A. Important data of the analyzed spiral bevel gear pair

<table><caption>Table 4. System data [20, 28]</caption><tr><td>Hand of pinion</td><td>Left</td></tr><tr><td>Shaft offset</td><td>0</td></tr><tr><td>Shaft angle</td><td>90 Deg</td></tr><tr><td>Loaded side of gear</td><td>Convex</td></tr><tr><td>Driver</td><td>Pinion</td></tr><tr><td>Coefficient of friction</td><td>0</td></tr><tr><td>Pinion speed [rpm]</td><td>100</td></tr><tr><td>Diametral pitch [1/in]</td><td>5.141</td></tr><tr><td>Damping ratio [38], ζ</td><td>0.01</td></tr></table>

<table><caption>Table 5. Pinion and gear data [20, 28]</caption><thead><tr><th></th><th>Pinion</th><th>Gear</th></tr></thead><tbody><tr><td>No. of teeth</td><td>12</td><td>36</td></tr><tr><td>Transverse Circular Tooth Thickness at Pitch Cone [in]</td><td>0.32</td><td>0.15</td></tr><tr><td>Outer cone distance [in]</td><td>3.691</td><td>3.691</td></tr><tr><td>Face Width [in]</td><td>1.0</td><td>1.0</td></tr><tr><td>Face Angle [deg]</td><td>22.31667</td><td>72.5</td></tr><tr><td>Back Angle [deg]</td><td>18.433</td><td>71.5666</td></tr><tr><td>Spiral Angle [deg]</td><td>35</td><td>35</td></tr><tr><td>Pitch Angle [deg]</td><td>18.433</td><td>71.5666</td></tr><tr><td>Diameter of cylinder at the base of the tooth [in]</td><td>1.138</td><td>61.5</td></tr><tr><td>Young's Modulus [psi]</td><td>30.0×10<sup>6</sup></td><td>30.0×10<sup>6</sup></td></tr><tr><td>Poisson's Ratio</td><td>0.3</td><td>0.3</td></tr></tbody></table>
<table><thead><tr><th>Machine settings</th><th>Concave</th><th>Convex</th><th>Concave</th><th>Convex</th></tr></thead><tbody><tr><td>Radial setting [in]</td><td>2.947802</td><td>2.801049</td><td>2.85995</td><td>2.85995</td></tr><tr><td>Blank offset [in]</td><td>0.1545759</td><td>-0.1742616</td><td>0</td><td>0</td></tr><tr><td>Root Angle [deg]</td><td>16.8666</td><td>16.8666</td><td>67.68333</td><td>67.68333</td></tr><tr><td>Machine Center to Back [in]</td><td>-0.04023062</td><td>0.05414291</td><td>0</td><td>0</td></tr><tr><td>Sliding Base [in]</td><td>0.01167273</td><td>-0.01570932</td><td>0</td><td>0</td></tr><tr><td>Cradle Angle [deg]</td><td>63.94203</td><td>53.92599</td><td>59.2342023</td><td>59.2342023</td></tr><tr><td>Ratio of Roll</td><td>3.242698536</td><td>3.105176807</td><td>1.051674445</td><td>1.051674445</td></tr><tr><th>Cutter geometry</th><th>Concave</th><th>Convex</th><th>Concave</th><th>Convex</th></tr><tr><td>Cutter type</td><td>Straight</td><td>Straight</td><td>Straight</td><td>Straight</td></tr><tr><td>Point radius [in]</td><td>2.965621</td><td>3.071306</td><td>3.0325</td><td>2.9675</td></tr><tr><td>Blade angle [deg]</td><td>18.6015</td><td>24.90</td><td>22.0</td><td>22.0</td></tr><tr><td>Edge radius [in]</td><td>0.045</td><td>0.045</td><td>0.01</td><td>0.01</td></tr></tbody></table>

## References

1. M. Buzzoni, G. D'Elia, E. Mucchi, G. Dalpiaz, A vibration-based method for contact pattern assessment in straight bevel gears, Mechanical Systems and Signal Processing. 120 (2019) 693-707. https://doi.org/10.1016/j.ymssp.2018.10.002.

2. M. Molaie, F.S. Samani, F. Pellicano, Spiral Bevel Gears Nonlinear Vibration Having Radial and Axial Misalignments Effects, Vibration. 4 (3) (2021) 666-678. https://doi.org/10.3390/vibration4030037.

3. H. Motahar, F.S. Samani, M. Molaie, Nonlinear vibration of the bevel gear with teeth profile modification, Nonlinear Dynamics. 83 (4) (2016) 1875-1884. https://doi.org/10.1007/s11071-015-2452-z

4. P. McFadden, Determining the location of a fatigue crack in a gear from the phase of the change in the meshing vibration, Mechanical Systems and Signal Processing. 2 (4) (1988) 403-409. https://doi.org/10.1016/0888-3270(88)90063-5.

5. H.N. Özgüven, D.R. Houser, Mathematical models used in gear dynamics—a review, Journal of sound and vibration. 121 (3) (1988) 383-411. https://doi.org/10.1016/S0022-460X(88)80365-1.

6. S.D. Yavuz, Z.B. Saribay, E. Cigeroglu, Nonlinear time-varying dynamic analysis of a spiral bevel geared system, Nonlinear Dynamics. 92 (4) (2018) 1901-1919. https://doi.org/10.1007/s11071-018-4170-9.

7. D. Yassine, H. Ahmed, W. Lassaad, H. Mohamed, Effects of gear mesh fluctuation and defaults on the dynamic behavior of two-stage straight bevel system, Mechanism and Machine Theory. 82 (2014) 71-86. https://doi.org/10.1016/j.mechmachtheory.2014.07.013

8. A. Masoumi, F. Pellicano, F.S. Samani, M. Barbieri, Symmetry breaking and chaos-induced imbalance in planetary gears, Nonlinear Dynamics. 80 (1) (2015) 561-582. https://doi.org/10.1007/s11071-014-1890-3.

9. T. Kiekbusch, D. Sappok, B. Sauer, I. Howard, Calculation of the combined torsional mesh stiffness of spur gears with two-and three-dimensional parametrical FE models, Strojniški vestnik-Journal of Mechanical Engineering. 57 (11) (2011) 810-818. https://doi.org/10.5545/sv-jme.2010.248

10. J.-y. Tang, Z.-h. Hu, L.-j. Wu, S.-y. Chen, Effect of static transmission error on dynamic responses of spiral bevel gears, Journal of Central South University. 20 (3) (2013) 640-647. https://doi.org/10.1007/s11771-013-1530-y.

11. X. Su, M.M. Tomovic, D. Zhu, Diagnosis of gradual faults in high-speed gear pairs using machine learning, Journal of the Brazilian Society of Mechanical Sciences and Engineering. 41 (4) (2019) 1-11. https://doi.org/10.1007/s40430-019-1701-3.

12. L. Gelman, N.H. Chandra, R. Kurosz, F. Pellicano, M. Barbieri, A. Zippo, Novel spectral kurtosis technology for adaptive vibration condition monitoring of multi-stage gearboxes, Insight-Non-Destructive Testing and Condition Monitoring. 58 (8) (2016) 409-416. https://doi.org/10.1784/insi.2016.58.8.409.

13. S. Peng, H. Ding, G. Zhang, J. Tang, Y. Tang, New determination to loaded transmission error of the spiral bevel gear considering multiple elastic deformation evaluations under different bearing supports, Mechanism and Machine Theory. 137 (2019) 37-52. https://doi.org/10.1016/j.mechmachtheory.2019.03.013.
14. Q. Wang, K. Chen, B. Zhao, H. Ma, X. Kong, An analytical-finite-element method for calculating mesh stiffness of spur gear pairs with complicated foundation and crack, Engineering Failure Analysis. 94 (2018) 339-353. https://doi.org/10.1016/j.engfailanal.2018.08.013.

15. M. Barbieri, A. Zippo, F. Pellicano, Adaptive grid-size finite element modeling of helical gear pairs, *Mechanism and Machine Theory*. 82 (2014) 17-32. https://doi.org/10.1016/j.mechmachtheory.2014.07.009.

16. F. Zheng, X. Han, H. Lin, W. Zhao, Research on the cutting dynamics for face-milling of spiral bevel gears, *Mechanical Systems and Signal Processing*. 153 (2021) 107488. https://doi.org/10.1016/j.ymssp.2020.107488.

17. J. Klingelnberg, Bevel gear, Springer Vieweg Berlin, Heidelberg, 2016. https://doi.org/10.1007/978-3-662-43893-0.

18. Q. Fan, Enhanced algorithms of contact simulation for hypoid gear drives produced by face-milling and face-hobbing processes, *Journal of Mechanical Design*. 129 (1) (2007) 31-37. https://doi.org/10.1115/1.2359475.

19. Q. Fan, Ease-Off and Application in Tooth Contact Analysis for Face-Milled and Face-Hobbed Spiral Bevel and Hypoid Gears, In: Goldfarb, V., Barmina, N. (eds) Theory and Practice of Gearing and Transmissions, Mechanisms and Machine Science. 34, 2016, pp.321-339. https://doi.org/10.1007/978-3-319-19740-1_15.

20. S. Vijayakar, S. Abad, R. Gunda, Multi-Body Dynamic Contact Analysis, Tool for Transmission Design SBIR Phase II Final Report, Advanced Numerical Solutions. Ohio, 2003.

21. M. Faggioni, F.S. Samani, G. Bertacchi, F. Pellicano, Dynamic optimization of spur gears, *Mechanism and machine theory*. 46 (4) (2011) 544-557. https://doi.org/10.1016/j.mechmachtheory.2010.11.005.

22. M. Molaie, F.S. Samani, H. Motahar, Nonlinear vibration of crowned gear pairs considering the effect of Hertzian contact stiffness, *SN Applied Sciences*. 1 (5) (2019) 414. https://doi.org/10.1007/s42452-019-0439-y.

23. S. Nie, J. Deng, X. Deng, L. Geng, A flank modification method for spiral bevel gears based on mismatch topography adjustment, *Journal of Advanced Mechanical Design, Systems, and Manufacturing*. 12 (2) (2018). https://doi.org/10.1299/jamdsm.2018jamdsm0057.

24. R. Handschuh, T. Kicher, A method for thermal analysis of spiral bevel gears, *Journal of Mechanical Design*. 118 (4) (1996) 580-585. https://doi.org/10.1115/1.2826932.

25. V. Simon, Optimal modifications of gear tooth surfaces, *Gear Technology*. 3 (4) (2011) 62-72.

26. R.F. Handschuh, Recent advances in the analysis of spiral bevel gears, in: *Proceedings of MTM'97 International Conference on Mechanical Transmissions and Mechanisms*. Tianjin-China, (1997) 635-641.

27. R.F. Handschuh, Testing of face-milled spiral bevel gears at high-speed and load, in: *Proceedings of the International Conference on Mechanical Transmissions. NASA/TM-2001-210743*, Chongqing-China, 2001.

28. R.F. Handschuh, F.L. Litvin, A method for determining spiral-bevel gear tooth geometry for finite element analysis. *NASA Technical Paper 3096*, AVSCOM Technical Report 91-C-020, 1991.

29. R.F. Handschuh, Comparison of experimental and analytical tooth bending stress of aerospace spiral bevel gears, in: *Proceedings of the 4th World Congress on Gearing and Power Transmission*. NASA/TM-1999-208903, Paris-France, 1999.

30. R.F. Handschuh, M. Nanlawala, J.M. Hawkins, D. Mahan, Experimental comparison of face-milled and face-hobbed spiral bevel gears, in: *Proceedings of the International Conference on Motion and Power Transmission. NASA/TM-2001-210940*, Fukuoka-Japan, 2001.

31. G.D. Bibel, K. Tiku, A. Kumar, Prediction of Contact Path and Load Sharing in Spiral Bevel Gears, Lewis Research Center, NASA CR-195305, ARL-CR-146, 1994.

32. F.S. Samani, M. Salajegheh, M. Molaie, Nonlinear vibration of the spiral bevel gear under periodic torque considering multiple elastic deformation evaluations due to different bearing supports, *SN Applied Sciences*. 3 (9) (2021) 1-9. https://doi.org/10.1007/s42452-021-04755-6.

33. F.S. Samani, M. Molaie, F. Pellicano, Nonlinear vibration of the spiral bevel gear with a novel tooth surface modification method, *Meccanica*. 54 (7) (2019) 1071-1081. https://doi.org/10.1007/s11012-019-00973-w.
34. G. Liu, R.G. Parker, Nonlinear dynamics of idler gear systems, Nonlinear Dynamics. 53 (4) (2008) 345-367. https://doi.org/10.1007/s11071-007-9317-z.

35. A. Kahraman, R. Singh, Non-linear dynamics of a spur gear pair, Journal of sound and vibration. 142 (1) (1990) 49-75. https://doi.org/10.1016/0022-460X(90)90582-K.

36. G. Bonori, Static and dynamic modeling of gear transmission error, Ph. D. Thesis, University of Modena and Reggio Emilia. Italy, 2006.

37. R.K. Mobley, Gears and Gearboxes, in: Plant Engineer's Handbook, Butterworth-Heinemann. 2001, pp. 629-637.

38. G. Blankenship, A. Kahraman, Steady state forced response of a mechanical oscillator with combined parametric excitation and clearance type non-linearity, Journal of Sound and Vibration. 185 (5) (1995) 743-765. https://doi.org/10.1006/jsvi.1995.0416.