© AFM, EDP Sciences 2018

1 Introduction

Gear transmission system is widely used in several mechanical engineering applications for the purpose of power and motion transmission. Practically, noise and vibration reduction of gear systems remains one of the main goals of researchers. Many excitation sources are responsible for the dynamic behaviour of such transmission. They can be divided into internal or external ones. Concerning internal excitations, time varying mesh stiffness is considered as the main cause since there will a periodic fluctuation of the number of teeth pairs in contact, many authors considered that this variation is the main source of observed noise and vibration [1–3]. For external excitations, variability of loading conditions and rotational speed are one of the main causes leading to both amplitude and frequency modulations [4,5]. Understanding the dynamic behaviour of gear transmission becomes more complicated when nonlinearity occurs [6]. They are induced basically by backlash leading to contact loss between mating gears or in bearings. Spur gear backlash is considered as the amount of tooth space allowing smooth meshing and reduces friction. Dynamic behaviour of gear transmissions in presence of backlash was extensively studied from decades. Several studies [7–9] showed that gear dynamics cannot be predicted with linear model. Munro et al. [10] developed a test rig to measure the dynamic transmission error of a spur gear pair and showed that the tooth separation takes place when the mean load is less than the design load. Kubo et al. [11] observed a jump in the frequency response of the gear pair with backlash. Authors modelled backlash between mating teeth by a discontinuous and non-differentiable function giving rise to a quite complicate nonlinear term in the equation of motion. Numerical techniques should be used in this case. Kahraman and Blankenship [12] used a digital simulation technique and the method of harmonic balance to solve the equation of motion. Wang [13] used lumped parameters models including backlash and solved the differential equation of motion using the piece-wise linear technique which gives only solutions for the equivalent linear systems. Cai and Hayashi [14] presented a nonlinear model in which backlash is considered as a time varying function, resulting in dynamic forces equals to zero when there is no contact between tooth pairs.

Yang et al. [15] developed a nonlinear time-varying dynamic model for right-angle gear pair systems, taking into account both backlash and asymmetric mesh effects. The resolution is done using an enhanced multi-term harmonic balance method with a modified discrete Fourier transform process and the numerical continuation method. They studied the stability of the steady-state harmonic balance solutions and the effects of the variation and asymmetry in mesh stiffness and directional rotation radius on the gear dynamic responses.

Walha et al. [16] investigated the nonlinear dynamics of a two-stage gear system including backlash and time varying mesh stiffness. The nonlinear dynamic response of the system is computed using a linearization technique. He showed clearly that the discontinuity of the motion induced contact loss.

Recently, Moradi and Salarieh [9] studied the nonlinear oscillations of spur gear pairs including the backlash nonlinearity is studied by using classical single degree of freedom model in terms of dynamic transmission. Backlash is considered as a displacement type nonlinear function approximated by a cubic function. They implemented multiple scale method in order to compute forced vibration responses of the gear system including primary, super-harmonic and sub-harmonic resonances.

In the work [17], we have applied an implicit algorithm based on an implicit high order technique and a homotopy transformation. This transformation is performed by introducing a variable change from the initial conditions.

In this paper, there are two original ideas that we will discussed. One is the regularization of expression of the mesh force and the other is the application of the algorithm used in [18–20] for the nonlinear problem of dynamic model of one-stage and two-stage spur gear systems with backlash. This algorithm is based on the association of a high order implicit technique and a homotopy transformation. This homotopy transformation is performed by introducing a variable change from the solution evaluated at a previous time. The proposed algorithm is adapted to strong nonlinearity and is numerically efficient since it requires a reduced computation time. The main idea of the proposed algorithm is to apply the homotopy transformation by introducing an arbitrary invertible matrix [K^*] and an homotopy parameter ϵ varying from zero to one [18–20]. When this parameter is zero, we obtain a simplified problem easy to solve and when the homotopy parameter is equal to one, we recover the initial problem without transformation. By applying the Taylor series technique, unknowns of the nonlinear problem are expanded in power series with respect to the homotopy parameter [21]. This technique allows one to transform the nonlinear problem into a sequence of linear ones which have the same tangent matrix [K^*]. Consequently only one tangent matrix triangulation is needed to compute all the terms of the series for many time steps. A key point of the proposed algorithm is the possibility of choosing the tangent matrix to be decomposed and the possibility of computing many time steps with single tangent matrix decomposition [18–20]. To illustrate the performance of the proposed algorithm, we give numerical comparisons with the classical implicit iterative algorithm of Newton–Raphson type [22].

2 Modelling of gear transmission including backlash

In this section, two dynamic models of gear transmission including tooth backlash, treated in [16], are presented. The objective is to derive the equation of motion of the two models in order to implement the proposed algorithm.

2.1 Dynamic model of a one-stage spur gear system with backlash

A lumped parameters model of a single stage spur gear transmission is presented in Figure 1. The transmission is divided into two parts. The first one is composed of a driving motor, a connecting shaft and the pinion and has a mass m₁. The second one is composed of the wheel and the loading machine connected by the output shaft and has a mass m₂. Motor, pinion, wheel and load are considered as rigid bodies, their respective inertias are I_m, I₁, I₂ and I_L. The two parts are supported by bearings modelled by linear springs $k_{x_1}$, $k_{y_1}$, $k_{x_2}$ and $k_{y_2}$. Connecting shafts are modelled by torsional stiffness $k_{{\theta }_1}$ and $k_{{\theta }_2}$. Each part has four degrees of freedom: two translations according to x and y and two rotations. The degrees of freedom vector can be written as: $$\left\{q\right\}{=}^t\langle x_1,y_1,{\theta }_m,{\theta }_1,x_2,y_2,{\theta }_2,{\theta }_L\rangle$$(1)

