Moslem Molaie
email: moslem_molaie@unimore.it

Associate Professor Farhad S Samani
email: farhad.samani@uk.ac.ir

Antonio Zippo
email: antonio.zippo@unimore.it

Giovanni Iarriccio
email: giovanni.iarriccio@unimore.it

Francesco Pellicano
email: francesco.pellicano@unimore.it

Spiral Bevel Gears: bifurcation and chaos analyses of pure torsional system

Keywords: Spiral bevel gear, Mesh stiffness, Torsional shaft stiffness, Bifurcation and chaos analyses, nonlinear time series analysis

The goal of this paper is to provide a detailed analysis of the complex dynamic scenario of spiral bevel gears by considering the torsional shaft stiffness. The dynamic model takes into account: timedependent stiffness and non-smooth nonlinearity due to the backlash, i.e., teeth contact loss. The time-varying meshing stiffness is evaluated by means of a nonlinear finite element model, which allowed an accurate evaluation of the global and local teeth deformation in both directions: forward and reverse motions. Due to intentional tooth profile errors, e.g., profile modification, manufacturing error, or assembly error, the loaded and unloaded tooth contact analyses were conducted. The dynamic model was validated by comparisons with a verified SDOF model in terms of linear natural frequencies and nonlinear dynamic response. The present study provides amplitude frequency diagrams and bifurcation diagrams; for specific regimes, periodic, quasiperiodic, and chaotic responses were found. The dynamic behavior of the systems was evaluated using various tools such as modal analysis, nonlinear time series analysis, spectra, 3D-phase diagrams, and Poincaré maps. This study proves that decreasing the DOFs of the system can lose some dynamics. Different phenomena have been found such as trapping and boom-and-bust cycle. The dynamic response of different cases is evaluated by estimating the largest Lyapunov exponent and the correlation dimension. It is found that the system undergoes complex dynamic phenomena including nT-periodic, trapping, and aperiodic motions, which are evidenced in 3D-bifurcation diagrams, spectra, Poincaré maps, and 3Dphase diagrams. The chaotic motion found in this study is an undesirable behavior from a practical point of view, which can be avoided by adjusting the designed parameters. The results provide a basis for parameter design and dynamic characteristic control of the spiral bevel gear drive system.

Introduction

The dynamic behavior of gear systems is a topic of great interest in mechanical transmissions.

Both the lifetime and the noise generated by a gear set can be influenced by its dynamic response.

The presence of spiral bevel gears (SBGs) in many examples of high-power transmission systems has drawn attention to the importance of comprehending their dynamic behavior. SBGs are often utilized when a high-power density is required between intersected shafts. The complex geometry of SBGs makes it necessary to carefully consider the parameters that have a significant impact on the system response, for example the backlash and the gear tooth mesh stiffness (MS) which may lead to a nonsmooth time-varying dynamic system.

It is interesting to note that the dynamic models for the gear systems may have great variety; however, it may still be possible to achieve similar predictions by employing completely different models for certain systems [START_REF] Özgüven | A non-linear mathematical model for dynamic analysis of spur gears including shaft and bearing dynamics[END_REF]. This depends on the dynamic characteristics of the systems under investigation; for instance, it may be sufficient to use a single-degree-of-freedom (SDOF) model to obtain accurate prediction, when the shaft torsional and flexural stiffnesses are decoupled from the other modes. In such cases, using a multi-degree-of-freedom model that accounts for the coupling between the torsional and transverse vibrations will not significantly improve the accuracy. On the contrary, when the dynamic properties lead to a strong coupling among meshing modes and other modes, the SDOF model cannot be significantly accurate. In such cases, no matter how sophisticated the SDOF is (considering non-linearity of the elements, excitation due to gear errors, time variation of MS, and damping), it cannot be accurate enough for response predictions.

Yassine et al. [START_REF] Yassine | Effects of gear mesh fluctuation and defaults on the dynamic behavior of two-stage straight bevel system[END_REF] developed a three-dimensional model of a two-stage straight bevel gear system and conducted an investigation of its dynamic response. The main source of excitation was deemed to be the periodic fluctuations of the gear MS. They made a comparison between the behavior of a defective gear system (i.e., runout, profile error, tooth crack) and that of a flawless system. Their results have implications in the design and maintenance of straight bevel gear transmissions. In 2018, Yavuz et al. [START_REF] Yavuz | Nonlinear time-varying dynamic analysis of a spiral bevel geared system[END_REF] carried out a study on a dynamic model of a drivetrain consisting of a spiral bevel gear pair, shafts, and bearings. The dynamic model included the flexibilities of the shaft bearings, the backlash, and the time-varying MS. The multi-term harmonic balance method was applied to solve the governing equations, in conjunction with Newton's method with arc-length continuation. They investigated the impact of backlash, fluctuation of gear MS, and variation of bearing stiffness on the dynamic scenario. Sun et al. [START_REF] Sun | On Dynamic Mesh Force Evaluation of Spiral Bevel Gears[END_REF] developed a 12-DOF model, besides a tooth pair mesh model proposed to improve describing mesh characteristics of individual tooth pairs. Three different approaches were compared to investigate the effect of mesh stiffness representation, including the average slope approach, local slope approach, and a quasi-statically defined interpolation function. In 2020, a dynamic model for the nonlinear time-varying behavior of a drivetrain consisting of parallel and intersecting axes was presented by Yavuz et al. [START_REF] Yavuz | Nonlinear dynamic analysis of a drivetrain composed of spur, helical and spiral bevel gears[END_REF]. Shaft and bearing flexibilities are included in the model through the use of finite element modeling and mesh models that account for backlash nonlinearity and fluctuating MS.

A significant part of the scientific literature has been dedicated to developing dynamic models for gear pairs in order to consider the influence of specific parameter parts, e.g., backlash [START_REF] Huang | Bifurcation and chaos analysis of a spur gear pair system with fractal gear backlash[END_REF],

manufacturing error [START_REF] Bonori | Non-smooth dynamics of spur gears with manufacturing errors[END_REF], external load [START_REF] Theodossiades | Periodic and chaotic dynamics of motor-driven gear-pair systems with backlash[END_REF], and shaft flexibility [START_REF] Litak | Vibration in gear systems[END_REF]. Shi and Li [START_REF] Shi | Nonlinear dynamics of hypoid gear with coupled dynamic mesh stiffness[END_REF] proposed two different dynamic models: 4-DOF and 10-DOF, for hypoid gears that consider MS dependent on the dynamic mesh force. The model takes into account the interaction between MS and dynamic response; a parametric surface approach was used to overcome the computational challenges of determining the dynamic mesh stiffness. The study showed that dynamic mesh stiffness has a significant impact on hypoid gear dynamic response, especially at lower loads, and provided a physicsbased tool for understanding the interaction between mesh stiffness and dynamic mesh force. In 2022, a torsional-bending-axial dynamic model (12-DOFs) of the straight bevel gear system was proposed by Gou et al. [START_REF] Gou | Dynamic characteristics of a straight bevel gear drive system considering multi-state meshing and time-varying parameters[END_REF] to investigate the dynamic characteristics. Their research provided a theoretical basis for the parametric design and dynamic characteristics control of straight bevel gear systems.

The main goal of the present study is to investigate the nonlinear dynamic scenario of spiral bevel gears, considering the effect of torsional shaft stiffness on the dynamic response of a real model of the SBG used in the transmission power system of helicopters. Shafts can have different lengths and diameters, which significantly affects the dynamic response of the system. In this study, the possibility of backside contact is considered, which leads to calculating the tooth-MS in both the forward and the reverse motions due to the geometry of the SBG. A nonlinear finite element method-based program, Transmission3D Calyx, is used to perform a precise static analysis to determine the system's MS.

Because of the mismatch between the tooth surface of gear and pinion, a loaded tooth contact analysis (LTCA) and an unloaded tooth contact analysis (UTCA) are conducted. The dynamic model is used for investigating the nonlinear dynamic scenario; the analysis is carried out through a three-DOF model considering: torsional shaft stiffness, the effect of backlash, which leads to the non-smooth equation, variable stiffness, i.e., time-varying coefficients equation and consequent parametric excitation and the effect of stiffness variation and phasing in the case of backside contact (reverse motion). For a deeper understanding of the torsional shaft stiffness effect, a reduced single degree of freedom is compared with a system with three-DOF. Besides, the effect of shaft stiffness was investigated from a design point of view. The dynamic behavior of the systems has been evaluated through different tools: modal analysis, amplitude-frequency diagrams, bifurcation diagrams, timehistory responses, spectra, 3D-phase diagrams, Poincaré maps, largest Lyapunov exponent, and correlation dimension.

Dynamic system

Fig. 1 shows a scheme of a spiral bevel gear system, where mechanical power is transmitted between two orthogonal axes. A purely torsional dynamic model is used, i.e., all translational DOFs are constrained and both the gear and the pinion rotate only around their axes.

System with three-DOF

