Issue 
Mechanics & Industry
Volume 19, Number 3, 2018



Article Number  310  
Number of page(s)  16  
DOI  https://doi.org/10.1051/meca/2017006  
Published online  11 October 2018 
Regular Article
An implicit algorithm for the dynamic study of nonlinear vibration of spur gear system with backlash
Laboratoire d’Ingénierie et Matériaux LIMAT, Faculté des Sciences Ben M’Sik, Université Hassan II de Casablanca,
Sidi Othman,
Casablanca, Morocco
^{*} email: b.braikat@gmail.com
Received:
22
February
2016
Accepted:
17
January
2017
In this work, we propose some regularization techniques to adapt the implicit high order algorithm based on the coupling of the asymptotic numerical methods (ANM) (Cochelin et al., Méthode Asymptotique Numérique, HermèsLavoisier, Paris, 2007; Mottaqui et al., Comput. Methods Appl. Mech. Eng. 199 (2010) 1701–1709; Mottaqui et al., Math. Model. Nat. Phenom. 5 (2010) 16–22) and the implicit Newmark scheme for solving the nonlinear problem of dynamic model of a twostage spur gear system with backlash. The regularization technique is used to overcome the numerical difficulties of singularities existing in the considered problem as in the contact problems (Abichou et al., Comput. Methods Appl. Mech. Eng. 191 (2002) 5795–5810; Aggoune et al., J. Comput. Appl. Math. 168 (2004) 1–9). This algorithm combines a time discretization technique, a homotopy method, Taylor series expansions technique and a continuation method. The performance and effectiveness of this algorithm will be illustrated on two examples of onestage and twostage gears with spur teeth. The obtained results are compared with those obtained by the Newton–Raphson method coupled with the implicit Newmark scheme.
Key words: spur gear system / high order algorithm / homotopy / vibration / backlash
© 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 nondifferentiable 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 piecewise 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 timevarying dynamic model for rightangle gear pair systems, taking into account both backlash and asymmetric mesh effects. The resolution is done using an enhanced multiterm harmonic balance method with a modified discrete Fourier transform process and the numerical continuation method. They studied the stability of the steadystate 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 twostage 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, superharmonic and subharmonic 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 onestage and twostage 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 onestage 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_{1}. The second one is composed of the wheel and the loading machine connected by the output shaft and has a mass m_{2}. Motor, pinion, wheel and load are considered as rigid bodies, their respective inertias are I_{m}, I_{1}, I_{2} and I_{L}. The two parts are supported by bearings modelled by linear springs , , and . Connecting shafts are modelled by torsional stiffness and . 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: (1)
Mesh stiffness is modelled by a function 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 is given as a function of an average component denoted by and of a fluctuating component denoted by . This fluctuating component is given by: (2) where n is an integer which designates the number of the period, i designates the stage number, the terms is the contact ratio at stage i, is the period and the average component , in the case of twostage; with Z_{2} and Z_{3} 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]: (3)where α is the pressure angle, and are respectively the base radius of pinion and wheel, {χ} is defined by: (4)and, here, the degrees of freedom θ_{1} and θ_{2} 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: (5)where 2b_{s} is the total backlash as shown in Figure 3.
Writing Lagrange formulation allows obtaining the equation of motion expressed by: (6) where [M] is the inertia matrix given by: (7)
[K_{b}] is the bearings and shafts matrix expressed by: (8)
The vector of external applied forces {F_{ext}} is expressed by: (9)
The damping matrix [C] in this model is proportional both to the mass and stiffness matrices of the system given by: (10)
where c_{1} and c_{2} are the damping coefficients which are expressed respectively in s^{−1} 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 .
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 twostage spur gear system with backlash
A model of a twostage 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_{1}) along x (horizontal) and y (vertical) directions. The corresponding DOF are x_{1} and y_{1}.

Translations of part 2 (having the mass m_{2}) x and y directions. The corresponding DOF are x_{2} and y_{2}.

Translations of part 3 (having the mass m_{3}) along x and y directions. The corresponding DOF are x_{3} and y_{3}.