Mesh stiffness is modelled by a function $k_{m_i}\left(t\right)$ spring acting along the line of action (see Fig. 2). This model has been used by Fakhfakh et al. in [23]. This mesh stiffness is fluctuating in rectangular form depending on the number of teeth in contact. In this case, the mesh stiffness $k_{m_i}\left(t\right)$ is given as a function of an average component denoted by ${\overline{k}}_{m_i}$ and of a fluctuating component denoted by $k_{v_i}\left(t\right)$. This fluctuating component is given by: $$\begin{array}{ccc}{k}_{v_1}(t)=\{\begin{array}{ccc}{\bar{k}}_{m_1}\frac{2-{\epsilon }_{{\alpha }_1}}{2{\epsilon }_{{\alpha }_1}({\epsilon }_{{\alpha }_1}-1)}\hfill & \mathit{if}\hfill & {\mathit{nT}}_{e_1}\le t\le (n+{\epsilon }_{{\alpha }_1}-1)T_{e_1}\hfill \\ \hfill & \hfill & \hfill \\ -\frac{{\bar{k}}_{m_1}}{2{\epsilon }_{{\alpha }_1}}\hfill & \mathit{if}\hfill & (n+{\epsilon }_{{\alpha }_1}-1)T_{e_1}\le t\le (n+1)T_{e_1}\hfill \end{array}\hfill & \mathit{for}\hfill & \mathit{one}\text{-}\mathit{stage}\hfill \\ k_{v_i}(t)=\{\begin{array}{ccc}{\bar{k}}_{m_i}\frac{2-{\epsilon }_{{\alpha }_i}}{2{\epsilon }_{{\alpha }_i}({\epsilon }_{{\alpha }_i}-1)}\hfill & \mathit{if}\hfill & (n+2-{\epsilon }_{{\alpha }_i})T_{e_i}\le t\le (n+1)T_{e_i}\hfill \\ \hfill & \hfill & \hfill \\ 2212;\frac{{\bar{k}}_{m_i}}{2{\epsilon }_{{\alpha }_i}}\hfill & \mathit{if}\hfill & {\mathit{nT}}_{e_i}\le t\le (n+2-{\epsilon }_{{\alpha }_i})T_{e_i}\hfill \end{array}\hfill & \mathit{for}\hfill & \mathit{two}\text{-}\mathit{stage}\hfill \end{array}$$(2) where n is an integer which designates the number of the period, i designates the stage number, the terms ${ϵ}_{{\alpha }_i}$ is the contact ratio at stage i, $T_{e_i}$ is the period and the average component ${\overline{\mathit{k}}}_{m_i}=k_{\mathit{max}}({ϵ}_{{\alpha }_i}-1)+k_{\mathit{min}}(2-{ϵ}_{{\alpha }_i})$, $T_{e_2}=\frac{Z_2}{Z_3}{T}_{e_1}$ in the case of two-stage; with Z₂ and Z₃ are the teeth numbers. The reliability of the physical results is dependent on the validity of the model and that definite validation is still to come. An input torque C_m is applied by the motor and a resisting torque C_L is opposed by the load. The expression of the transmission error describing the displacement along the line of action of mating teeth can be expressed as [24]: $$\delta \left(t\right)=\left(x_1-x_2\right)\sin \left(\alpha \right)+\left(y_1-y_2\right)\cos (\alpha )+{\theta }_1{r}_{b_1}+{\theta }_2{r}_{b_2}{=}^T\{\chi \}\{q\}$$(3)where α is the pressure angle, $r_{b_1}$ and $r_{b_2}$ are respectively the base radius of pinion and wheel, {χ} is defined by: $${}^T\left\{\chi \right\}=\langle \sin (\alpha ),\cos (\alpha ),0,r_{b_1},-\sin (\alpha ),-\cos (\alpha ),r_{b_2},0\rangle$$(4)and, here, the degrees of freedom θ₁ and θ₂ are algebraic quantities. A backlash between teeth is considered. To take into account of this backlash, one can write the expression of the mesh force F_s as following: $$F_s=\{\begin{array}{ccc}{k}_{m_i}(t)(\delta (t)+b_s)\hfill & ,\mathit{if}\hfill & \delta (t)<-b_s\hfill \\ \hfill & \hfill & \hfill & \hfill \\ 0,\hfill & \mathit{if}\hfill & -b_s\le \delta (t)\le b_s\hfill \\ \hfill & \hfill & \hfill & \hfill \\ k_{m_i}(t)(\delta (t)-b_s)\hfill & ,\mathit{if}\hfill & \delta (t)>b_s\hfill \end{array}$$(5)where 2b_s is the total backlash as shown in Figure 3.

Writing Lagrange formulation allows obtaining the equation of motion expressed by: $$\left[M\right]\left\{\ddot{q}\right\}+\left[C\right]\left\{\stackrel{}{q}\right\}+\left[K_b\right]\left\{q\right\}+F_s\left(\delta \left(t\right)\right)\left\{\chi \right\}=\left\{F_{\mathit{ext}}\right\}$$(6) where [M] is the inertia matrix given by: $$[M]=\left[\begin{array}{cccccccc}\hfill m_1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill m_1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill I_m\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill I_1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill m_2\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill m_2\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill I_2\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill I_L\hfill \end{array}\right]$$(7)

[K_b] is the bearings and shafts matrix expressed by: $$[K_b]=\left[\begin{array}{cccccccc}\hfill k_{x_1}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill k_{y_1}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill k_{{\theta }_1}\hfill & \hfill -k_{{\theta }_1}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -k_{{\theta }_1}\hfill & \hfill k_{{\theta }_1}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill k_{x_2}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill k_{y_2}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill k_{{\theta }_2}\hfill & \hfill -k_{{\theta }_2}\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -k_{{\theta }_2}\hfill & \hfill k_{{\theta }_2}\hfill \end{array}\right]$$(8)

The vector of external applied forces {F_ext} is expressed by: $$\left\{F_{\mathit{ext}}\right\}{=}^t\langle 0,0,C_m,0,0,0,0,C_Z\rangle$$(9)

The damping matrix [C] in this model is proportional both to the mass and stiffness matrices of the system given by: $$\left[C\right]=c_1\left[M\right]+c_2\left[K_{\mathit{moy}}\right]$$(10)

where c₁ and c₂ are the damping coefficients which are expressed respectively in s⁻¹ and s which can be determined from two damping ratios specified as in [25] and [K_moy] is the average stiffness matrix of the system [16]. The initial conditions used for this problem are: {q} = {0} and $\left\{\stackrel{}{q}\right\}=\left\{0\right\}$.

Fig. 1 Single stage spur gear system model.

Fig. 2 Spring acting along the line of action.

Fig. 3 Contact loss modelling.

2.2 Dynamic model of a two-stage spur gear system with backlash

A model of a two-stage spur gear transmission is presented in Figure 4. It can be divided into three parts. Part 1 is composed of driving unit, input shaft and pinion 1. Part 2 includes wheel 1 and pinion 2 connected by the intermediate shaft. Part 3 is composed of wheel 2, output shaft and load.

The system has 12 degrees of freedom (DOF) which can be detailed as follows:

• Translations of part 1 (having the mass m₁) along x (horizontal) and y (vertical) directions. The corresponding DOF are x₁ and y₁.

• Translations of part 2 (having the mass m₂) x and y directions. The corresponding DOF are x₂ and y₂.

• Translations of part 3 (having the mass m₃) along x and y directions. The corresponding DOF are x₃ and y₃.

• Rotations of motor θ_m, pinion 1 θ₁, wheel 1 θ₂, pinion 2 θ₃, wheel 2 θ₄ and load θ_L.

The following inertia are considered: I_m for the motor, I₁ for the pinion 1, I₂ for wheel 1, I₃ for pinion 2, I₄ for wheel 2 and I_L for load. Input, output and intermediate shafts elasticity is modelled by torsional stiffness respectively $k_{{\theta }_1}$, $k_{{\theta }_2}$ and $k_{{\theta }_3}$. Bearings supporting input shaft are modelled by linear springs $k_{x_1}$ and $k_{y_1}$ acting along x and y axis. Bearings supporting intermediate shaft are modelled by linear springs $k_{x_2}$ and $k_{y_2}$. Bearings supporting intermediate shaft are modelled by linear spring $k_{x_3}$ and $k_{y_3}$.

Two mesh stiffness functions $k_{m_1}\left(t\right)$ and $k_{m_2}\left(t\right)$ fluctuating in a rectangular shape are introduced between pinion 1 and wheel 1 and between pinion 2 and wheel 2 (see Fig. 5). If we consider backlash for the two-stage, two mesh forces should be considered $F_{s_1}$ and $F_{s_2}$. One should be aware here for phasing that can occur between the two mesh functions depending on the position of stage 1 with respect to stage 2.

The vector defining the degrees of freedom {q} is expressed by: $$\left\{q\right\}{=}^t\langle x_1,y_1,{\theta }_m,{\theta }_1,x_2,y_2,{\theta }_2,{\theta }_3,x_3,y_3,{\theta }_4,{\theta }_L\rangle$$(11)

Two transmission errors δ₁(t) and δ₂(t) for the two-stage are expressed by [16]: $$\{\begin{array}{ccc}{\delta }_1(t)\hfill & =\hfill & (x_1-x_2)\hspace{0.17em}\text{sin}\hspace{0.17em}({\alpha }_1)+(y_1-y_2)\hspace{0.17em}\text{cos}\hspace{0.17em}({\alpha }_1)\hfill \\ & & +r_{b_1}{\theta }_1-r_{b_2}{\theta }_2={}^T\{_{{\chi }_1}^{\}}\{q\}\hfill \\ {\delta }_2(t)\hfill & =\hfill & (x_2-x_3)\hspace{0.17em}\text{sin}\hspace{0.17em}({\alpha }_2)-(y_2-y_3)\hspace{0.17em}\text{cos}\hspace{0.17em}({\alpha }_2)\hfill \\ & & +r_{b_3}{\theta }_3-r_{b_4}{\theta }_4={}^T\{_{{\chi }_2}^{\}}\{q\}\hfill \end{array}$$(12) where α₁ and α₂ are the pressure angles, {χ₁} and {χ₂} are given by: $$\{\begin{array}{c}{}^T\{{\chi }_1\}=<\hspace{0.17em}\text{sin}\hspace{0.17em}({\alpha }_1),\hspace{0.17em}\text{cos}\hspace{0.17em}({\alpha }_1),0,r_{b_1},-\hspace{0.17em}\text{sin}\hspace{0.17em}({\alpha }_1),-\hspace{0.17em}\text{cos}\hspace{0.17em}({\alpha }_1),-r_{b_2},0,0,0,0,0>\hfill \\ \hfill \\ {}^T\{{\chi }_2\}=<0,0,0,0,\hspace{0.17em}\text{sin}\hspace{0.17em}({\alpha }_2),-\hspace{0.17em}\text{cos}\hspace{0.17em}({\alpha }_2),0,r_{b_3},-\hspace{0.17em}\text{sin}\hspace{0.17em}({\alpha }_2),\hspace{0.17em}\text{cos}\hspace{0.17em}({\alpha }_2),-r_{b_4},0>\hfill \end{array}$$(13)and, here, the degrees of freedom θ₁, θ₂, θ₃ and θ₄ are positive quantities. Taking into account of the backlashes in the two mesh zones, the expression of the mesh forces $F_{s_1}$ and $F_{s_2}$ respectively measured on the meshes in stage 1 and stage 2 are as following: $$F_{s_1}=\{\begin{array}{ccc}{k}_{m_1}(t)({\delta }_1(t)+b_1)\hfill & ,\mathit{if}\hfill & {\delta }_1(t)<-b_1\hfill \\ \hfill & \hfill & \hfill & \hfill \\ 0,\hfill & \mathit{if}\hfill & -b_1\le {\delta }_1(t)\le b_1\hfill \\ \hfill & \hfill & \hfill & \hfill \\ k_{m_1}(t)({\delta }_1(t)-b_1)\hfill & ,\mathit{if}\hfill & {\delta }_1(t)>b_1\hfill \end{array}$$(14) $$F_{s_2}=\{\begin{array}{ccc}{k}_{m_2}(t)({\delta }_2(t)+b_2)\hfill & ,\mathit{if}\hfill & {\delta }_2(t)<-b_2\hfill \\ \hfill & \hfill & \hfill & \hfill \\ 0,\hfill & \mathit{if}\hfill & -b_2\le {\delta }_2(t)\le b_2\hfill \\ \hfill & \hfill & \hfill & \hfill \\ k_{m_2}(t)({\delta }_2(t)-b_2)\hfill & ,\mathit{if}\hfill & {\delta }_2(t)>b_2\hfill \end{array}$$(15)

By the same Lagrange formulation it is possible to obtain the equation of motion of the two-stage gear system by: $$\left[M\right]\left\{\ddot{q}\right\}+\left[C\right]\left\{\stackrel{}{q}\right\}+{\left[K\right]}_b\left\{q\right\}+F_{s_1}\left({\delta }_1\right)\left\{{\chi }_1\right\}+F_{s_2}\left({\delta }_2\left(t\right)\right)\left\{{\chi }_2\right\}=\left\{F_{\mathit{ext}}\right\}$$(16) where [M] is the inertia matrix given by: $$[M]=\left[\begin{array}{cccccccccccc}\hfill m_1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill m_1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill I_m\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill I_1\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill m_2\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill m_2\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill I_2\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill I_3\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill m_3\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill m_3\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill I_4\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill I_L\hfill \end{array}\right]$$(17)

[K_b] is the bearing and shafts matrix expressed by: $$[K_b]=\left[\begin{array}{cccccccccccc}\hfill k_{x_1}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill k_{y_1}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill k_{{\theta }_1}\hfill & \hfill -k_{{\theta }_1}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill -k_{{\theta }_1}\hfill & \hfill k_{{\theta }_1}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill k_{x_2}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill k_{y_2}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill k_{{\theta }_2}\hfill & \hfill -k_{{\theta }_2}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -k_{{\theta }_2}\hfill & \hfill k_{{\theta }_2}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill k_{x_3}\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill k_{y_3}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill k_{{\theta }_3}\hfill & \hfill -k_{{\theta }_3}\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill & \hfill -k_{{\theta }_3}\hfill & \hfill k_{{\theta }_3}\hfill \end{array}\right]$$(18)

The vector of external applied forces{F_ext} is expressed by: $$\left\{F_{\mathit{ext}}\right\}{=}^t\left\{0,0,C_m,0,0,0,0,0,0,0,0,C_L\right\}$$(19)

The initial conditions used for this problem are: {q} = {0} and $\left\{\stackrel{}{q}\right\}=\left\{0\right\}$. So, the initial conditions of transmission errors deduced from equations (3) and (12) are given by: $$\begin{array}{cc}\left\{\begin{array}{c}\delta (0)=0\hfill \end{array}\right\}\hfill & \mathit{for\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\; one}\text{-}\mathit{stage}\hfill \\ \left\{\begin{array}{c}{\delta }_1(0)=0\hfill \\ {\delta }_2(0)=0\hfill \end{array}\right\}\hfill & \mathit{for\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\; two}\text{-}\mathit{stage}\hfill \end{array}$$(20)

Fig. 4 Two-stage spur gear system model.

Fig. 5 Spring acting along the line of action.

3 The proposed algorithm

The proposed algorithm for solving the problem (6) or (16) is based on four steps: a regularization technique, a time discretization, the introduction of a perturbation from the solution at a previous time, a homotopy transformation and a Taylor series technique. The presented algorithm includes a generic procedure and it is a generic solver that can be applied to a wide class of nonlinear unsteady equations as equations (6) and (16). The Taylor series technique allows us to calculate all terms of the series with a single inversion of the tangent matrix.

3.1 Regularization technique

In the case of problems (6) and (16), non-smooth functions appear F_s(δ(t)) defined by (5), $F_{s_1}\left({\delta }_1\left(t\right)\right)$ defined by (14) and $F_{s_2}\left({\delta }_2\left(t\right)\right)$ defined by (15). The idea is to replace these non smooth functions by a smooth function denoted by: $${\tilde{F}}_s\left(\delta \left(t\right)\right)=k_{m_i}\delta \left(t\right)+\frac{1}{2}\left(\sqrt{k_{m_i}^2{\left(\delta \left(t\right)-b_s\right)}^2+4{\eta }^2{b}_s^2}-\sqrt{k_{m_i}^2{\left(\delta \left(t\right)+b_s\right)}^2+4{\eta }^2{b}_s^2}\right)$$(21)

The function defined here involves a regularization parameter η. This parameter is chosen significantly small in order that the modified function is close to the non smooth function defined by (5) (see Fig. 6). To make easy the management of the regularization, we have defined dimensionless regularization parameter.

To obtain a quadratic problem, we introduce the new variables h(δ(t)) and g(δ(t)) defined by: $$\{\begin{array}{c}h(\delta (t))=\sqrt{k_{m_i}^2{(\delta (t)-b_s)}^2+4{\eta }^2{b}_s^2}\hfill \\ \hfill \\ g(\delta (t))=\sqrt{k_{m_i}^2{(\delta (t)+b_s)}^2+4{\eta }^2{b}_s^2}\hfill \end{array}$$(22)

Taking into account the equation (22), the modified function (21) is set in the following form: $${\tilde{F}}_s\left(\delta \left(t\right)\right)=k_{m_i}\delta \left(t\right)+\frac{1}{2}\left(h\left(\delta \left(t\right)\right)-g\left(\delta \left(t\right)\right)\right)$$(23)

Fig. 6 The approximated of the $({\tilde{F}}_s\delta (t))$ based on a regularization technique for different values of η.

3.2 Time discretization

For solving the nonlinear unsteady problem (16), we use here the implicit Newmark scheme widely used in the resolution of dynamic problems [26]. This very popular implicit scheme requires less, but more expensive time steps to follow the equilibrium path of a dynamic system [26]. This scheme gives the expression of speed $\left\{{\stackrel{}{q}}^{n+1}\right\}$ and acceleration $\left\{{\ddot{q}}^{n+1}\right\}$ at time (n + 1)Δt as follows: $$\{\begin{array}{l}\left\{{\ddot{q}}^{n+1}\right\}=a_0(\{q^{n+1}\}-\{q^n\})-\left\{{\psi }^n\right\}\hfill \\ \left\{{\ddot{q}}^{n+1}\right\}=\Delta t\beta a_0(\{q^{n+1}\}-\{q^n\})+\left\{w^n\right\}\hfill \end{array},$$(24) where $\left\{{\omega }^n\right\}=a_1\left\{{\stackrel{}{q}}^n\right\}+a_2\left\{{\ddot{q}}^n\right\}$, $\left\{{\psi }^n\right\}=\left(1-a_1\mathrm{\Delta }t\beta \right)\left\{{\stackrel{}{q}}^n\right\}+\left(\mathrm{\Delta }t\left(1-\beta \right)-a_2\mathrm{\Delta }t\beta \right)\left\{{\ddot{q}}^n\right\}$, $a_0=\frac{1}{\gamma \mathrm{\Delta }{t}^2}$, $a_1=\frac{1}{\gamma \mathrm{\Delta }t}$, $a_2=\frac{1}{2\gamma }-1$, n denotes the time level, Δt is the time step, the constants β and γ have the range 0 ≤ β ≤ 1, $0\le \gamma \le \frac{1}{2}$ with typical values being $\beta =\frac{1}{2}$, $\gamma =\frac{1}{4}$. Using equation (24), the problem (6) at time (n + 1)Δt is written as: $$\left(a_0\left[M\right]+\beta a_0\mathrm{\Delta }t\left[C\right]\right)(\{q^{n+1}\}-\{q^n\})+\left[K_b\right]\{q^{n+1}\}+{\tilde{F}}_s({\delta }^{n+1})\left\{\chi \right\}=\left\{F_{\mathit{ext}}\right\}+\left[M\right]\left\{{\psi }^n\right\}-\left[C\right]\left\{{\omega }^n\right\}$$(25)

where {qⁿ⁺¹}, δⁿ⁺¹, hⁿ⁺¹ and gⁿ⁺¹ are the unknowns at time t = (n + 1)Δt.

3.3 Introduction of a perturbation from the previous time

Instead of searching the unknowns {qⁿ⁺¹}, δⁿ⁺¹, hⁿ⁺¹ and gⁿ⁺¹ one seeks the perturbations {Δq}, ΔD, ΔH and ΔG from the solutions {qⁿ}, δⁿ, hⁿ and gⁿ at time t = nΔt in the following manner. $$\{\begin{array}{c}\{q^{n+1}\}=\{q^n\}+\{\mathrm{\Delta }Q\}\hfill \\ \hfill \\ {\delta }^{n+1}={\delta }^n+\mathrm{\Delta }D\hfill \\ \hfill \\ h^{n+1}=h^n+\mathrm{\Delta }H\hfill \\ \hfill \\ g^{n+1}=g^n+\mathrm{\Delta }G\hfill \end{array}$$(26)

The expressions (26) are introduced into equations (3), (22), and (25), we obtain a nonlinear problem: $$\{\begin{array}{c}[K_T^n]\{\mathrm{\Delta }Q\}+\{\mathit{FQ}\}=\{S^n\}\hfill \\ \mathrm{\Delta }D={}^T\{\chi \}\{\mathrm{\Delta }Q\}\hfill \\ 2h^n\mathrm{\Delta }H=k_{m_i}^{{n+1}^2}{({\delta }^n-b_s)}^2+4{\eta }^2{b}_s^2-h^{n^2}+2k_{m_i}^{{n+1}^2}({\delta }^n-b_s)\mathrm{\Delta }D+k_{m_i}^{{n+1}^2}\mathrm{\Delta }{D}^2-\mathrm{\Delta }{H}^2\hfill \\ 2g^n\mathrm{\Delta }G=k_{m_i}^{{n+1}^2}{({\delta }^n+b_s)}^2+4{\eta }^2{b}_s^2-g^{n^2}+2k_{m_i}^{{n+1}^2}({\delta }^n+b_s)\mathrm{\Delta }D+k_{m_i}^{{n+1}^2}\mathrm{\Delta }{D}^2-\mathrm{\Delta }{G}^2\hfill \end{array}$$(27) where $\left[K_T^n\right]$ is the tangent matrix given by $$\left[K_T^n\right]=a_0\left[M\right]+\mathrm{\Delta }t\beta a_0\left[C\right]+\left[K_b\right]+\left[K_{\mathit{en}}^n\right]$$(28)with $$\left[K_{\mathit{en}}^n\right]=\left(k_{m_i}^{n+1}+\frac{{k_{m_i}^{n+1}}^2}{2}\left(\frac{({\delta }^n-b_s)}{h^n}-\frac{({\delta }^n+b_s)}{g^n}\right)\right){\left\{\chi \right\}}^T\left\{\chi \right\},$$(29)

{FQ} is a quadratic form defined by: $$\left\{\mathit{FQ}\right\}=\frac{1}{4}(k_{m_i}^{{n+1}^2}(\frac{1}{h^n}-\frac{1}{g^n})\mathrm{\Delta }{D}^2-\frac{\mathrm{\Delta }{H}^2}{h^n}+\frac{\mathrm{\Delta }{G}^2}{g^n})\left\{\chi \right\}$$(30) and {Sⁿ} is the second hand side expressed by: $$\{S^n\}=\{F_{\mathit{ext}}\}+[M]\{{\psi }^n\}-[C]\{{\omega }^n\}-[K_b]\{q^n\}-(k_{m_i}^{n+1}{\delta }^n+\frac{1}{2}(h^n-g^n))\{\chi \}-\frac{1}{4}(k_{m_i}^{{n+1}^2}(\frac{{({\delta }^n-b_s)}^2}{h^n}-\frac{{({\delta }^n+b_s)}^2}{g^n})+4{\eta }^2{b}_s^2(\frac{1}{h^n}-\frac{1}{g^n})+g^n-h^n)\{\chi \}$$(31)

3.4 Homotopy transformation

It is well known that the homotopy transformations have been developed to overcome the local convergence of the Newton-like methods and may provide a reliable way to solve the nonlinear equations [21]. In this work, a homotopy transformation is used while introducing an invertible matrix [K^*] in the problem (27)a and the parameter ϵ of this transformation is considered as an expanding parameter [18–20]: $$[K^{*}]\left\{\mathrm{\Delta }W\right\}+ϵ((\left[K_T^n\right]-[K^{*}])\left\{\mathrm{\Delta }W\right\}+\left\{\mathit{FQ}\right\})=ϵ\left\{S^n\right\}$$(32)

By this way, the solutions ΔW(ϵ) of the nonlinear equation (32) passe continuously from 0 for ϵ = 0 to the solution of the nonlinear problem (27) for ϵ = 1.

3.5 Taylor series technique

The basic idea of this technique consists in searching the solution branches of the nonlinear problem (32) in the form of a truncated Taylor expansion: $$\{\begin{array}{c}\{\mathrm{\Delta }W(\epsilon )\}=\sum _{i=0}^{i=p}{\epsilon }^i\{\mathrm{\Delta }{W}_i\}\hfill \\ \mathrm{\Delta }D(\epsilon )=\sum _{i=0}^{i=p}{\epsilon }^i\mathrm{\Delta }{D}_i\hfill \\ \mathrm{\Delta }H(\epsilon )=\sum _{i=0}^{i=p}{\epsilon }^i\mathrm{\Delta }{H}_i\hfill \\ \mathrm{\Delta }G(\epsilon )=\sum _{i=0}^{i=p}{\epsilon }^i\mathrm{\Delta }{G}_i\hfill \end{array}$$(33)

By introducing the expansions (33) into equation (32) and by equating like powers of ϵ, we obtain a set of linear problems satisfied by the terms of the series (33) given by: $$p=1:\{\begin{array}{c}[K^{*}]\{\mathrm{\Delta }{W}_1\}=([K_T^{n+1}]-[K^{*}])\{\mathrm{\Delta }{W}_0\}-\{\mathit{FQ}(\mathrm{\Delta }{W}_0,\mathrm{\Delta }{W}_0)\}+\{S^{n+1}\}\hfill \\ \mathrm{\Delta }{D}_1={}^T\{\chi \}\{\mathrm{\Delta }{W}_1\}\hfill \\ \mathrm{\Delta }{H}_1=\frac{k_{m_i}^{{n+1}^2}({\delta }^n+\mathrm{\Delta }{D}_0-b_s)}{h^n+\mathrm{\Delta }{H}_0}\mathrm{\Delta }{D}_1\hfill \\ \mathrm{\Delta }{G}_1=\frac{k_{m_i}^{{n+1}^2}({\delta }^n+\mathrm{\Delta }{D}_0+b_s)}{g^n+\mathrm{\Delta }{G}_0}\mathrm{\Delta }{D}_1\hfill \end{array}$$(34) $$p\ge 2:\{\begin{array}{c}[K^{*}]\{\mathrm{\Delta }{W}_p\}=-([K_T^{n+1}]-[K^{*}])\{\mathrm{\Delta }{W}_{p-1}\}-\sum _{r=0}^{p-1}\{\mathit{FQ}(\mathrm{\Delta }{W}_r,\mathrm{\Delta }{W}_{p-1-r})\}\hfill \\ \mathrm{\Delta }{D}_p={}^T\{\chi \}\{\mathrm{\Delta }{W}_p\}\hfill \\ \mathrm{\Delta }{H}_p=\frac{2k_{m_i}^{{n+1}^2}({\delta }^n+\mathrm{\Delta }{D}_0-b_s)\mathrm{\Delta }{D}_p+k_{m_i}^{{n+1}^2}{\sum }_{r=1}^{p-1}\mathrm{\Delta }{D}_r\mathrm{\Delta }{D}_{p-r}-{\sum }_{r=1}^{p-1}\mathrm{\Delta }{H}_r\mathrm{\Delta }{H}_{p-r}}{2(h^n+\mathrm{\Delta }{H}_0)}\hfill \\ \mathrm{\Delta }{G}_p=\frac{2k_{m_i}^{{n+1}^2}({\delta }^n+\mathrm{\Delta }{D}_0+b_s)\mathrm{\Delta }{D}_p+k_{m_i}^{{n+1}^2}{\sum }_{r=1}^{p-1}\mathrm{\Delta }{D}_r\mathrm{\Delta }{D}_{p-r}-{\sum }_{r=1}^{p-1}\mathrm{\Delta }{G}_r\mathrm{\Delta }{G}_{p-r}}{2(g^n+\mathrm{\Delta }{G}_0)}\hfill \end{array}$$(35)

The solution is obtained step by step. Each step is represented by the power series (33). The length of each step is computed by requiring that the relative difference between the solutions with two consecutive truncation orders must remain small enough by comparison to a given accuracy parameter κ. The step length is then computed by defining a maximal value of the parameter ϵ_max [18–20] as follows: $${\epsilon }_{\mathit{max}}={\left(\kappa \frac{\parallel \{\mathrm{\Delta }{W}_1\}\parallel }{\parallel \{\mathrm{\Delta }{W}_p\}\parallel }\right)}^{\frac{1}{p-1}}$$(36)

The solution of problem (25) is obtained while the parameter ϵ_max ≥ 1.

4 Numerical applications

In this work, we present two typical examples of the gear system with one-stage and two-stage. The two nonlinear problems are solved by the detailed algorithm above. In both examples, the matrix [K^*] has the same form as $\left[K_T^n\right]$ defined at the initial conditions. If $K^{*}=K_T^n$, our algorithm requires the inversion of this matrix each time step which demands a considerable computation time. This way of doing is similar to the iterative and incremental method type Newton–Raphson.

4.1 Example 1: one-stage gear system

Let us consider a spur gear system driven by a squirrel cage motor. Its characteristics are given in Table 1.

In this application, we have chosen a pre-conditioner [K^*] which is taken equal to the tangent matrix $\left[K_T^n\right]$ evaluated at the initial time, a time step Δt = 3.7 × 10⁻⁶ s, a tolerance parameter κ = 10⁻⁴ and various values for the truncation order for our proposed algorithm, a tolerance parameter κ = 10⁻⁴ for Newton–Raphson algorithm and a regularization parameter η = 10⁻⁶. The time interval for this study is [0, 0.09 s]. In a first step, we studied the influence of the truncation order on the solution quality which is defined from the residual vector: $$\left\{R\right\}=\left[K_T^n\right]\left\{\mathrm{\Delta }Q\right\}+\left\{\mathrm{FQ}\left(\mathrm{\Delta }Q,\mathrm{\Delta }Q\right)\right\}-\left\{S^n\right\}$$(37)

The solution quality is measured by the decimal logarithm of the residual vector R. In Figure 7, we represent the evolution of the solution quality versus the time. We remark that this quality increases with the truncation order of the series for p = 3, 4, 5, 6. The curve of Figure 7 is obtained by the used algorithm in a single step of continuation for orders greater or equal than 3. So the stage of the continuation of procedure is not used.

In this numerical test, we study the influence of several values of the regularization parameter η on the solution of the problem. In Figure 8, we represent the transmission error versus to time for η = 10⁻², η = 10⁻³, η = 10⁻⁴, η = 10⁻⁵ and η = 10⁻⁶. According to this test, we noticed that the solution of the problem is stabilized starting from the value η = 10⁻³ (see Fig. 8).

In the following, we fix a truncation order p = 3, a tolerance κ = 10⁻⁴ for our algorithm and a tolerance κ = 10⁻⁴ for the Newton–Raphson method for a quality solution less than 10⁻⁴ and a regularization parameter η = 10⁻³. In Figure 9, we represent the transmission error obtained by the Newton–Raphson method and the proposed algorithm versus time. From this figure we see that the solution reveals the contact losses at startup. We can say that the clearance between teeth has any influence at the beginning of the movement.

In Table 2, we present the number of triangulations of tangent matrix for the both methods, i.e. the proposed algorithm and the Newton–Raphson method.

The complete solution along the considered time interval obtained by our algorithm requires only one inversion of the matrix [K^*] while that obtained by Newton–Raphson has requested 73 709 inversions of the tangent matrix.

Figure 10 represents the temporary evolution of the linear displacement of the first and second bearing. As shown in Figure 10, the detailed comparison between the used algorithm and the Newton–Raphson results are made. The used algorithm results are in good agreement with the Newton–Raphson results.

Fig. 7 Evolution of the solution quality log10||{R}|| versus the time t.

Fig. 8 Influence of regularization parameter on the evolution of the transmission error δ versus the time t for η = 10−2, η = 10−3, η = 10−4, η = 10−5 and η = 10−6.

Fig. 9 Time evolution of transmission error, Proposed algorithm truncated at order p = 3 and κ = 10−4 (ANM order 3), Newton–Raphson method with tolerance κ = 10−4 (Newton–Raphson).

Fig. 10 Temporary evolution of the linear displacement of the first and second bearing, Proposed algorithm truncated at order p = 3 and κ = 10−4 (ANM order 3), Newton–Raphson method with κ = 10−4 (Newton–Raphson).

4.2 Example 2: two-stage gear system

In this example, the nonlinear behaviour of a two-stage gear reducer is surveyed while supposing that the backlash is localized to the two trains [16]. In this example, we will take into account the presence of the clearance between teeth as well as changes of the number of teeth into contact in the modeling problem. In order to see the impact of the introduction of functional clearance and its location, we will compare the behavior of the mechanism in both situations, the presence of functional clearance in the first meshing or both meshing. The characteristics of this example are as in [16]. In Table 3, it brings together all the parameters of the two-stage gear system.

4.2.1 Case of one backlash at the first stage: b₁ = b/2 and b₂ = 0

We chose a preconditioner [K^*] Corresponding to the tangent matrix evaluated at the initial time. The time interval for this study is [0, 0.02 s] and the time steps is 4 × 10⁻⁶ s. For this, we fix the tolerance parameter of the two algorithms, i.e. Newton–Raphson algorithm and the proposed algorithm to 10⁻⁴ and η₁ = η₂ = 10⁻⁶. In order to choose the optimal truncation order p that will be used in the following, we evaluate the evolution of the solution quality according to the truncation order. Figure 11 shows that as the truncation order p increases, the solution quality becomes best. Beyond the order p = 5, the solution quality becomes insensitive.

In the following, we choose a truncation order p = 3 and a tolerance parameter κ = 10⁻⁴ for the proposed algorithm and a tolerance parameter κ = 10⁻⁴ for the Newton–Raphson method.

In Figure 12a and b, we represent the time evolution of transmission errors δ₁ and δ₂ at the first and second meshing respectively in the presence of backlash on the first stage. The second teeth deflection fluctuates around zeros but the first nonlinear one fluctuates with the top of the half of backlash.

In Figure 13, we represent the linear and nonlinear angular velocities fluctuation of the input wheel, the middle gear and the output wheel with respect to time. After a small time due to the first backlash, the angular velocities fluctuate around the linear responses.

In Figures 14 and 15, we present the linear (b₁ = 0, b₂ = 0) and nonlinear (b₁ = b/2, b₂ = 0) speeds ẋ₁ and ẋ₃ of the first and third block along the axis x and the corresponding accelerations ẍ₁ and ẍ₃ versus time.

In Figure 16, we represent the time relative teeth velocity fluctuation following teeth displacement in the case of one backlash located on the first stage gear.

In Table 4, we present the number of triangulations of tangent matrix for the both methods, i.e. the proposed algorithm and the Newton–Raphson method.

The complete solution along the considered time interval obtained by our algorithm requires only one inversion of the matrix [K^*] while that obtained by Newton–Raphson has requested 7586 inversions of the tangent matrix.

Fig. 11 Decimal logarithm of the norm the residual vector versus time, Proposed algorithm with different truncation orders p = 3, 4, 5, 6, 7, 8 (ANM order p) and κ = 10−4, Newton–Raphson algorithm with κ = 10−4 (NR).

Fig. 12 Temporal evolution of transmission errors, proposed algorithm with truncation order p = 3 (ANM) and κ = 10−4, Newton–Raphson algorithm with κ = 10−7 (NR) for (b1 = 0, b2 = 0) and (b1 = 10−4, b2 = 0).

Fig. 13 Temporal evolution of linear and nonlinear angular velocities fluctuation.

Fig. 14 Temporal evolution of bearings ẋ1 and ẋ3.

Fig. 15 Temporal evolution of bearings ẍ1 and ẍ3.

Fig. 16 Relative teeth velocity fluctuation following teeth displacement in the case of one backlash located on the first stage.

4.2.2 Case of two backlashes at the first and second stages: b₁ = b/2 and b₂ = b/2

In this second part, we take into account the presence of two backlashes at both stages (b₁ = b/2, b₂ = b/2). As previously, we choose the solution quality represented in Figure 17 versus the time for the truncation order p = 3 and the tolerance κ = 10⁻⁴ for both algorithms.

In Figure 18, we present the evolution of transmission error versus the time at the first and second stages in the presence of two backlashes. From these solutions curves, we note that during the transitional regime, there is occurrence of several states of contact loss in the first meshing. While the transitional regime of the second meshing is characterized on the one hand by a delay of 0.004 s and on the other hand by the appearance of a larger contact loss which manifests at time 0.007 s and it lasts until time 0.01 s.

Figure 19 shows the evolution of the angular velocities of each meshing versus the time. For both angular velocities ${\stackrel{}{\theta }}_1$ and ${\stackrel{}{\theta }}_3$, we notice the same comments about the delay and aspect of the transitional regime that in the case of a single backlash. While the transitional regime of angular velocity ${\stackrel{}{\theta }}_4$ tripped with a delay of 0.004 s. Similarly, the solution curve of this angular velocity continues to oscillate around the curve obtained without contact until the stabilization which starts from the time 0.02 s.

Figure 20 shows clearly the influence of the presence of two backlashes on the solution curves plotted in the phase space. This influence is characterized by the change of the position of the attractors centers and by the instability during the transient regime.

In Table 5 and in the same manner that in the previous numerical test, we present the number of triangulations of tangent matrix for the both methods.

The complete solution along the considered time interval obtained by our algorithm requires only one inversion of the matrix [K^*] while that obtained by Newton–Raphson has requested 8073 inversions of the tangent matrix.

Fig. 17 Decimal logarithm of the norm the residual vector versus the time, Proposed algorithm with truncation order p = 3 (ANM order 3) and κ = 10−4, Newton–Raphson algorithm with κ = 10−4 (NR).

Fig. 18 Temporal evolution of transmission errors, proposed algorithm with truncation order p = 3 (ANM) and κ = 10−4, Newton–Raphson algorithm with κ = 10−7 (NR) for (b1 = 0, b2 = 0) and (b1 = 10−4, b2 = 10−4).

Fig. 19 Temporal evolution of linear and nonlinear angular velocities fluctuation.

Fig. 20 Relative teeth velocity fluctuation following teeth displacement in the case of two backlashes located on the first and second stages.

5 Conclusion

In this work, we are interested to the numerical simulation of nonlinear dynamic response of a one-stage and two-stage gears used in rotating machines. For this, we propose to apply the algorithm based on modeling using a new regularization technique of functions representing the contact, a Newmark time scheme, a change of variables from solution at the previous time, a homotopy with pre-conditioner and a Taylor series expansions. The aim of our work is to reduce the computation time estimated by number of inversions of the tangent matrix with respect to the iterative incremental Newton–Raphson method. To show the efficiency of the proposed algorithm, we have chosen two examples of gear with one-stage and two-stage and we have discussed the different numerical solutions of the problems with or without backlash. Similarly, a preliminary study has been made through our modeling to find the optimal truncation order of the polynomial approximation, used in the Taylor series expansions, which allows a higher quality of the solution of the posed problem.

References

All Tables

All Figures

	Fig. 1 Single stage spur gear system model.
In the text

	Fig. 2 Spring acting along the line of action.
In the text

	Fig. 3 Contact loss modelling.
In the text

	Fig. 4 Two-stage spur gear system model.
In the text

	Fig. 5 Spring acting along the line of action.
In the text

	Fig. 6 The approximated of the $({\tilde{F}}_s\delta (t))$ based on a regularization technique for different values of η.
In the text

	Fig. 7 Evolution of the solution quality log10||{R}|| versus the time t.
In the text

	Fig. 8 Influence of regularization parameter on the evolution of the transmission error δ versus the time t for η = 10−2, η = 10−3, η = 10−4, η = 10−5 and η = 10−6.
In the text

	Fig. 9 Time evolution of transmission error, Proposed algorithm truncated at order p = 3 and κ = 10−4 (ANM order 3), Newton–Raphson method with tolerance κ = 10−4 (Newton–Raphson).
In the text

	Fig. 10 Temporary evolution of the linear displacement of the first and second bearing, Proposed algorithm truncated at order p = 3 and κ = 10−4 (ANM order 3), Newton–Raphson method with κ = 10−4 (Newton–Raphson).
In the text

	Fig. 11 Decimal logarithm of the norm the residual vector versus time, Proposed algorithm with different truncation orders p = 3, 4, 5, 6, 7, 8 (ANM order p) and κ = 10−4, Newton–Raphson algorithm with κ = 10−4 (NR).
In the text

	Fig. 12 Temporal evolution of transmission errors, proposed algorithm with truncation order p = 3 (ANM) and κ = 10−4, Newton–Raphson algorithm with κ = 10−7 (NR) for (b1 = 0, b2 = 0) and (b1 = 10−4, b2 = 0).
In the text

	Fig. 13 Temporal evolution of linear and nonlinear angular velocities fluctuation.
In the text

	Fig. 14 Temporal evolution of bearings ẋ1 and ẋ3.
In the text

	Fig. 15 Temporal evolution of bearings ẍ1 and ẍ3.
In the text

	Fig. 16 Relative teeth velocity fluctuation following teeth displacement in the case of one backlash located on the first stage.
In the text

	Fig. 17 Decimal logarithm of the norm the residual vector versus the time, Proposed algorithm with truncation order p = 3 (ANM order 3) and κ = 10−4, Newton–Raphson algorithm with κ = 10−4 (NR).
In the text

	Fig. 18 Temporal evolution of transmission errors, proposed algorithm with truncation order p = 3 (ANM) and κ = 10−4, Newton–Raphson algorithm with κ = 10−7 (NR) for (b1 = 0, b2 = 0) and (b1 = 10−4, b2 = 10−4).
In the text

	Fig. 19 Temporal evolution of linear and nonlinear angular velocities fluctuation.
In the text

	Fig. 20 Relative teeth velocity fluctuation following teeth displacement in the case of two backlashes located on the first and second stages.
In the text