To consider the effect of torsional shaft stiffness on the dynamic response of the system, the governing equations of motion are derived for the system with three-DOF, see Fig. 2. The dynamic equations of motion in "𝜃" and "𝜑" direction for the gear (𝜃 ! ) and the pinion (𝜑 " ) are:

𝜃 ! -DOF: 𝐼 ! " 𝜃 ̈! = -𝑟 ! 𝐹 # -𝐾 $ ! /𝜃 ! -𝜃 % 0 -𝐶 $ ! /𝜃 ̇! -𝜃 ̇%0 (1) 𝜑 & -DOF: 𝐼 & ' 𝜑̈& = 𝑟 & 𝐹 # -𝐾 $ & /𝜑 & -𝜑 ( 0 -𝐶 $ & /𝜑̇& -𝜑̇(0 (2) Pinion Gear φ 𝑇 ! 𝜃 𝑇 # x y z 𝐾 ! (𝑡) Pinion Gear 𝑒(𝑡) 𝑟 $ 𝑟 $ 𝐶 ! 𝜃 % 𝜃 # 𝜑 & 𝜑 ! Motor Load 2b K ' ( K ' ) C ' ) C ' (
The torsional shaft stiffness is:

8 8 𝐾 $ )) = 𝐺 𝐽 𝑙 𝜏 (*' = 𝑇 • 𝑑 2𝐽 %+,- @AAB 𝐾 $ )) = 𝑇 ) • 𝐺 10 • 𝜏 , (𝑖 = 𝑚, 𝑙; 𝑖𝑖 = 𝑝, 𝑔) (3)
Where, 𝜏 𝑚𝑎𝑥 is the maximum torsional stress of the shaft, 𝐺 is the shear modulus, 𝐽 is the second polar moment of shaft area, 𝑟 " and 𝑟 ! are the mean radii at the meshing points, 𝐶 ' " and 𝐶 ' ! are the material damping of the pinion and gear shafts. The gear system is simulated taking the motor and load masses into consideration. Their relative governing equations in the "𝜑" and "𝜃" directions are respectively given for the motor (𝜑 ( ) and load (𝜃 # ) by:

𝜑 ( -DOF: 𝐼 ( 𝜑̈( = 𝑇 ( -𝐾 $ & /𝜑 ( -𝜑 & 0 -𝐶 $ & /𝜑̇( -𝜑̇&0 (4)
𝜃 % -DOF:

𝐼 % 𝜃 ̈% = -𝑇 % -𝐾 $ ! /𝜃 % -𝜃 ! 0 -𝐶 $ ! /𝜃 ̇% -𝜃 ̇!0 (5)
Where, 𝑇 % = /N / N 0 ⁄ 0 • 𝑇 1 , N 0 and N / are the pinion and gear teeth number respectively. Due to mounting, manufacturing error (Fig. 3), or intentional teeth profile modifications, a local gap between mating teeth can appear, it is called geometric transmission error (GTE), which is here calculated through a UTCA for forward and reverse motions, 𝑒 * (𝑡) and 𝑒 + (𝑡) respectively:

𝑒(𝑡) = ⎩ ⎪ ⎨ ⎪ ⎧ 𝑒 𝐹 (𝑡) = 𝑒 2 3 + V 𝑒𝑎
The linear dynamic transmission error (DTE), 𝜆, along the line of action is defined as:

𝜆 = /𝑟 & 𝜑 & -𝑟 ! 𝜃 ! 0𝑎 (7)
where, 𝑎 = cos 𝛼 × cos 𝛽, β is the spiral angle and α is the normal pressure angle. The dynamic load of pinion and its component along the line of action is 𝐹 , and along the 𝑧-axis is 𝐹 -which are defined as:

𝐹 # = -𝐹 8 𝑎; 𝐹 8 = 𝐶 ( /𝜆 ̇-𝑒0 + 𝐾 ( (𝑡)𝑓(𝜆 -𝑒) (8)
In ( 8) the nominal force is given by a linear and an elastic nonlinear term, where:

𝑓(𝜆 -𝑒) = d 𝜆 -𝑏 -𝑒, 𝜆 -𝑒 > 𝑏 0, -𝑏 ≤ 𝜆 -𝑒 ≤ 𝑏 𝜆 + 𝑏 -𝑒, 𝜆 -𝑒 < -𝑏 ; 𝐾 ( (𝑡) = i 𝐾 ( 3 (𝑡), 𝜆 -𝑒 > 𝑏 𝐾 ( 7 (𝑡), 𝜆 -𝑒 < -𝑏 (9)
𝑓(𝜆 -𝑒) is the backlash function that is piecewise linear; see Fig. 4. This function multiplied by the stiffness returns the elastic restoring force [START_REF] Samani | Nonlinear vibration of the spiral bevel gear with a novel tooth surface modification method[END_REF]. Whenever 𝜆 -𝑒 is between -𝑏 and +𝑏, the contact loss happens [START_REF] Liu | Nonlinear dynamics of idler gear systems[END_REF], [START_REF] Kahraman | Non-linear dynamics of a spur gear pair[END_REF]. For 𝜆 -𝑒 > 𝑏, the mesh is expected to be in the forward contact (desired situation), where 𝐾 𝑚 𝐹 (𝑡) is extracted from static nonlinear FEA during the forward motion.

While if 𝜆 -𝑒 < -𝑏, undesired backside contact happens (double-sided impact), where 𝐾 𝑚 𝑅 (𝑡) is obtained from static simulation during the reverse motion; see Ref. [START_REF] Kahraman | Non-linear dynamics of a spur gear pair[END_REF]. The results of the experiments given in Ref. [START_REF] Al-Shyyab | Non-linear dynamic analysis of a multi-mesh gear train using multiterm harmonic balance method: sub-harmonic motions[END_REF] revealed that tooth separation is likely to take place, which suggests that gear backlash non-linearity needs to be considered in the dynamic models. Furthermore, it was clarified that the mesh stiffness fluctuation is one of the predominant sources of excitation from the gear mesh and must therefore be considered in the governing equations. 𝑡) [START_REF] Shi | Nonlinear dynamics of hypoid gear with coupled dynamic mesh stiffness[END_REF] Where 𝜔 ( is the fundamental frequency of 𝐾 ( , technically meshing frequency, 𝑁 5 is the number of samples considered within a meshing period when the 𝐾 ( is numerically analyzed by nonlinear FEA.

In order to normalize the governing equation, the following parameters are introduced. Consequently, Eqs. ( 1)-( 5) can be rewritten as follows:

𝜏 = 𝜔 8 𝑡 𝜆 ̅ = /𝑟̅ & 𝜑 & -𝑟̅ ! 𝜃 ! 0𝑎 : + 𝑒̅ (𝑡) 𝜃 ̈= 𝜔 8 ; 𝜃 << 𝜃 ̇= 𝜔 8 𝜃 < (11) 𝑟̅ & = 𝑟 & 𝑏 s 𝑟̅ ! = 𝑟 ! 𝑏 s 𝜑̈= 𝜔 8 ; 𝜑 << 𝜑̇= 𝜔 8 𝜑 < 𝐾 t $ ! = 𝐾 $ ! 𝐼 ! " 𝜔 8 ; 𝐶 ̅ $ ! = 𝐶 $ ! 2𝐼 ! " 𝜔 8 ( 12
) 𝐾 t $ & = 𝐾 $ & 𝐼 & ' 𝜔 8 ; 𝐶 ̅ $ & = 𝐶 $ & 2𝐼 & ' 𝜔 8 ( 13
) 𝐾 t $ ( = 𝐾 $ & 𝐼 ( 𝜔 8 ; 𝐶 ̅ $ ( = 𝐶 $ & 2𝐼 ( 𝜔 8 𝑇 u ( = 𝑇 ( 𝐼 ( 𝜔 8 ; ( 14
) 𝐾 t $ % = 𝐾 $ ! 𝐼 % 𝜔 8 ; 𝐶 ̅ $ % = 𝐶 $ ! 2𝐼 % 𝜔 8 𝑇 u % = 𝑇 % 𝐼 % 𝜔 8 ; (15
𝜃 ! << = 𝑚 => 𝑏𝑟 ! 𝑎 : 𝐼 ! " /𝐾 t ( (𝑡)𝑓/𝜆 ̅ -𝑒̅ 0 + 2𝜉𝜆 ̅ < 0 + 𝐾 t $ ! /𝜃 % -𝜃 ! 0 + 2𝐶 ̅ $ ! /𝜃 % < -𝜃 ! < 0 (16) 𝜑 & << = 𝐾 t $ & /𝜑 ( -𝜑 & 0 + 2𝐶 ̅ $ & /𝜑 ( < -𝜑 & < 0 - 𝑚 => 𝑏𝑟 & 𝑎 : 𝐼 & ' /𝐾 t ( (𝑡)𝑓/𝜆 ̅ -𝑒̅ 0 + 2𝜉𝜆 ̅ < 0 (17) 𝜑 ( << = 𝑇 u ( + 𝐾 t $ ( /𝜑 & -𝜑 ( 0 + 2𝐶 ̅ $ ( /𝜑 & < -𝜑 ( < 0 (18) 𝜃 % << = 𝐾 t $ % /𝜃 ! -𝜃 % 0 + 2𝐶 ̅ $ % /𝜃 ! < -𝜃 % < 0 -𝑇 u % (19)
where,

𝑓(λ u -e u) = i λ u -e u -1, λ u -e u > 1 0, -1 ≤ λ u -e u ≤ 1 𝜆 ̅ -e u + 1, λ u -e u < -1 (20) 𝐾 ( (𝑡) = 1 + V 𝑎 4 𝑚 ?@ 𝜔 8 ; cos (𝑗𝜔 ( 9 4+6 𝑡) + V 𝑏 4 𝑚 ?@ 𝜔 8 ; sin (𝑗𝜔 ( 9 4+6 𝑡) (21)
Typically, in gearboxes the driven shaft is connected to big inertial masses; therefore, 𝜃 𝑙 can be neglected in Eqs. ( 16)- [START_REF] Molaie | Spiral Bevel Gears: nonlinear dynamic model based on accurate static stiffness evaluation[END_REF]. The four-DOF system is then reduced to a three-DOF system as represented in Eqs. ( 22)-( 24):

𝜃 ! << = 𝑚 => 𝑏𝑟 ! 𝑎 : 𝐼 ! " /𝐾 t ( (𝑡)𝑓/𝜆 ̅ -𝑒̅ 0 + 2𝜉𝜆 ̅ < 0 -𝐾 t $ ! /𝜃 ! 0 -2𝐶 ̅ $ ! /𝜃 ! < 0 (22) 𝜑 & << = 𝐾 t $ & /𝜑 ( -𝜑 & 0 + 2𝐶 ̅ $ & /𝜑 ( < -𝜑 & < 0 - 𝑚 => 𝑏𝑟 & 𝑎 : 𝐼 & ' /𝐾 t ( (𝑡)𝑓/𝜆 ̅ -𝑒̅ 0 + 2𝜉𝜆 ̅ < 0 (23)
𝜑 ( << = 𝑇 u ( + 𝐾 t $ ( /𝜑 & -𝜑 ( 0 + 2𝐶 ̅ $ ( /𝜑 & < -𝜑 ( < 0 (24)

System with single DOF

To investigate the effect of shaft stiffness on the dynamic behavior of the system, the three-DOF model is compared with the SDOF system. The dynamic model of the SDOF system is represented in Fig. 5. If the shaft compliance can be neglected, i.e., they are extremely stiff, then 𝜃 ! -𝜃 # ≈ 0 and 𝜑 " -𝜑 ( ≈ 0. Then the three-DOF system can be reduced to SDOF system in Fig. 5.

𝑟 $ 𝑟 $ 𝐶 ! 𝜃 % 𝜑 & 2b
Considering the aforementioned assumption on the shaft stiffness Eqs. (16-(19) can be reduced to:

𝐼 ! " 𝜃 ̈! = -𝑟 ! 𝐹 # -𝑇 % (25) 𝐼 & ' 𝜑̈& = 𝑟 & 𝐹 # + 𝑇 ( (26)
Following Refs. [START_REF] Samani | Nonlinear vibration of the spiral bevel gear with a novel tooth surface modification method[END_REF], [START_REF] Samani | Nonlinear vibration of the spiral bevel gear under periodic torque considering multiple elastic deformation evaluations due to different bearing supports[END_REF]- [START_REF] Molaie | Spiral Bevel Gears: nonlinear dynamic model based on accurate static stiffness evaluation[END_REF], we can introduce the parameter 𝜆 of Eq. ( 11) and the following parameters:

𝑇 u => = 1 𝑏𝑚 => 𝜔 8 ; y 𝑇 ( 𝑟 & 𝑎 : z 𝜔 8 = { 𝐾 2 𝑚 => s (27) 𝜉 = 𝐶 ( 2𝑚 => 𝜔 8 𝑚 => = | /𝑟 & 𝑎0 ; 𝐼 & ' + /𝑟 ! 𝑎0 ; 𝐼 ! " } A6
Then the rigid body motion is removed to get the SDOF model, in terms of the DTE (𝜆). The governing equations ( 25)and ( 26) are combined to obtain the following:

𝜆 ̅ << + 𝐾 ( (𝑡)𝑓/𝜆 ̅ -𝑒̅ 0 + 2𝜉𝜆 ̅ < = 𝑇 u => (28)
It must be noted that the backlash function 𝑓G𝜆 ̅ -𝑒̅ I, the mesh stiffness 𝐾 ( (𝑡), and DTE 𝜆 ̅ are the same for both the SDOF and the three-DOF systems.

Mesh stiffness calculation

One of the fundamental steps needed to carry out the vibration analysis, is the evaluation of the mesh stiffness through the loaded and unloaded contact analyses, namely LTCA and UTCA. LTCA and UTCA provide the transmission error under nominal and negligible loading conditions, respectively.

The difference between the transmission error obtained from the LTCA and the UTCA is the part of the transmission error related exclusively to the elastic deformation. This value, known as static transmission error, is used for calculation of the mesh stiffness.

SBGs are typically manufactured with an intentional tooth profile error, i.e., the mismatch between two teeth surfaces of gear and pinion [START_REF] Simon | Optimal modifications of gear tooth surfaces[END_REF]- [START_REF] Handschuh | A Method for Thermal Analysis of Spiral Bevel Gears[END_REF]. Such intentional profile error leads to effectively performs gear mesh even in the presence of misalignments, manufacturing errors, or high torque levels which cause the gears to deflect to imperfect locations. Although conjugate teeth profiles provide an ideal transmission with a constant transmission ratio, which is not practical for SBGs. Consequently, SBGs operate with a certain amount of geometric transmission error [START_REF] Simon | Optimal modifications of gear tooth surfaces[END_REF], [START_REF] Handschuh | A Method for Thermal Analysis of Spiral Bevel Gears[END_REF], [START_REF] Handschuh | Recent Advances in the Analysis of Spiral Bevel Gears[END_REF], and therefore LTCA and UTCA are of utmost importance in obtaining the correct mesh stiffness.

The LTCA and UTCA could be conducted using different approaches: theoretical models, experiments, and FEM:

1. Theoretical models provide solutions for simplified cases while, for complex gear geometry, such as the case of SBGs, closed form or approximate solutions are not accurate enough.

2. Experimental methods rely on well-equipped testing facilities. Despite their effectiveness, tests are time-demanding, expensive, and measurements on rotating parts are limited to non-contact surfaces.

3. FEM is well-established for providing affordable and reliable results within a reasonable time. More specifically, nonlinear FEM has become increasingly popular thanks to the possibility of including contact nonlinearities throughout the analysis [START_REF] Molaie | Nonlinear vibration of crowned gear pairs considering the effect of Hertzian contact stiffness[END_REF], [START_REF] Molaie | Spiral Bevel Gears: nonlinear dynamic model based on accurate static stiffness evaluation[END_REF], [START_REF] Motahar | Nonlinear vibration of the bevel gear with teeth profile modification[END_REF], [START_REF] Faggioni | Dynamic optimization of spur gears[END_REF].

In the present study, a nonlinear finite element software was used to perform tooth contact analysis. The kinematics of gear cutting was implemented in ANSOL-Calyx (Calyx-HypoidFaceMilled product form ANSOL company -Version 2011) to generate the spiral bevel gear geometry. A mesh sensitivity analysis was carried out to determine the appropriate mesh size (see details given in Ref.

[19]), and the simulation output was validated considering NASA experimental data [START_REF] Handschuh | A method for determining spiral-bevel gear tooth geometry for finite element analysis[END_REF].

The results of the loaded and unloaded tooth contact analysis are displayed in Figs. 6789., both the forward and the reverse motions, where red dots represent the mesh stiffness and transmission error evaluated from FEM analyses, and the continuous line represents the Fourier series used to fit the FEM results.

Considering the normalized meshing period, it is worth noticing that, in the LTCA, the teeth undergo deflection due to a combination of elastic deformation and rigid transmission error, see Fig. 6, while, during the UTCA, the teeth are subjected to a pure rotational displacement, see Fig. The static transmission error obtained from the difference between the LTCA and UTCA is shown in Fig. 8. Since the LTCA and UTCA are performed under constant torque, the knowledge of the static transmission error allows for the evaluation of the mesh stiffness over the meshing period, see Fig. 9.

Once the mesh stiffness is known, its approximated expression, obtained from interpolation, can be used in Eqs. ( 16)- [START_REF] Molaie | Spiral Bevel Gears: nonlinear dynamic model based on accurate static stiffness evaluation[END_REF].

Modal analysis

Modal analysis is one of the effective tools for identifying the mechanical properties within gear systems; for such analysis, a simplified linear time-invariant form of equations is employed (i.e., Eqs.

(1)-(5) for the system with three-DOF, and Eqs. ( 25)-(26) for the system with SDOF). The dynamic load along the line of action is 𝐹 , = 𝐾 ( 𝜆 + 𝐶 ( 𝜆 ̇, where 𝐾 ( is the mean value of the mesh stiffness, (𝐾 ( = 𝑘 3 ), and the backlash is considered exactly zero. Eqs. ( 1)-( 5) can be rewritten as follows:

𝜽 𝒈 -𝐃𝐎𝐅: 𝐼 ! " 𝜃 ̈! -𝑟 ! 𝑎 ; 𝐾 ( 𝑟 & 𝜑 & -𝑟 ! 𝑎 ; 𝐶 ( 𝑟 & 𝜑̇& + '𝐾 $ ! + 𝑟 ! 𝑎 ; 𝐾 ( 𝑟 ! ƒ𝜃 ! + '𝐶 $ ! + 𝑟 ! 𝑎 ; 𝐶 ( 𝑟 ! ƒ𝜃 ̇! -𝐾 $ ! 𝜃 % -𝐶 $ ! 𝜃 ̇% = 0 (29)
𝝋 𝒑 -𝐃𝐎𝐅:

𝐼 & ' 𝜑̈& + '𝐾 $ & + 𝑟 & 𝑎 ; 𝐾 ( 𝑟 & ƒ𝜑 & + '𝐶 $ & + 𝑟 & 𝑎 ; 𝐶 ( 𝑟 & ƒ𝜑̇& -𝑟 & 𝑎 ; 𝐾 ( 𝑟 ! 𝜃 ! -𝑟 & 𝑎 ; 𝐶 ( 𝑟 ! 𝜃 ̇! -𝐾 $ & 𝜑 ( -𝐶 $ & 𝜑̇( = 0 (30)
𝝋 𝒎 -𝐃𝐎𝐅:

𝐼 ( 𝜑̈( + 𝐾 $ & 𝜑 ( + 𝐶 $ & 𝜑̇( -𝐾 $ & 𝜑 & -𝐶 $ & 𝜑̇& = 𝑇 ( (31)
𝝋 𝒍 -𝐃𝐎𝐅:

𝐼 % 𝜃 ̈% + 𝐾 $ ! 𝜃 % + 𝐶 $ ! 𝜃 ̇% -𝐾 $ ! 𝜃 ! -𝐶 $ ! 𝜃 ̇! = -𝑇 % ( 32
)
Mesh stiffness [N/m] Normalized time steps Normalized time steps Mesh stiffness [N/m] Finally, we can obtain the stiffness and mass matrices:

[𝑀] = ⎣ ⎢ ⎢ ⎢ ⎡ 𝐼 & ' 0 0 0 0 𝐼 ! " 0 0 0 0 𝐼 ( 0 0 0 0 𝐼 % ⎦ ⎥ ⎥ ⎥ ⎤ ; [𝐾] = ⎣ ⎢ ⎢ ⎢ ⎡ 𝐾 $ & + 𝑟 & 𝑎 ; 𝐾 ( 𝑟 & -𝑟 & 𝑎 ; 𝐾 ( 𝑟 ! -𝐾 $ 0 -𝑟 ! 𝑎 ; 𝐾 ( 𝑟 & 𝐾 $ ! + 𝑟 ! 𝑎 ; 𝐾 ( 𝑟 ! 0 -𝐾 $ ! -𝐾 $ & 0 𝐾 $ & 0 0 -𝐾 $ ! 0 𝐾 $ ! ⎦ ⎥ ⎥ ⎥ ⎤ ; 𝑞 ⃑ = • 𝜑 & 𝜃 ! 𝜑 ( 𝜃 % ' (33)
Likewise, for the reduced system, Eqs. ( 25)-( 26) can be written as follows:

𝜽 𝒈 -𝐃𝐎𝐅: 𝐼 ! " 𝜃 ̈! -𝑟 ! 𝑎 ; 𝐾 ( 𝑟 & 𝜑 & -𝑟 ! 𝑎 ; 𝐶 ( 𝑟 & 𝜑̇& + '𝑟 ! 𝑎 ; 𝐾 ( 𝑟 ! ƒ𝜃 ! + '𝑟 ! 𝑎 ; 𝐶 ( 𝑟 ! ƒ𝜃 ̇! = -𝑇 % (34)
𝝋 𝒑 -𝐃𝐎𝐅:

𝐼 & ' 𝜑̈& + '𝐾 $ & + 𝑟 & 𝑎 ; 𝐾 ( 𝑟 & ƒ𝜑 & + '𝐶 $ & + 𝑟 & 𝑎 ; 𝐶 ( 𝑟 & ƒ𝜑̇& -𝑟 & 𝑎 ; 𝐾 ( 𝑟 ! 𝜃 ! -𝑟 & 𝑎 ; 𝐶 ( 𝑟 ! 𝜃 ̇! = 𝑇 ( (35)
Therefore, the stiffness and mass matrices for the reduced system are given by:

[𝑀] = ' 𝐼 & ' 0 0 𝐼 ! " " ; [𝐾] = ' 𝑟 & 𝑎 ; 𝐾 ( 𝑟 & -𝑟 & 𝑎 ; 𝐾 ( 𝑟 ! -𝑟 ! 𝑎 ; 𝐾 ( 𝑟 & 𝑟 ! 𝑎 ; 𝐾 ( 𝑟 ! " ; 𝑞 ⃑ = {𝜑 & 𝜃 ! } F (36)
Under the hypothesis of synchronous motion, i.e., 𝑞 ⃑ = ℎ(𝑡) • 𝑢 R⃑, Eqs. ( 33) and ( 36) can be expressed in a compact form as:

([𝐾] -𝜔 8 ; [𝑀]) • 𝑢 -⃑ 8 = 0 (37)
Where 𝑢 --⃑ 𝑛 represents the n-th eigenvector associated with the n-th eigenvalue 𝜔 , . By solving Eq.

(37), we can obtain four natural frequencies, one of which must be zero due to the rigid body motion.

Table 1 shows the non-zero natural frequencies for two different systems. Fig. 10 shows all three rotational mode shapes associated with the natural frequencies. Four nodes represent the rotational modes of motor (node no. 1), pinion (node no. 2), gear (node no. 3), and load (node no. 4) when the excitation frequency comes to the natural frequencies, see Table 2. The normalized mode shapes are represented in Table 2. For each mode, the maximum parameter is considered to normalize the mode.

Generic viscous damping

A classical method for considering a dissipation in a dynamical system is the Rayleigh proportional damping: the damping matrix C is a linear combination of the mass matrix M and the stiffness matrix

K, i.e.: [𝐶] = 𝛼[𝑀] + 𝛽[𝐾].
It should be kept in mind that this damping model is a kind of approximation of unknown dissipation mechanisms, while it does not change the modes of vibration (eigenvectors). It means that, for a generic system, i.e., non-proportionally damped systems, the damping matrix C cannot be diagonalized by the modes of the undamped system. Indeed, the presence of dissipation changes both the eigenvalues and the eigenvectors, which are in general complex. In order to generalize the modal analysis to generic damping, working in the phase space (𝑞̇, 𝑞) is required; in particular, in this paper, the Duncan approach is applied [START_REF] Foss | Co-Ordinates Which Uncouple the Equations of Motion of Damped Linear Dynamic Systems[END_REF]. what can be expressed in a general form governing equations is:

˜[𝑀]𝑞̇-[𝑀]𝑞̇= 0 [𝑀]𝑞̈+ [𝐶]𝑞̇+ [𝐾]𝑞 = f (38)
Where, 𝑞 ⃑ = {𝜑 & 𝜃 ! } F for the SDOF system and 𝑞 ⃑ = {𝜑 & 𝜃 ! 𝜑 ( 𝜃 % } F for the three-DOF system. Similar to the undamped systems for the free vibration problem, it can be written as:

[𝐴]{𝑦̇} + [𝐵]{𝑦} = 0 ( 39
)
(39) is the governing equation in the state space. Where:

[𝐴] = • 0 𝑀 𝑀 𝐶 ž ; {𝑦̇} = d 𝑞̈ 𝑞̇Ÿ ; [𝐵] = • -𝑀 0 0 𝐾 ž ; {𝑦} = d 𝑞̇ 𝑞 Ÿ ( 40
)
Using a general form of the solution, i.e., 𝑦(𝑡) = 𝑣 𝑒 8' , into Eq. ( 39), one obtains an eigenvalue problem; (𝐴𝜆 + 𝐵)𝑣 = 0. The characteristic equation is a polynomial of degree 2N: det(𝐴𝜆 + 𝐵) = 0.

It provides 2N complex conjugate eigenvalues and 2N complex conjugate eigenvectors. The introduction of damping changes the nature of the eigenvalue problem associated with the free vibration of a mechanical system: without damping one has N real eigenvalues and eigenvectors; with damping, one has 2N pairs of complex conjugate eigenvalues and eigenvectors. Table 3 represents the eigenvalues and eigenvectors for damped and undamped systems with three-DOF. As seen from ( 41), extracted damping ratios are represented:

ξ 6 = 0.045 ξ ; = 0.024 ξ : = 0.020 (41)
Eigenvalues discrepancy for the damped and undamped systems is not observable as considered damping for the system are small so their effects on the eigenvectors are not remarkable.

Numerical results

Nonlinear differential Eqs. ( 16)- [START_REF] Molaie | Spiral Bevel Gears: nonlinear dynamic model based on accurate static stiffness evaluation[END_REF] and Eq. ( 28) with time-varying parameters represent a nonsmooth and non-autonomous dynamical system with three and single DOFs, respectively after removing the rigid body motion. These equations are solved numerically through an implicit fifthorder Runge-Kutta scheme (RADAU) coded in FORTRAN language, the algorithm is enough stable and accurate, see e.g., Refs. [START_REF] Samani | Nonlinear vibration of the spiral bevel gear with a novel tooth surface modification method[END_REF], [START_REF] Molaie | Nonlinear vibration of crowned gear pairs considering the effect of Hertzian contact stiffness[END_REF], [START_REF] Motahar | Nonlinear vibration of the bevel gear with teeth profile modification[END_REF], [START_REF] Faggioni | Dynamic optimization of spur gears[END_REF]. In addition to the external constant torque (power transmitted), the time-varying mesh stiffness excites the system parametrically.

The time responses calculated by direct simulation are used for building amplitude-frequency and bifurcation diagrams of Poincaré maps, where the control parameter is the pinion rotation speed. For each simulation, the transient is separated from the steady-state part, i.e., the first 10000 periods are assumed to be transient, and beyond that, the influence of excitation frequency variation is considered expired; therefore, the first 10000 periods of simulation are never recorded and analyzed.

The Poincaré map diagrams are derived from the last 1000 periods of time responses. Each Poincaré map is obtained by sampling the time histories with the same period of meshing frequency. It must be pointed out that the governing equations are written in dimensionless form for using the natural frequencies for dimensionless time; in particular, 𝜔 , GHIJ is used for the SDOF, and the second natural frequency, 𝜔 , KL.0 is used for the three-DOF. All necessary data for dynamic simulation are given in Appendix A. and Appendix B.

SDOF vs three-DOF systems: effect of the low torsional shaft stiffness

In this section, the three-DOF system, see Fig. 2, is analyzed considering low torsional stiffness: in particular, when the shaft flexibility is extremely low, 𝐾 ' " = 𝐾 ' ! ≈ 0, Eqs. ( 1)-( 5) can be reduced to a single degree of freedom. The analyses are carried out on the linearized system and the nonlinear one.

For the linear analysis, natural frequencies and modes are investigated. For the nonlinear analysis, the forced response is studied.

The following stiffness is considered for the three-DOF model, 𝐾 ' " = 2.76 × 10 M N • m and 𝐾 ' ! = 8.27 × 10 M N • m, and zero for the SDOF model to obtain the natural frequencies of the systems (Table 1); the differences between meshing natural frequency and second natural frequency is 9%. In order to understand the relationship between the natural frequencies of and torsional shaft stiffness, the respectively. The accuracy of the system with a SDOF is analyzed in Refs. [START_REF] Samani | Nonlinear vibration of the spiral bevel gear with a novel tooth surface modification method[END_REF], [START_REF] Molaie | Spiral Bevel Gears Nonlinear Vibration Having Radial and Axial Misalignments Effects[END_REF]- [START_REF] Molaie | Spiral Bevel Gears: nonlinear dynamic model based on accurate static stiffness evaluation[END_REF], so by this comparison, we can assume that the three-DOF model is validated as well. Fig. 13 depicts the rootmean-square (RMS) of the system response vs the excitation frequency; simulations are carried out for both increasing and decreasing the excitation frequency (forward and backward). As the excitation frequency depends on the rotation speed of the input shaft for backward simulation, the rotation speed starts from the highest value to the lowest value; on the contrary, in the forward simulation, the speed starts from the lowest to the highest speed (from low to high excitation frequency). In agreement with the linear analysis, in these simulations, the damping ratio for the SDOF is 0.024, while for the three-DOF system, the damping ratios are ξ 5 = 0.045, ξ 0 = 0.024, and ξ V = 0.020. The backlash function is the main source of nonlinearity, as it introduces a discontinuity; the system may experience three different states: drive-side contact (desirable contact), tooth separation, and backside contact (coast-side contact). As it was explained in section 2, when the amplitude of 𝜆 (DTE) is between -1 and 1, tooth separation takes place. The discontinuity and presence of torque give rise to a softening-type resonance (vibration amplitude bent to the left).

Both the SDOF and the three-DOF models experience three different jumps due to different resonances: the primary resonance, the super-harmonic resonance, and the parametric resonance:

• Super-harmonic resonance, $%.' = 0.505 for the three-DOF system; see Table 4 -case 𝑪.

• Primary resonance,

N ! N " = 1:
Ø Forward: jump occurs at

N ! N "
$%.' = 0.859 for the three-DOF system, and

N ! N "
()*+ = 0.849 for the SDOF system; see Table 4 -case 𝑬.

Ø Backward: jump occurs at

N ! N "
$%.' = 1.132 for the three-DOF system, and

N ! N "
()*+ = 0.936 for the SDOF system; see Table 4 -case 𝑩.

• Sub-harmonic resonance,

N ! N " = 2:
Ø Forward: jump occurs at

N ! N "
$%.' = 1.907 for the three-DOF system, and

N ! N "
()*+ = 1.9 for the SDOF system; see Table 4 -case 𝑭.

Ø Backward: jump occurs at

N ! N "
$%.' = 2.087 for the three-DOF system, and

N ! N "
()*+ = 2.095 for the SDOF system; see Table 4 -case 𝑨.

Additionally, instability is observed in both the SDOF and the three-DOF systems at different regimes:

• The three-DOF system: • The SDOF system: forward simulation;

Ø Backward simulation; N ! N " $%.' ∈ [2.
N ! N " ()*+ ∈ [1.078, 1.082].
The tooth separation in the system represents itself as a jumping phenomenon: when the system is excited from low to high frequencies, the response grows until the point that any further increase in excitation frequency would cause a spontaneous jump in the amplitude of the dynamic transmission error, after the jump a further increase in the frequency causes a reduction of amplitude. The superharmonic peaks ( As soon as the nondimensional amplitude of the response falls below -1, the undesired backside contact occurs; this phenomenon takes place at different frequencies for the SDOF system. Once the backside contact is developed, increasing excitation frequency causes an increment of amplitude, indicating a hardening behavior. The response exhibits backside contacts at the end of the softening branch for

N ! N " = 6 ; , 6
\ C \ D
HIJK ∈ [0.937, 0.973 ] (see Table 4 -Case 𝐴 _ ). On the contrary, backside contact does not appear in the model of the system with three-DOF.

By analyzing the amplitude-frequency diagram, the main and remarkable difference between SDOF and three-DOF systems comes to the possibility of backside contact phenomena. In the current case, considering torsional shaft stiffness, i.e., the three-DOF model, prevents the system from experiencing undesirable backside contact. Now the bifurcation diagrams of the Poincaré maps are analyzed in order to complete the analysis.

The maps are built by sampling the time histories at the same frequency of excitation. When the frequency is changed, the initial condition is the final state of the system at the previous step. In order to eliminate the transient due to change of frequency, the first 10000 periods are not recorded. From the bifurcation diagram, the differences between the SDOF and three-DOF systems are highlighted in Fig. 15 (backward) and Fig. 16 (forward). Comparing the bifurcation and amplitude-frequency diagrams explain better the dynamic phenomena. Indeed, by considering amplitude-frequency diagrams, it seems that SDOF and three-DOF systems experience the same behavior, conversely there are remarkable differences in the bifurcation diagrams, see Fig. 15 and Fig. 16 (light blue islands).

Forward simulation: results shown in Fig. 15. First, both systems experience a periodic behavior, until the excitation frequency approaches the fundamental frequency, where the system experiences a jump phenomenon due to the primary resonance followed by a tooth separation; at

N ! N " = 0.86
(three-DOF) and 0.85 (SDOF). At the beginning of tooth separation, the SDOF system remains periodic, but the three-DOF system exhibits chaotic behavior. After that, both systems first experienced a chaotic response, then at

N ! N "
≈ 1.08 4-T responses take place, and again the SDOF and three-DOF systems become periodic at

N ! N "
= 1.13, and 1.15 respectively, see zone-I. By increasing the speed and approaching the sub-harmonic range, both systems move from periodic response to the 2-T responses for the SDOF system, and non-periodic response for the three-DOF system, see zone-II. Therefore, we can claim that using a SDOF system, some dynamics are not correctly captured. Backward simulation: Fig. 16 shows the bifurcation diagrams of the Poincaré maps vs the mesh frequency ratio in the case of backward speed variations. For each frequency, the final condition of each response is considered the initial condition of the subsequent frequency. The dynamic behavior initially (high frequency) exhibits a periodic response until the excitation frequency attains (

N ! N " = 2),
where instability due to parametric resonance occurred. The dynamic response of the SDOF system results in the attainment of a 2-T response; while an non-periodic response is evidenced from the three-DOF system at the onset ( N ! N " $%.' ∈ [2.087, 1.975]), see Fig. 16: zone-I. Subsequently, the dynamic response of the three-DOF system shifts to 2-T ( [1.328, 1.268]) and 4-T ( [1.414, 1.329], [1.130, 1.089]) responses, Fig. 16: zone-II. The phenomenon, visible in Fig. 16: zone-II, is called a boom-and-bust cycle. It is a phenomenon where a system or a process experiences a period of strong expansion or growth, often characterized by high levels of activity, followed by a sudden and severe contraction or decline. Indeed, During the boom phase, oscillators experience high levels of vibration. However, this growth is often unsustainable and eventually leads to a plummet to a low level of vibration [START_REF] Alligood | CHAOS: An Introduction to Dynamical Systems[END_REF]. Finally, it leads to an aperiodic behavior, Fig. 16: zone-III.

N ! N " $%.' ∈ [1.955, 1.415],
N ! N " $%.' ∈ [1.974, 1.956],
At the super-harmonic frequency regime, the three-DOF system experiences an unsteady behavior due to the second super-harmonic resonance, see Fig. 16: zone-IV.

Nonlinear dynamic scenario: complex analysis

From the analysis of Fig. 15 and Fig. 16, it can be proved that the three-DOF system experienced some phenomena that do not occur in the SDOF system although in both systems just rotational DOFs are considered. Based on the frequency at which the system is excited, we can observe a specific type of state or phenomena: periodic, period doubling route to chaos, quasiperiodic, chaos, trapping, sideband, and modulation. The trapping phenomenon refers to the phenomenon of a vibrating object becoming stuck or confined to a particular region of motion. The motion of the system is often complex and can involve multiple frequencies and modes of vibration. In some cases, the system may encounter non-linearities or energy dissipation mechanisms that cause it to become trapped in a particular region of the phase space. One way to understand trapping in vibrating systems is to consider the concept of resonance. Resonance occurs when the frequency of the external excitation matches one of the natural frequencies of the vibrating object. At resonance, the amplitude of the vibration can increase dramatically, which can lead to instability and even damage to the object. However, in some cases, the vibration may not be able to escape the region of the resonance due to damping or other factors. This can result in trapping, where the object remains stuck in a particular region of the phase space despite the presence of external excitation. For further information about trapping phenomena see Refs. [START_REF] Galias | Short periodic orbits for the Lorenz system[END_REF]- [START_REF] Lozi R | Coexisting Chaotic Attractors in Chua's Circuit[END_REF].

All Four cases, plotted in Fig. 17 and Fig. 18 with distinctive characteristics belong to states where the system response is non-periodic and the Poincaré map presents a complex shape. To illustrate the nonlinear dynamic behavior of the three-DOF system, amplitude-frequency, and bifurcation diagrams are shown in Fig. 17 and Fig. 18, for the backward and forward simulations; in order to improve the graphical impact, bifurcation, and amplitude-frequency diagrams are plotted with different vertical axes: non-dimensional amplitude and root-mean-square of the amplitude, respectively; the same horizontal axis is considered: frequency ratio

\ C \ D
. The challengeable zones are recognized and some cases at specified excitation frequencies are chosen for further investigation. There are several regions where the system experienced non-periodic behavior, the most remarkable of which are represented in subfigures: details-I, II, and III for the backward simulation, and details-I and II for the forward simulation. Among all ranges of frequency ratios in these three highlighted zones, four cases are chosen to analyze the behavior of the system in detail by means of time-histories, spectra, 3D-phase portrait, and Poincaré maps, see cases-I, II, III, and IV in Fig. 17 and Fig. 18.

By analyzing and investigating the results of Fig. 17 and Fig. 18, it can be claimed that the system represented the same behavior in both the forward and the backward simulations although with different amplitudes at different frequency regimes. Each case of Fig. 17 and Fig. 18, case-I, II, III, and IV, represent different dynamic scenarios. In the following parts, a comprehensive explanation is provided on each case: case-I, quasiperiodic behavior; case-II, trapped response; case-III, aperiodic response; and case-IV, a combination of aperiodic and modulation behavior [START_REF] Moon | Chaotic Vibrations: An Introduction for Applied Scientists and Engineers[END_REF], [START_REF] Masoumi | Symmetry breaking and chaos-induced imbalance in planetary gears[END_REF].

A) Case-I: quasiperiodic behavior

Quasiperiodic and modulated behaviors are observed in the case-I of Fig. 17 and Fig. 18, backward and forward simulations respectively, since we can see a closed curve in the Poincaré map and sideband in the spectrum. The spectrum shows a series of peaks at discrete frequencies, indicating the presence of periodic behavior with slight fluctuations over time.

The strength of periodic behavior is an important key point. If the periodic behavior is strong and there is a clear and prominent repetition of a particular frequency, the spectrum will exhibit pronounced peaks at those frequencies. Indeed, narrow and distinct peaks indicate that the system has a well-defined periodicity. However, if the periodic behavior is weak and the repetition is not as clear or prominent, the spectrum will not exhibit pronounced peaks.

Instead, the spectrum will appear more like a continuous band with some level of variability, indicating that the system has a weaker or more complex periodicity.

By investigating the spectrum, some peaks may be observed with a constant frequency shift, 𝛽 ; , which represented a beat frequency. Therefore, the approximate period of the torus cycle can be defined as: 𝑇 torus ≈ 2π/𝛽 ; . The distribution of the remaining incommensurate frequency peaks, in the spectrum is governed by the following formula [START_REF] Borkowski | FFT Bifurcation Analysis of Routes to Chaos via Quasiperiodic Solutions[END_REF]:

𝛺 8 = 𝛺 2 + 𝑛 • 𝛽 ; ( 42
)
Where 𝛺 𝑛 is known also as the second disproportionate frequency, and 𝑛 is the number of tori frequency; 𝑛D-torus or 𝑛-frequency quasiperiodic solution. Besides, the third disproportionate frequency can be described as:

𝛺 8,( = 𝛺 8 + 𝑚 • 𝛽 6 ( 43
)
Where 𝑛 and 𝑚 are the number of frequency peaks. One of the applications of the FFT spectrum diagram is to help us to understand more about modulation in the dynamic behavior of the system by observing the sideband in the plotted results. There are various types of modulation, including amplitude modulation (AM), frequency modulation (FM), phase modulation (PM), and pulse width modulation (PWM). Each type of modulation is used to convey information in different ways, depending on the specific requirements of the application. The FFT spectrum is a useful tool that helps us to recognize types of modulation, for instance, Fig. 19 represents an amplitude modulation for case-I for backward simulation.

Let us show this parameter for the case-I of Fig. 17

𝛽 ] N 𝛽 ] 6 𝛽 ] O 𝛽 ] 4 𝛽 ] P 𝛽 ] P 𝛽 ] N 𝛽 ] 4 𝛽 ] 6 𝛽 ] O 𝑇 6 = 0.001772 s 𝝀 ' 𝑭𝑭𝑻 𝒔𝒑𝒆𝒄𝒕𝒓𝒖𝒎

𝒕 [𝒔] B) Case-II: trapping phenomenon

Fourier spectrum analyses associated with cases-II of Fig. 17 and Fig. 18 are shown in Fig. 22 for the backward and forward simulations respectively, where

\ C \ D
= 0.98. We can define the excitation frequency ratio,

\ C \ D
, as the control parameter of the three-DOF system. For values of the control parameter less than the critical value, the dynamical system has an attracting limit cycle. Thus, the system oscillates in a regular fashion and is stable. As

\ C \ D
slightly exceeds the threshold value

\ C \ D
, the system response appears to be regular and closely resembles the oscillatory behavior for the

\ C \ D
less than the critical excitation frequency ratio (i.e., where the system has a periodic behavior). With increasing

\ C \ D
, the transient phase between two periodic and aperiodic responses passed and it becomes more and more difficult to recognize any regularity in the system response. As

\ C \ D
is increased further, eventually the response becomes fully irregular (chaotic). Therefore, case-II belongs to the zone where the system experiences a route to chaos; see Fig. 21, blue dots.

By observing the 2D and 3D-Bifurcation diagram (see Fig. 17 and Fig. 21), during the transient phase, 6 different branches (red dots) are recorded. Besides, the Poincaré map represents six trapped regions and the Fourier spectrum shows a weak modulation with a complex periodic behavior. Indeed, the chaotic attractors become a chaotic sea with islands representing semiquasiperiodic behavior and the semi-quasiperiodic islands surrounded by the chaotic orbit. As shown in Fig. 22, due to the weak modulation, there is a weak regularity in the time-history response of the system although the periodicity is not strong enough to be followed all over time; after a while, some irregularity is observed in the system behavior.

Ω 3 β ] N 2 • β ] N 3 • β ] N 𝝀 ' 𝑭𝑭𝑻 𝒔𝒑𝒆𝒄𝒕𝒓𝒖𝒎 𝑭𝑭𝑻 𝒔𝒑𝒆𝒄𝒕𝒓𝒖𝒎

𝒕 [𝒔] 𝒕 [𝒔]

A) Backward B) Forward 𝝀 '

Ω 3 β ] N 2 • β ] N 3 • β ] N 𝛽 Š 0 = 0.167 𝑇 N ≈ 0.00214 𝑠 𝑇 3 𝑇 3

C) Case-III: aperiodic response:

By analyzing the results of case-III in Fig. 17 Case-III belongs to the zone where there are no remarkable peaks at the FFT spectrum except at frequencies that are simple multiples of the system's meshing frequency. The results of case-III in Fig. 17 and Fig. 18 are obtained at excitation frequency before the dynamic response evolved into four cloudy branches; for this reason, the Poincarè map clearly shows three branches instead of four branches. These irregularities in the response lead to the spreading out of the energy of system over a wide range of frequency. It can be expressed that there is a relation between results that are shown in Poincarè map, where there are three branches for the case-III for instance, and Fourier spectrum, where an increase of the energy (concentration energy) in the spectrum at 1/3 and 2/3 of the meshing frequency is observed.

Therefore, by considering the 𝛽 Š 0 ≈ 0.333, a weak regularity is found in the time-history response of the system, see Fig. 24.

𝛽 ] N 𝛽 ] N 𝛽 ] N 𝛽 ] N 𝛽 ] N 𝛽 ] N 𝛽 ] N 𝛽 ] N 𝝀 ' 𝑭𝑭𝑻 𝒔𝒑𝒆𝒄𝒕𝒓𝒖𝒎 𝒕 [𝒔] 𝒕 [𝒔] 𝝀 ' 𝑭𝑭𝑻 𝒔𝒑𝒆𝒄𝒕𝒓𝒖𝒎 𝑇 3 𝑇 3

D) Case-IV: combination of irregularity and period-doubling behavior

Considering the frequency range

\ C \ D EF.G ∈ [0.85, 1.1]
, a coincidence of 2T and aperiodic behavior are observed in the bifurcation diagram, see Fig. 17 and Fig. 18 for the backward and forward simulations respectively. Case-IV in Fig. 17 and Fig. 18, where

\ C \ D EF.G = 2.02, is a
sample that is chosen among the specific range of frequencies to analyze the behavior of the system in detail. From Fourier spectrum, 2T responses are conducted through the clear spike at

\ \ C = 0.5.
The broadband in the FFT spectrum with two remarkable peaks appeared, and also in the Poincaré map we can see that plotted points are separated into two distinct parts.

However, due to the strong irregularity and the weak regularity, it cannot be observable by analyzing the time-history responses. Indeed, the pattern that we expect to see due to 2Tresponses does not regularly occur.

Largest Lyapunov exponent and Correlation dimension estimations:

Investigating complex nonlinear dynamical systems, computational techniques are utilized such as nonlinear time series analysis. Chaos arises from the exponential growth of infinitesimal perturbations, together with global folding mechanisms to guarantee boundedness of the solutions. This exponential instability is characterized by the spectrum of Lyapunov exponents [START_REF] Eckmann | Ergodic theory of chaos and strange attractors[END_REF]. If one assumes a local decomposition of the phase space into directions with different stretching or contraction rates, the spectrum of exponents represents the average of these rates across the entire invariant set, and thus such spectrum includes an exponent for each spatial direction [START_REF] Stoop | Calculation of Lyapunov exponents avoiding spurious elements[END_REF].

In terms of chaos theory, attractors, which are geometrical objects formed by the trajectories of the system, are characterized by fractal dimensions. The fractal dimension of attractors arising from regular deterministic systems, such as limit cycles or tori, is equal to their topological dimension.

However, in the case of attractors derived from chaotic systems, called strange attractors, their fractal structure is usually characterized by a non-integer fractal dimension. Among various measures of the fractal dimension the most common is the correlation dimension due to its computational simplicity [START_REF] Mccue | Use of Lyapunov Exponents to Predict Chaotic Vessel Motions[END_REF].

Before estimating the largest Lyapunov exponent and correlation dimension, the phase space reconstruction approach must be defined. One of the main techniques used to reconstruct a phase space from a time series is the Method of Delays, which is based on the evaluation of the embedding dimension and the delay time. It has become quite familiar in the analysis of "observe time series" from nonlinear systems to make a time-delay reconstruction of a phase space (Method of Delay approach) to illustrate the dynamics [START_REF] Kennel | Determining embedding dimension for phase-space reconstruction using a geometrical construction[END_REF]. This point must be considered that the shape and size of the neighborhoods are adjusted based on the Lyapunov exponents. It considers that embedding in an 𝑚-dimensional space where 𝑚 < 𝑚 3 . As a result of this projection, the original topological structure is no longer maintained. It means that points are mapped into neighbors' points which they would not belong in higher dimensions. These points are called false neighbors. If now the dynamics is applied, these false neighbors are usually not mapped into the vicinity of their original neighbors but rather to different locations, causing the average diameter to increase significantly [START_REF] Hegger | Practical implementation of nonlinear time series methods: The TISEAN package[END_REF], [START_REF] Rosenstein | A practical method for calculating largest Lyapunov exponents from small data sets[END_REF], [START_REF] Deshmukh | Using scaling-region distributions to select embedding parameters[END_REF].

Largest Lyapunov exponent:

Lyapunov exponents, which associated with the fast divergence or convergence of neighboring trajectories in phase space, provide us with this ability to characterize the dynamical behavior of a system, both qualitatively and quantitatively [START_REF] Wolf | Determining Lyapunov exponents from a time series[END_REF]:

• The signs of the Lyapunov exponents (positive or negative) provide a qualitative understanding of the dynamic response of a system. As an example, one-dimensional maps are characterized by a single Lyapunov exponent which is positive for chaos, zero for a marginally stable orbit, and negative for a periodic orbit.

• The magnitudes of the Lyapunov exponents quantify an attractor's dynamics in terms of information theory. The exponents measure the rate at which system processes generate or dissipate information.

Lyapunov exponents are convenient for categorizing steady-state behavior. Attractors are classified as follows [START_REF] Parker | Chaos: A Tutorial for Engineers[END_REF], [START_REF] Pincus | Approximate entropy as a measure of system complexity[END_REF]:

• Positive Lyapunov Exponent (𝜆 g > 0) indicates a chaotic system: It means that nearby trajectories in phase space diverge from each other exponentially over time. It is a characteristic of chaotic behavior, as a small perturbation in initial conditions result in a bounded aperiodic response over time.

• Zero Lyapunov Exponent (𝜆 g = 0) indicates periodicity in system: It means that nearby trajectories neither converge nor diverge but stay close together over time. Systems with zero Lyapunov exponents are often associated with periodic or quasiperiodic behavior.

For instance, for a stable limit cycle, 𝜆 5 = 0 and 𝜆 g < 0 for i = 2, … , 𝑛, and for a stable torus, 𝜆 5 = 𝜆 0 = 0 and 𝜆 g < 0 for i = 3, … , 𝑛.

• Negative Lyapunov Exponent (𝜆 g < 0) indicates the existence of a stable equilibrium points: It means that nearby trajectories in phase space converge toward each other exponentially. Negative exponents are associated with stable or asymptotically converging behavior.

The maximal Lyapunov exponent can be determined without the explicit construction of a model for the time series. [START_REF] Rosenstein | A practical method for calculating largest Lyapunov exponents from small data sets[END_REF], [START_REF] Kantz | A robust method to estimate the maximal Lyapunov exponent of a time series[END_REF]. 𝑥(𝑡) is the time evolution of initial condition, 𝑥(0), in an appropriate state space. Then the maximal Lyapunov exponent is found with "probability one", which is:

𝜆 1UV = lim W→Y lim Z→2 1 t ln y |x(t) -x Z (t)| 𝜖 z ( 44
)
for all difference vectors 𝜖 = x(0) -x h (0).

The first algorithm to compute Lyapunov exponents for a time series was introduced in 1985 by Wolf et al. [START_REF] Wolf | Determining Lyapunov exponents from a time series[END_REF]. In delay coordinates of appropriate dimension, one seeks a point within the time series that is closest to its initial point. This is considered as the beginning of a neighboring trajectory, determined by the consecutive delay vectors. Therefore, one calculates the growth or increase in the distance between these two trajectories over time. When the distance exceeds some threshold, a new neighboring trajectory is sought with the closest distance and similar direction as the old one. This is the direction of the local eigenvector associated to 𝜆 max . The logarithms of the stretching factors of the difference vectors are averaged over time to calculate the maximal Lyapunov exponent [START_REF] Kantz | A robust method to estimate the maximal Lyapunov exponent of a time series[END_REF].

Similar to Wolf et al. [START_REF] Wolf | Determining Lyapunov exponents from a time series[END_REF], Kantz et al. [START_REF] Hegger | Practical implementation of nonlinear time series methods: The TISEAN package[END_REF], [START_REF] Kantz | A robust method to estimate the maximal Lyapunov exponent of a time series[END_REF] used of the fact that the distance between two trajectories typically increases with a rate given by the maximal Lyapunov exponent.

Correlation dimension

A dynamical system with 𝑁 degrees of freedom may choose to evolve on a lower-dimensional manifold, so that only a fraction of the total number of degrees of freedom remains active. In such scenarios, it becomes valuable to determine the number of active degrees of freedom, and it is obvious that this information can be deduced from the dimension of the attractor associated with the system [START_REF] Galka | Topics in Nonlinear Time Series Analysis[END_REF]. If it looks like an open subset of ℝ , in a neighborhood of every point, an attractor could be defined to be 𝑛-dimensional, i.e., it is diffeomorphic to an open subset of ℝ , . This is the approach to define the dimension of a manifold in differential topology. For instance, a limit cycle is one dimensional since it looks, locally, like an interval. A torus is two-dimensional since, locally, it resembles an open subset of ℝ 0 . An equilibrium point is considered to have zero dimension. The neighborhood of any point of a strange attractor, however, has a fine structure and does not resemble any Euclidean space. Hence, strange attractors are not manifolds and do not have integer dimensions.

This point worth to mention that an attractor of a Poincare map with dimension d corresponds to an attractor of the underlying flow with dimension d + 1. For example, if the attractor of the Poincare map is a closed curve (dimension, d =1), the attractor of the flow is a torus (dimension 2).

There are several ways to generalize dimension to the fractional case [START_REF] Parker | Chaos: A Tutorial for Engineers[END_REF] such as capacity, information dimension, Lyapunov dimension and correlation dimension. In this study, the correlation dimension, which is introduced in 1983 by Grassberger et al. [START_REF] Grassberger | Characterization of Strange Attractors[END_REF], is chosen to characterize the dynamic behavior of system in different working conditions.

Validation: periodic case,

𝝎 𝒎 𝝎 𝒏 𝐍𝐨.𝟐 = 𝟏. 𝟐

As mentioned before, the parameters used in chaotic recognition, are the largest Lyapunov exponents (𝜆 lYm ) and correlation dimensions. In the present study, to conduct the nonlinear time series analysis, "nonlinear time series analysis package (TISEAN-V3.0.1)" [START_REF] Hegger | Practical implementation of nonlinear time series methods: The TISEAN package[END_REF], and a script written in MATLAB [START_REF] Mccue | Use of Lyapunov Exponents to Predict Chaotic Vessel Motions[END_REF], [START_REF] Rosenstein | A practical method for calculating largest Lyapunov exponents from small data sets[END_REF] are used. In appendix C, the script to carry out the calculations by TISEAN and MATLAB is defined. The TISEAN package was employed to validate the accuracy of the data extracted from the MATLAB algorithm in the case of

\ C \ D EF.G = 1.
2, (backward simulation, Fig. 26-a) where the behavior of system has already recognized, and then for other cases the MATLAB code written is used to calculate the 𝜆 lYm and correlation dimension. As it is clear from Fig. 26-a and Fig. 26-b, where time history and Poincaré map of the system is represented, the system experiences a periodic response. As explained, the delay time is determined as the first local minimum of the mutual information function (Fig. 27-a), and embedding dimension is defined by using a 10% threshold [START_REF] Deshmukh | Using scaling-region distributions to select embedding parameters[END_REF] of fraction of false neighbors (Fig. 27-b). Extracted results from Fig. 28 are listed in Table 5. The estimated correlation dimension is equal to 1, and the calculated λ lYm for both approaches is weakly negative, i.e., the system exhibits a stable

Estimation of Largest Lyapunov exponent and correlation dimension

The calculations have been performed using the results extracted during the backward simulation for all the four investigated cases, where a different excitation frequency was considered as shown in Fig. 17. The parameters required for the calculation are given in Table 6, such as delay time, embedding dimension, total considered samples, and sampling intervals. By analyzing the extracted results in Table 6 and Fig. 29, the following key points is concluded:

• For case-I, where the behavior of the system is quasiperiodic, a flat plot for the Average Log Divergence function is expected [START_REF] Rosenstein | A practical method for calculating largest Lyapunov exponents from small data sets[END_REF] as we can see in Fig. 29

Torsional shaft stiffness effects on the dynamic response

The effect of shaft stiffness on dynamic behavior is critical from a design perspective. This analysis offered guidelines for the design of gear systems in order to choose suitable shafts and bearings.

Indeed, it would be interesting if we can change the dynamic response of the system from unstable behavior to periodic response. In this section, the dynamic behavior of the system is evaluated by varying the torsional shaft stiffness. In order to perform this analysis, it is first necessary to determine the excitation frequency at which the system must operate. Therefore, a fixed-initial condition simulation is carried out. Then at a specific frequency ratio, the new simulation is conducted to extract the bifurcation diagram vs the torsional shaft stiffness. This diagram can help to determine the torsional shaft stiffness value at which the system transitions to either an aperiodic or periodic response. In other words, we can prevent the system from chaotic behavior by considering the suitable stiffness for the shaft, i.e., considering different shapes or materials for the shaft. ) and gear (𝐾 ' ' ! ) as the control parameter of the system, see Fig. 30, Analyses 1-4.

Fig. 30-Analysis-1 (where the 𝐾 ' ' " is the control parameter) and Fig. 30-analysis-2 (where the 𝐾 ' ' ! is the control parameter) are associated with the case-I where the system experiences an aperiodic behavior that can be chaos. By analyzing the bifurcation diagrams, it is conducted that by increasing the 𝐾 ' ' ! and reaching the highlighted zone with green color in Fig. The same analysis can be carried out for the case-II, see Fig. 30-Analysis-3 (where the 𝐾 ' ' " is the control parameter) and Fig. 30-analysis-4 (where the 𝐾 ' ' ! is the control parameter). Therefore, to skip the period doubling responses, the behavior that the system experienced at case-II, there are two different ways: decrease 𝐾 ' ' " or increase 𝐾 ' ' ! . Indeed, based on the analyses which are done for the analysis-3 and analysis-4, it can be concluded that to have a periodic response for the system, 𝐾 ' ' " must be less than 0.367, or 𝐾 ' ' ! must be higher than 0.146.

Conclusion

This paper described the nonlinear dynamic behavior of a spiral bevel gear pair by considering the effects of torsional shaft stiffness. The dynamic simulation of the system was conducted considering the possibility of backside contact. This consideration led to extracting the mesh stiffness in forward and reverse motions. A modal analysis of the linear system was conducted, and validations were carried out to ensure the accuracy of the present model. The dynamical model, which is governed by a set of non-smooth and non-autonomous equations of motion, was numerically solved by means of an implicit Runge-Kutta algorithm. Nonlinear time series analysis was carried out to investigate more the dynamic behavior of different cases where we observe different complexity. To understand the effect of the shaft stiffness on the system's dynamics, two different models were represented: three-DOF and SDOF. Dynamic simulations were carried out both increasing and decreasing the excitation frequency. The main outcomes of this investigation are listed below:

• At first glance, it might be concluded that the dynamic behavior of SDOF and three-DOF systems are identical; however, by analyzing the dynamic response and delving into the possibility of phenomena occurrence, it is proved that by decreasing the DOF of the system, some phenomena might be lost.

• Different phenomena appeared in the system by considering the shaft stiffness, such as the trapping phenomenon and boom-and-bust cycle. By estimating the largest Lyapunov exponent and the correlation dimension, the behavior of the phenomena is evaluated.

• There were some situations where the SDOF system represents a periodic or 2T-periodic responses, whereas the behavior of the three-DOF model was markedly different in the same scenarios, even an aperiodic or quasiperiodic responses might be recorded.

Consequently, designing the system based on the results of the SDOF model might not be suitable for manufacturing a gear set.

• The results provide a foundation for designing parameters and controlling the dynamic characteristics of spiral bevel gear drive systems.

1. Dimensional parameters:

Fig. 1 .

 1 Fig. 1. The dynamic model of the gear system with torsional-DOF.

Fig. 2 .

 2 Fig. 2. Torsional dynamic model of the gear system with three-DOF.

Fig. 3 .

 3 Fig. 3. Geometric transmission error: lack of material along the line of action

Fig. 4 . 23 𝑁" 9 4+6

 4239 Fig. 4. Equivalent gear model and backlash function. The meshing stiffness of the gear pair is a time-varying periodic function, with the mesh frequency, 𝜔 ( = 01 23 𝑁 " 𝑛 4 , where 𝑛 4 is the input shaft speed and 𝑁 " is the number of pinion teeth;the equivalent meshing stiffness Fourier series is given by:

Fig. 5 .

 5 Fig. 5. Dynamic model of gear system with SDOF.

Fig. 6 .Fig. 7 .Fig. 8 .

 678 Fig. 6. Transmission error (gear rotation) -LTCA: torque 537.13 N⋅m: a) forward motion, b) reverse motion; ▬ Fitted curve, • FEM results. a) forward motion b) reverse motion

Fig. 9 .

 9 Fig. 9. Mesh stiffness diagram: a) forward motion, b) reverse motion; ▬ Fitted curve by Fourier series, • FEM results.

Fig. 10 .

 10 Fig. 10. Normilized mode shapes for all three natural frequencies: First natural frequency (•), second natural frequency (•), and third natural frequency(•).

a)

  First mode: 𝜔 /012.4 = 2704 rad s N b) Second mode: 𝜔 /012.5 = 17957 rad s N c) Third mode: 𝜔 /012.6 = 35140 rad s N Node No. 1 represents rotation of Motor. Node No. 2 represents rotation of Pinion. Node No. 3 represents rotation of Gear. Node No. 4 represents rotation of Load.

U

  Fig. 11. Non-zero natural frequency of the system with three-DOF. If a sufficiently low value for the 𝐾 ' " and 𝐾 ' ! is considered (𝐾 ' " = 27.6 N • m and 𝐾 ' ! = 82.7 N • m), the mesh natural frequency for both systems gradually approaches each other, 𝜔 ,OKL.V NWXQQOHIJ = 19595 XYZ [ ≈ 𝜔 , GHIJ , while the first and second natural frequencies for the three-DOF system tend toward zero. Besides, a linear simulation is carried out to extract the DTE time-history of both systems,

Fig. 12 . 6 '

 126 Fig. 12. Extracted time-history response from linear solutions with the low torsional shaft stiffness; system with SDOF (•), and three-DOF (•). Now the fully nonlinear models are considered for investigating the differences between the SDOF and three-DOF models, while a low value for the 𝐾 ' " and 𝐾 ' ! is considered, 27.6 N • m and 82.7 N • m

Fig. 13 . 5 M 3 . 2 .Fig. 14 .

 1353214 Fig. 13. RMS response extracted from nonlinear dynamic simulations with the low torsional shaft stiffness; system with SDOF (•), and three-DOF (•).As the comparison revealed, both the SDOF and the three-DOF models experienced backside

  ' = 0.49 for the three-DOF system, and N ! N " ()*+ = 0.489 for the SDOF system; seeTable 4 -case 𝑫.

: , 6 O= 1 = 2 ,

 612 , …) are resonances caused by the super-harmonics of the meshing stiffness; the one at the fundamental frequency close to N ! N " is the primary resonance, which is due to the fundamental harmonic. Close to N ! N " instability due to the principal parametric resonance (Mathieu type) is present, and the resonance branch shows a softening character due to loosing of contact.

Fig. 15 .

 15 Fig. 15. Bifurcation diagram extracted from forward simulation; comparison between the systems with SDOF (•), and three-DOF (•).

Fig. 16 .

 16 Fig. 16. Bifurcation diagram extracted from backward simulation; comparison between the systems with one-DOF (•), and three-DOF (•).

  in Fig. 19 (backward simulation) and the case-I of Fig. 18 in Fig. 20 (forward simulation), where

Fig. 19 .Fig. 20 .

 1920 Fig. 19. FFT spectrum of case-I for backward simulation

Fig. 21 .Fig. 22 .

 2122 Fig. 21. 3-D bifurcation diagram for the backward simulation, N ! N " $%.' ∈ [0.88, 1.1]

= 1 ,

 1 and Fig.18, for backward and forward simulations, the system experiences an aperiodic behavior. The bifurcation diagram can provide a clearer exhibition of the behavior of the system, particularly the 3D-diagram, see Fig.23where the simulation is carried out for the backward simulation. As shown in Fig.23, there is a zone with the cloud of points (Fig.23, blue dots) which evolves from having three branches into four four cloudy branches of the dynamic responses of the system turn into 4T-periodic responses (Fig.23, pink dots).

Fig. 23 .Fig. 24 .

 2324 Fig. 23. 3-D bifurcation diagram for the backward simulation, N ! N " $%.' ∈ [0.9, 1]

Fig. 25 . 3 β 4 =

 2534 Fig. 25. FFT spectrum of case-I for backward simulation

Fig. 26 .

 26 Fig. 26. Dynamical response of the three-DOF model at N ! N " $%.' = 1.2: a) time history, and b) Poincaré map

Fig. 27 . 2 a

 272 Fig. 27. Extracting a) delay time and b) embedding dimension for the periodic case in backward simulation at N ! N " $%.' = 1.2

Fig. 28 .

 28 Fig. 28. Extracted largest Lyapunov exponent and correlation dimension from simulation done by TISEAN and MATLAB Fig. 28-c allows for the estimation of the correlation dimension. Considering a neighborhood

  steady state response. By comparing the extracted results between two techniques, TISEAN and TISEAN MATLAB a) average log divergence vs the expansion step d) Slope of the Correlation integral versus the neighborhood size c) Correlation integral versus the neighborhood size a) average log divergence vs the expansion step 51 MATLAB code, 1% difference in the calculated correlation dimension and 2% difference in the calculated λ lYm is observed.

Fig. 29

 29 Fig. 29-a2 Fig. 29-a1

Fig. 29

 29 Fig. 29-b2Fig.29-b1

Fig. 29

 29 Fig. 29-c2Fig.29-c1

Fig. 29

 29 Fig. 29-d2Fig.29-d1

Fig. 29 .

 29 Fig. 29. largest Lyapunov exponent and correlation dimension for four cases extracted from backward simulation: 1) largest Lyapunov exponent diagram: x-axis: expansion Step; y-axis: Average Log Divergence 2) correlation dimension diagram: x-axis: neighborhood Radius; y-axis: Correlation Integral

Fig. 30 -

 30 Fig. 30-a shows the dynamic response of the three-DOF system with the values which are defined in Appendix B., where 𝐾 ' ' ! = 0.0227 and 𝐾 ' ' " = 0.421. To delve into the effect of torsional shaft stiffness, two different parts of the excitation frequency range are selected, see Fig. 30-b. Now, there

Table 1 .

 1 Non-zero natural frequencies, [rad/s]. System with three-DOF (+1 rigid body mode)

	𝜔 8 GH.6 2704 System with SDOF (+1 rigid body mode) 𝜔 8 GH.; GH.: (meshing natural frequency) 𝜔 8 17957 35140 JKLM 𝜔 8 19592

Table 2 .

 2 Normalized mode shapes of three-DOF system.

	No. 1 2 3	Natural frequency [ rad 𝑠 ] 2704 17957 35140	motor 1 1 -1	Mode shape pinion gear 0.993 0.967 0.680 0.140 0.222 -0.001	load 0 0 0

Table 3 .

 3 extracted eigenvalues and eigenvectors for the three-DOF system.

		Damped system	Undamped system		Damped system	Undamped system
	Eigenvalues [rad/s]	0 Eigenvectors	0	Eigenvalues [rad/s]	55.1 + 2703.7i Eigenvectors:	2704i
	Motor	1	1	Motor	0.692 + 4.0 × 10 !" i	0.692 + 0i
	Pinion	1	1	Pinion	0.687 + 4.2 × 10 !" i	0.687 + 0i
	Gear	1	1	Gear	0.669 + 4.7 × 10 !" i	0.669 + 0i
	Load	1	1	Load	0	0
		Damped system	Undamped system		Damped system	Undamped system
	Eigenvalues [rad/s]	431.1 + 17952.3i Eigenvectors:	17957i	Eigenvalues [rad/s]	1571.3 + 35105.2i Eigenvectors:	35140i
	Motor	0.826 -9.7 × 10 !# i	0.826 + 0i	Motor	-0.976 -0.6 × 10 !# i	-0.976 + 0i
	Pinion	0.562 -0.2 × 10 !# i	0.562 + 0i	Pinion	0.217 -2.5 × 10 !# i	0.217 + 0i
	Gear	-0.116 + 7.5 × 10 !# i	-0.116 + 0i	Gear	-0.01 + 0.7 × 10 !# i	-0.001 + 0i
	Load	0	0	Load	0	0

Table 4 .

 4 The properties of starting and ending points of separation and backside contact for the forward

		motion
		SDOF	Three-DOF
	starting point	ending point	starting point	ending point

Table 5 .

 5 extracted data from nonlinear time series analysis for the periodic case-backward simulation

	Case study	Correlation dimension	Largest Lyapunov exponent	Delay time	Embedding dimension	Total samples	Sampling interval
	MATLAB TISEAN	1.01 1	-4.2094e-07 -4.1199e-07	6	3	39270	36

Table 6 .

 6 extracted data from nonlinear time series analysis for all four cases of backward simulation

	Case study		Correlation dimension	Largest Lyapunov exponent	Delay time	Embedding dimension	Total samples	Sampling interval
	Case I: 𝜔 𝜔 ! = 0.64	value Figure's number	2.08	-0.00013	9	3	587547	43

Table 9 .

 9 Values of considered parameters

	parameter	unit	value	parameter	unit	value
	𝑎 𝐶 ' " 𝐶 ( 𝑟 " 𝐼 " n 𝐼 ( 𝐾 ' " 𝐾 ( = 𝑘 3 𝑇 ( 𝑠 𝛼 c𝑟 & 𝑎d 𝑚 -. = b 5 𝐼 & Z + c𝑟 % 𝑎d 5 𝐼 % [ f 04	-N • s/m N • s/m mm kg • m 0 kg • m 0 N • m rad › N/m N • m -deg kg	0.76975 0.07 380 25.6 1.7 × 10 OM 2.7 × 10 Oe 2.76 × 10 M 1.4478 × 10 p 179.04 5 20 0.377182	𝑏 𝐶 ' ! G 𝑟 ! 𝐼 ! o 𝐼 # 𝐾 ' ! 𝑁 5 𝑇 # 𝜏 (qn 𝛽 𝜔 8A9ab3 = { 𝑘 2 𝑚 => s	mm N • s/m GPa mm kg • m 0 kg • m 0 N • m rad › -N • m MPa deg rad s ›	0.35 1.25 77 76.9 9.5 × 10 OV 2.84 8.27 × 10 M 21 537.12 50 35 19592
		Table 10. definition of considered parameters
	Half of the gear backlash Material damping of the pinion shaft Mesh damping Rotational inertia of pinion Rotational inertia of motor Torsional shaft stiffness of the pinion mean value of mesh stiffness mean radius at meshing point of the pinion Constant input torque Number of harmonics Normal pressure angle	𝑏 𝐶 ' " 𝐶 ( 𝐼 " n 𝐼 ( 𝐾 ' " 𝐾 ( = 𝑘 3 mean radius at meshing point of the gear Material damping of the gear shaft Shear Modulus Rotational inertia of gear Rotational inertia of load Torsional shaft stiffness of the gear Number of samples for UTCA and LTCA 𝑟 " Constant output torque 𝑇 ( maximum torsional stress 𝑠 The spiral angle 𝛼	! 𝐶 ' G o 𝐼 ! 𝐼 # ! 𝐾 ' 𝑁 5 𝑟 ! 𝑇 # 𝜏 (qn 𝛽

Table 11 .

 11 GTE and mesh stiffness values, for both the forward and the reverse motions

	Normalize	Time	Forward motion	Reverse motion
	step	step	GTE [µm]	MS [N/m]	GTE [µm]	MS [N/m]
	0	0.0025	57.4199	167077194	57.9535	65259864
	0.05	0.0050	66.7241	161191031	68.2185	64663942
	0.10	0.0075	76.1594	154213691	78.7352	63708319
	0.15	0.0100	85.7161	146394176	089.5220	62397179
	0.20	0.0125	95.1830	138138892	100.3302	60805706
	0.25	0.0150	102.8398	130738422	108.8252	59282854
	0.30	0.0175	108.3649	124235413	114.4459	57890915
	0.35	0.0200	111.6688	119726715	117.3800	56589953
	0.40	0.0225	112.7251	117679839	117.4634	55386935
	0.45	0.0250	111.5323	117253906	114.7359	54497509
	0.50	0.0275	108.0861	118398290	109.3410	54508453
	0.55	0.0300	102.3880	121170911	101.2252	55045242
	0.60	0.0325	94.4110	125810846	90.3533	56113207
	0.65	0.0350	84.1505	132729405	76.7455	57757463
	0.70	0.0375	71.5909	142560957	60.4187	60055301
	0.75	0.0400	56.7153	156481910	41.4372	63134970
	0.80	0.0425	39.5033	176556405	19.7876	67214271
	0.85	0.0450	30.3975	190867745	28.4547	66492226
	0.90	0.0475	39.2444	181912241	38.0885	65800291
	0.95	0.0500	48.2521	173705334	47.9154	65496346
	1	0.0525	57.4199	167077194	57.9535	65259864

Table 13 .

 13 Fourier series coefficients of GTE and MS for both the forward and the reverse motions

	𝑒̅ 2 𝐹 8.013 × 10 !&	Forward motion 𝑒𝑎 uuu 6 𝐹 𝑒𝑏 uuu 6 𝐹 -3.413 1.436 × 10 !& × 10 !&	𝑒𝑎 uuu ; 𝐹 -1.924 × 10 !"	𝑒̅ 𝑚 𝐹 (𝑡) = 𝑒𝑎 uuu : 𝑒 2 𝐹 𝑏 ³ =̅ . 𝐹 𝐹 1.852 × 10 !" 𝑒𝑏 uuu ; 𝐹 6.107 × 10 !"	+ V 𝑆=5 𝑗=1 2.550 𝑒𝑎 4 𝐹 𝑏 =* dddd / 𝐹 𝑒𝑏 uuu : 𝐹 × 10 !" cos(𝑗𝜔 𝑚 𝑒𝑎 uuu O 𝑡) + V 𝑆=5 𝑗=1 𝐹 1.695 5.212 𝑒𝑏 4 𝐹 𝑏 𝐹 =e dddd / 𝑒𝑏 uuu O 𝐹 × 10 !" × 10 !# sin(𝑗𝜔 𝑚 𝑒𝑎 uuu , 𝑡) 𝐹 9.095 × 10 !#	𝐹 -4.632 𝑒𝑏 uuu , × 10 !#
	𝑒̅ 2 7 7.908 × 10 !&	Reverse motion 𝑒𝑎 uuu 6 7 𝑒𝑏 uuu 6 7 -3.689 2.249 × 10 !& × 10 !&	𝑒𝑎 uuu ; 7 1.090 × 10 !"	𝑒̅ 𝑚 7 (𝑡) = 𝑒𝑎 uuu : 𝑒 2 7 𝑏 ³ =̅ . 0 7 3.939 × 10 !" 𝑒𝑏 uuu ; 7 7.205 × 10 !"	+ V 𝑆=5 𝑗=1 1.464 𝑒𝑎 4 7 𝑏 =* dddd / 0 𝑒𝑏 uuu : 7 × 10 !" cos(𝑗𝜔 𝑚 𝑒𝑎 uuu O 𝑡) + V 𝑆=5 𝑗=1 7 2.059 -1.043 𝑒𝑏 4 7 0 =e dddd / sin(𝑗𝜔 𝑚 𝑏 𝑒𝑏 uuu O 7 𝑒𝑎 uuu , 𝑡) 7 5.000 × 10 !" × 10 !" × 10 !'	𝑒𝑏 -1.490 7 uuu , × 10 !"
	𝑘 u 2 𝐹 1	Forward motion 𝑎 u 6 𝐹 𝑏 u 6 𝐹 2.173 -3.721 × 10 !( × 10 !&	𝑎 u ; 𝐹 2.277 × 10 !&	𝐾 t 𝑚 𝐹 -4.032 𝑏 u ; 𝐹 × 10 !&	(𝑡) = 𝑎 u : 𝐹 -1.166 𝑘 2 𝐹 𝑘 2 𝐹 ³ 𝐹 f d . × 10 !&	* d / 𝐹 𝑎 4 𝐹 𝑘 2 𝐹 ³ 𝑏 u : + V 𝑆=5 𝑗=1 𝐹 -2.675 cos(𝑗𝜔 𝑚 𝑎 u O 𝑡) + V 𝑆=5 𝑗=1 𝐹 -1.425 -7.900 e d / 𝐹 𝑏 4 𝐹 𝑘 2 𝐹 ³ sin(𝑗𝜔 𝑚 𝑏 u O 𝐹 𝑎 u , 𝑡) 𝐹 -1.063 × 10 !& × 10 !& × 10 !" × 10 !&	𝐹 𝑏 u , 1.992 × 10 !"
	𝑘 u 2 7 4.181 × 10 !(	Reverse motion 𝑎 u 6 7 𝑏 u 6 7 4.104 2.144 × 10 !& × 10 !"	𝑎 u 6 7 -4.970 × 10 !#	𝐾 𝑚 𝑅 (𝑡) = 𝑎 u 6 7 -7.286 𝑏 u 6 7 -3.221 𝑘 2 7 𝑘 2 𝐹 ³ 0 f d . × 10 !" × 10 !"	* d / 0 𝑎 4 7 𝑘 2 𝐹 ³ 𝑏 u 6 + V 𝑆=5 𝑗=1 7 -2.033 cos(𝑗𝜔 𝑚 𝑎 u 6 𝑡) + V 𝑆=5 𝑗=1 7 -3.184 4.325 e d / 0 𝑏 4 7 𝐹 ³ sin(𝑗𝜔 𝑚 𝑘 2 𝑏 u 6 7 𝑎 u 6 𝑡) 7 -7.483 × 10 !" × 10 !" × 10 !# × 10 !#	7 1.767 𝑏 u 6 × 10 !"

Acknowledgments:

The authors thank Dr. Giorgo Bonori and Advanced Numerical Solutions, Inc.

for providing Transmission3D software and their support.

Funding Source Declaration: The financial support provided by the "SUSTAINABLE MOBILITY CENTER (Centro Nazionale per la Mobilità Sostenibile -CNMS)" Spoke 13 "Electric Traction and Batteries (ETB)" project code CN00000023-CUP E93C22001070001

Method of Delays

Method of Delays is a technique of phase space reconstruction. Vectors in a new space, the embedding space, are formed from delay time values of the scalar measurements. The literature on the delay coordinate technique [START_REF] Kennel | Determining embedding dimension for phase-space reconstruction using a geometrical construction[END_REF]- [START_REF] Liu | Chaos characteristics and least squares support vector machines based online pipeline small leakages detection[END_REF] is vast, and for the sake of brevity are not explained in this paper. To reconstruct the trajectories by means of time response, one needs to obtain the embedding dimension, 𝑚, and the delay time or lag, 𝜏. Depending on the type of structure we want to explore a suitable delay time. Small delays yield strongly correlated vector elements, large delays lead to uncorrelated vectors and components, i.e., the data are randomly distributed in the embedding space [START_REF] Hegger | Practical implementation of nonlinear time series methods: The TISEAN package[END_REF]. In practice, the appropriate values for delay time 𝜏, and embedding dimension 𝑚 to use in this reconstruction have a variety of answers. In this study, embedding dimension is defined by "false nearest neighbors" approach [START_REF] Kennel | Determining embedding dimension for phase-space reconstruction using a geometrical construction[END_REF] and the delay time obtained by "mutual information" method [START_REF] Fraser | Independent coordinates for strange attractors from mutual information[END_REF].

• Delay time: The delay time mutual information was suggested by Fraser and Swinney [START_REF] Fraser | Independent coordinates for strange attractors from mutual information[END_REF] as a tool to determine a reasonable delay. Unlike the autocorrelation function, mutual information also takes into account nonlinear correlations. Indeed, the mutual information method measures the general dependence of two variables, regardless of whether it is linear or not; the value of τ giving the firs minimum in the mutual information function can be considered as a good candidate delay time value.

• Embedding dimension: Kennel et al [START_REF] Kennel | Determining embedding dimension for phase-space reconstruction using a geometrical construction[END_REF] introduced a method to determine the minimal sufficient embedding dimension 𝑚. It is called the "false nearest neighbor" method. In an embedding space with a dimension of 𝑚 3 , this implies that the reconstructed attractor perfectly mirrors the attractor in the original phase space in a one-to-one manner.

Especially, the topological properties are preserved. Thus, the neighbors of a given point are mapped onto neighbors in the embedding space. Due to the assumed smoothness of the dynamics, neighborhoods of the points are mapped onto neighborhoods repeatedly.

• For case-II, where the trapping phenomenon occurred, the Poincaré map represents a discontinuous curve with six separate lines, Fig. 17-Case-II. A non-integer correlation dimension of about 1.86 reveals the fractal nature of the attractor. A weak modulation with a complex periodic behavior is observed; therefore, the sign of largest Lyapunov exponent supposed to be positive which is, 𝜆 1UV = 0.054.

• For case-III and case-IV, where the system experiences a strong irregularity in the timehistories, the correlation dimensions are not integer, 2.6 and 2.3 respectively; besides, positive Lyapunov exponents are obtained for both cases.

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.

Supplementary data: Supplementary data and codes associated with this article can be found, in the online version, at: https://doi.org/10.17632/fcpgtmgxsf.1

Appendix A.

geometric parameters of the considered SBG.

Appendix C. a) TISEAN script:

The TISEAN package is available on the website, see Ref. [START_REF] Hegger | TISEAN: Nonlinear Time Series Analysis[END_REF]. Here, the procedure to install the package on Linux-Ubuntu is explained, and the comments to calculate the largest Lyapunov exponent and the correlation dimension is written. The default path for installation is home; therefore, at the end of installation the main files of package will be at home/bin. To do the calculation, the data on timeseries is required. This point should be mentioned that the data must contain only time-history response (y-axis) in one column not the time information as x-axis. The folderPath1 is where the time series response is located (in the comments, "TH.dat" is considered as the DAT file that contain the time series data), and the folderPath2 is where the output files, contained extracted results, must be created. To plot the diagram, "gnuplot" must be installed.

Installation: sudo apt update && sudo apt upgrade sudo apt install gcc dpkg -L gcc gcc --version sudo apt-get install gfortran mkdir -p HOME/bin ./configure --prefix=/home make make install

Maximal Lyapunov exponent:

Step 1: Obtaining key parameters: delay to use


Francesco Bignoli

Andrea Brognara

Philippe Djemia

Damien Faurie

Andrea Li Bassi

Gerhard Dehm

Matteo Ghidelli

Nanoengineering the glassy state: toward novel thin film metallic glasses with outstanding combination of mechanical properties

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Thematic Session: Nanomechanics Keywords: Metallic glasses, Thin films, Micro-mechanics. Disciplinary fields involved: Materials science, Physics, Engineering.

Nanoengineering the glassy state: toward novel thin film metallic glasses with outstanding combination of mechanical properties

Francesco Bignoli Thin film metallic glasses (TFMGs) are emerging materials characterized by a unique combination of mechanical involving large yield strength (~3 GPa) and large ductility (>10%) [1]. Nevertheless, the synthesis of advanced TFMGs with engineered microstructure and the understanding of their mechanical properties are barely tackled. Here, I will present recent results involving two (2) strategies to develop nanoengineered TFMGs with a controlled microstructure down to the atomic scale, resulting in outstanding and tunable mechanical properties. Firstly, I will show the potential of Pulsed Laser Deposition (PLD) as a novel technique to synthetize nanostructured Zr 50 Cu 50 (%at.) TFMGs. I will show how the control of PLD process parameters enables to synthetize a variety of film microstructures, e.g. compact fully amorphous and amorphous nanogranular, showing large free volume (Fig. a) [2]. This results in unique mechanical behavior as shown by in situ TEM/SEM tensile/compression tests, reporting homogeneous deformation for nanogranular TFMGs combined with large yield strength (>3 GPa) and ductility (>9%) [2].

In the second case, I will focus on the fabrication of multilayers with nanoscale period alternating either fully amorphous or amorphous/crystalline sublayers (Fig. b). I will show how the control of the sublayer thickness influences the deformation behavior affecting shear bands formation and tuning the mechanical properties. As an example, alternating CrCoNi (crystalline)/TiZrNbHf (amorphous) nanolayers results in an ultrahigh compressive yield strength (3.6 GPa) and large homogeneous deformation (~15%) [3].