Rotations of motor θ_{m}, pinion 1 θ_{1}, wheel 1 θ_{2}, pinion 2 θ_{3}, wheel 2 θ_{4} and load θ_{L}.
The following inertia are considered: I_{m} for the motor, I_{1} for the pinion 1, I_{2} for wheel 1, I_{3} for pinion 2, I_{4} for wheel 2 and I_{L} for load. Input, output and intermediate shafts elasticity is modelled by torsional stiffness respectively , and . Bearings supporting input shaft are modelled by linear springs and acting along x and y axis. Bearings supporting intermediate shaft are modelled by linear springs and . Bearings supporting intermediate shaft are modelled by linear spring and .
Two mesh stiffness functions and 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 twostage, two mesh forces should be considered and . 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: (11)
Two transmission errors δ_{1}(t) and δ_{2}(t) for the twostage are expressed by [16]: (12) where α_{1} and α_{2} are the pressure angles, {χ_{1}} and {χ_{2}} are given by: (13)and, here, the degrees of freedom θ_{1}, θ_{2}, θ_{3} and θ_{4} are positive quantities. Taking into account of the backlashes in the two mesh zones, the expression of the mesh forces and respectively measured on the meshes in stage 1 and stage 2 are as following: (14) (15)
By the same Lagrange formulation it is possible to obtain the equation of motion of the twostage gear system by: (16) where [M] is the inertia matrix given by: (17)
[K_{b}] is the bearing and shafts matrix expressed by: (18)
The vector of external applied forces{F_{ext}} is expressed by: (19)
The initial conditions used for this problem are: {q} = {0} and . So, the initial conditions of transmission errors deduced from equations (3) and (12) are given by: (20)
Fig. 4 Twostage 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), nonsmooth functions appear F_{s}(δ(t)) defined by (5), defined by (14) and defined by (15). The idea is to replace these non smooth functions by a smooth function denoted by: (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: (22)
Taking into account the equation (22), the modified function (21) is set in the following form: (23)
Fig. 6 The approximated of the 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 and acceleration at time (n + 1)Δt as follows: (24) where , , , , , n denotes the time level, Δt is the time step, the constants β and γ have the range 0 ≤ β ≤ 1, with typical values being , . Using equation (24), the problem (6) at time (n + 1)Δt is written as: (25)
where {q^{n+1}}, δ^{n+1}, h^{n+1} and g^{n+1} 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^{n+1}}, δ^{n+1}, h^{n+1} and g^{n+1} one seeks the perturbations {Δq}, ΔD, ΔH and ΔG from the solutions {q^{n}}, δ^{n}, h^{n} and g^{n} at time t = nΔt in the following manner. (26)
The expressions (26) are introduced into equations (3), (22), and (25), we obtain a nonlinear problem: (27) where is the tangent matrix given by (28)with (29)
{FQ} is a quadratic form defined by: (30) and {S^{n}} is the second hand side expressed by: (31)
3.4 Homotopy transformation
It is well known that the homotopy transformations have been developed to overcome the local convergence of the Newtonlike 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]: (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: (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: (34) (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: (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 onestage and twostage. The two nonlinear problems are solved by the detailed algorithm above. In both examples, the matrix [K^{*}] has the same form as defined at the initial conditions. If , 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: onestage 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 preconditioner [K^{*}] which is taken equal to the tangent matrix evaluated at the initial time, a time step Δt = 3.7 × 10^{−6} s, a tolerance parameter κ = 10^{−4} and various values for the truncation order for our proposed algorithm, a tolerance parameter κ = 10^{−4} for Newton–Raphson algorithm and a regularization parameter η = 10^{−6}. 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: (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^{−2}, η = 10^{−3}, η = 10^{−4}, η = 10^{−5} and η = 10^{−6}. According to this test, we noticed that the solution of the problem is stabilized starting from the value η = 10^{−3} (see Fig. 8).
In the following, we fix a truncation order p = 3, a tolerance κ = 10^{−4} for our algorithm and a tolerance κ = 10^{−4} for the Newton–Raphson method for a quality solution less than 10^{−4} and a regularization parameter η = 10^{−3}. 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.
Parameters of the onestage spur gear transmission.
Fig. 7 Evolution of the solution quality log_{10}{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). 
Comparison between the number of triangulation tangent matrices for the both methods.
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: twostage gear system
In this example, the nonlinear behaviour of a twostage 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 twostage gear system.
Parameters of the twostage spur gear transmission.
4.2.1 Case of one backlash at the first stage: b_{1} = b/2 and b_{2} = 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^{−6} s. For this, we fix the tolerance parameter of the two algorithms, i.e. Newton–Raphson algorithm and the proposed algorithm to 10^{−4} and η_{1} = η_{2} = 10^{−6}. 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^{−4} for the proposed algorithm and a tolerance parameter κ = 10^{−4} for the Newton–Raphson method.
In Figure 12a and b, we represent the time evolution of transmission errors δ_{1} and δ_{2} 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_{1} = 0, b_{2} = 0) and nonlinear (b_{1} = b/2, b_{2} = 0) speeds ẋ_{1} and ẋ_{3} of the first and third block along the axis x and the corresponding accelerations ẍ_{1} and ẍ_{3} 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 (b_{1} = 0, b_{2} = 0) and (b_{1} = 10^{−4}, b_{2} = 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. 
Comparison between the number of triangulation tangent matrices for the both methods.
4.2.2 Case of two backlashes at the first and second stages: b_{1} = b/2 and b_{2} = b/2
In this second part, we take into account the presence of two backlashes at both stages (b_{1} = b/2, b_{2} = 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^{−4} 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 and , 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 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 (b_{1} = 0, b_{2} = 0) and (b_{1} = 10^{−4}, b_{2} = 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. 
Comparison between the number of triangulation tangent matrices for the both methods.
5 Conclusion
In this work, we are interested to the numerical simulation of nonlinear dynamic response of a onestage and twostage 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 preconditioner 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 onestage and twostage 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
 A. Kahraman, R. Singh, Nonlinear dynamics of a spur gear pair, J. Sound Vib. 142 (1990) 49–75 [CrossRef] [Google Scholar]
 A. Kahraman, R. Singh, Interactions between timevarying mesh stiffness and clearance nonlinearities in a geared system, J. Sound Vib. 146 (1991) 135–156 [CrossRef] [Google Scholar]
 J. Lin, R.G. Parker, Mesh stiffness variation instabilities in twostage gear systems, J. Vib. Acoust. 124 (2002) 68–76 [CrossRef] [Google Scholar]
 F. Chaari, W. Bartelmus, R. Zimroz, T. Fakhfakh, M. Haddar, Effect of load shape in cyclic load variation on dynamic behavior of spur gear system, Key Eng. Mater. 518 (2012) 119–126 [CrossRef] [Google Scholar]
 F. Chaari, R. Zimroz, W. Bartelmus, T. Fakhfakh, M. Haddar, Model based investigation on a two stages gearbox dynamics under nonstationary operations, in: Condition monitoring of machinery in nonstationary operations, Springer, Berlin, Heidelberg, 2012, pp. 133–142 [CrossRef] [Google Scholar]
 R.W. Gregory, S.L. Harris, R.G. Munro, Dynamic behaviour of spur gears, Proc. Inst. Mech. Eng. 178 (1963–1964) 207–218 [CrossRef] [Google Scholar]
 S. Yongjun, Y. Shaopu, P. Cunzhi, L. Xiandong, Nonlinear dynamics of a spur gear pair with timevarying stiffness and backlash, J. Low Freq. Noise Vib. Act. Control 23 (2004) 179–187 [CrossRef] [Google Scholar]
 R.G. Parker, S.M. Vijayakar, T. Imajo, Nonlinear dynamic response of a spur gear pair: modelling and experimental comparisons, J. Sound Vib. 237 (2000) 435–455 [CrossRef] [Google Scholar]
 H. Moradi, H. Salarieh, Analysis of nonlinear oscillations in spur gear pairs with approximated modelling of backlash nonlinearity, Mech. Mach. Theory 51 (2012) 14–31 [CrossRef] [Google Scholar]
 R.G. Munro, L. Morrish, D. Palmer, Gear transmission error outside the normal path of contact due to corner and top contact, J. Mech. Eng. Sci. 213 (1999) 389–400 [CrossRef] [Google Scholar]
 A. Kubo, K. Yamada, T. Aida, S. Sato, Research on ultra high speed gear devices, Trans. Jpn. Soc. Mech. Eng. 38 (1972) 2692–2715 [CrossRef] [Google Scholar]
 A. Kahraman, G.W. Blankenship, Effect of involute contact ratio on spur gear dynamics, J. Mech. Des. 1 (2007) 112–118 [Google Scholar]
 S.M. Wang, Analysis of nonlinear transient motion of a geared torsional, J. Eng. Ind. 96 (1974) 51–59 [CrossRef] [Google Scholar]
 Y. Cai, T. Hayashi, The linear approximated equation of vibration of a pair of spur gears: theory and experiment, J. Mech. Des. 116 (1994) 558–570 [CrossRef] [Google Scholar]
 Y. Junyi, P. Tao, C.L. Teik, An enhanced multiterm harmonic balance solution for nonlinear periodβ dynamic motions in rightangle gear pairs, Nonlinear Dyn. 76 (2014) 1237–1252 [CrossRef] [Google Scholar]
 L. Walha, T. Fakhfakh, M. Haddar, Nonlinear dynamics of a twostage gear system with mesh stiffness fluctuation, bearing flexibility and backlash, Mech. Mach. Theory 44 (2009) 1058–1069 [CrossRef] [Google Scholar]
 Z. El Ouehabi, F. Chaari, H. Lahmam, B. Braikat, H. Mottaqui, T. Fakhfakh, M. Haddar, N. Damil, Un algorithme implicite d'ordre eleve pour la dynamique des systemes d'engrenages utilises dans les machines tournantes, Rev. Mec. Appl. Theor. 2 (2011) 455–464 [Google Scholar]
 M. Jamal, B. Braikat, S. Boutmir, N. Damil, M. PotierFerry, A high order implicit algorithm for solving instationary nonlinear problems, Comput. Mech. 28 (2002) 375–380 [CrossRef] [Google Scholar]
 S. Boutmir, B. Braikat, M. Jamal, N. Damil, B. Cochelin, M. PotierFerry, Des solveurs implicites d'ordre superieur pour les problemes de dynamique non lineaire des structures, Rev. Eur. Elem. Finis 13 (2004) 449–460 [Google Scholar]
 A. Timesli, B. Braikat, H. Lahmam, H. Zahrouni, A new algorithm based on moving least square method to simulate material mixing in friction stir welding, Eng. Anal. Bound. Elem. 50 (2015) 372–380 [CrossRef] [Google Scholar]
 E.L. Allgower, K. Georg, Numerical continuation methods: an introduction, in: Springer series in computational mathematics, vol. 13, 1990 [CrossRef] [Google Scholar]
 M.A. Crisfield, Faster modified Newton–Raphson iteration, Comput. Methods Appl. Mech. Eng. 20 (1979) 267–278 [CrossRef] [Google Scholar]
 T. Fakhfakh, F. Chaari, M. Haddar, Numerical and experimental analysis of a gear system with teeth defects, Int. J. Adv. Manuf. Technol. 25 (2005) 542–550 [CrossRef] [Google Scholar]
 F. Chaari, W. Baccar, M.S. Abbes, M. Haddar, Effect of spalling or tooth breakage on gearmesh stiffness and dynamic response of a onestage spur gear transmission, Eur. J. Mech. A/Solids 27 (2008) 691–705 [CrossRef] [Google Scholar]
 A.K. Chopra, Dynamics of structures: theory and applications to earthquake engineering, Prentice Hall, New Jersey, 2001 [Google Scholar]
 N.M. Newmark, A method of computation for structural dynamics, in: Proceedings of ASCE, vol. 85, 1959 [Google Scholar]
Cite this article as: Y. Hilali, B. Braikat, H. Lahmam, N. Damil, An implicit algorithm for the dynamic study of nonlinear vibration of spur gear system with backlash, Mechanics & Industry 19, 310 (2018)
All Tables
Comparison between the number of triangulation tangent matrices for the both methods.
Comparison between the number of triangulation tangent matrices for the both methods.
Comparison between the number of triangulation tangent matrices for the both methods.
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 Twostage 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 based on a regularization technique for different values of η. 

In the text 
Fig. 7 Evolution of the solution quality log_{10}{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 (b_{1} = 0, b_{2} = 0) and (b_{1} = 10^{−4}, b_{2} = 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 (b_{1} = 0, b_{2} = 0) and (b_{1} = 10^{−4}, b_{2} = 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.