An implicit algorithm for the dynamic study of nonlinear vibration of spur gear system with backlash

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., Methode Asymptotique Numerique, Hermes-Lavoisier, 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 non-linear problem of dynamic model of a two-stage 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 one-stage and two-stage gears with spur teeth. The obtained results are compared with those obtained by the Newton–Raphson method coupled with the implicit Newmark scheme.


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][2][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][8][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 equal 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][19][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 e varying from zero to one [18][19][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][19][20]. To illustrate the performance of the proposed algorithm, we give numerical comparisons with the classical implicit iterative algorithm of Newton-Raphson type [22].

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 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 k x 1 , k y 1 , k x 2 and k y 2 . Connecting shafts are modelled by torsional stiffness k u 1 and k u 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: Mesh stiffness is modelled by a function k m i ðtÞ spring acting along the line of action (see Fig. 2). This model has been used by Fakhfakh et al. [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 ðtÞ is given as a function of an average component denoted by k m i and of a fluctuating component denoted by k v i ðtÞ. This fluctuating component is given by: see equation (2) below where n is an integer which designates the number of the period, i designates the stage number, the terms e a i is the contact ratio at stage i, T e i is the period and the average component k m i ¼ k max ðe a i À 1Þ þ k min ð2 À e a i Þ, T e 2 ¼ Z 2 Z 3 T e 1 in the case of two-stage; 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]: where a is the pressure angle, r b 1 and r b 2 are respectively the base radius of pinion and wheel, {x} is defined by: T fxg ¼ 〈sinðaÞ; cosðaÞ; 0;r b 1 ;ÀsinðaÞ;ÀcosðaÞ;r b 2 ; 0〉 ð4Þ and, here, the degrees of freedom u 1 and u 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: where 2b s is the total backlash as shown in Figure 3. Writing Lagrange formulation allows obtaining the equation of motion expressed by: where [M] is the inertia matrix given by: [K b ] is the bearings and shafts matrix expressed by: if ðn þ 1 À 1ÞT e 1 t ðn þ 1ÞT e 1 for one-stage if nT e i t ðn þ 2 À i ÞT e i for two-stage ð2Þ The vector of external applied forces {F ext } is expressed by: The damping matrix [C] in this model is proportional both to the mass and stiffness matrices of the system given by: 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 ḟqg ¼ f0g.

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 1 ) along x (horizontal) and y (vertical) directions. The corresponding DOF are x 1 and y 1 . 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 k u 1 , k u 2 and k u 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 ðtÞ and k m 2 ðtÞ 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: Two transmission errors d 1 (t) and d 2 (t) for the twostage are expressed by [16]: where a 1 and a 2 are the pressure angles, {x 1 } and {x 2 } are given by: T fx 1 g ¼< sin ða 1 Þ; cos ða 1 Þ; 0; r b 1 ; À sin ða 1 Þ; À cos ða 1 Þ; Àr b 2 ; 0; 0; 0; 0; 0 > T fx 2 g ¼< 0; 0; 0; 0; sin ða 2 Þ; À cos ða 2 Þ; 0; r b 3 ; À sin ða 2 Þ; and, here, the degrees of freedom u 1 , u 2 , u 3 and u 4 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: By the same Lagrange formulation it is possible to obtain the equation of motion of the two-stage gear system by: where [M] is the inertia matrix given by: see equation (17)   [K b ] is the bearing and shafts matrix expressed by: see equation (18) below The vector of external applied forces{F ext } is expressed by: The initial conditions used for this problem are: {q} = {0} and ḟqg ¼ f0g. So, the initial conditions of transmission errors deduced from equations (3) and (12) are given by: 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.

Regularization technique
In the case of problems (6) and (16), non-smooth functions appear F s (d(t)) defined by (5), F s 1 d 1 ðtÞ defined by (14) and F s 2 d 2 ðtÞ defined by (15). The idea is to replace these non smooth functions by a smooth function denoted by:F The function defined here involves a regularization parameter h. 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(d(t)) and g(d(t)) defined by: Taking into account the equation (22), the modified function (21) is set in the following form:

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 ḟq nþ1 g and acceleration fq nþ1 g at time (n + 1)Dt as follows: f _ q nþ1 g ¼ a 0 ðfq nþ1 g À fq n gÞ À fc n g f€ q nþ1 g ¼ Dtba 0 ðfq nþ1 g À fq n gÞ þ fw n g ; where fv n g ¼ a 1 f _ q n g þ a 2 f€ q n g;fc n g ¼ 1 À a 1 Dtb ð Þ f_ q n gþ Using equation (24), the problem (6) at time (n + 1)Dt is written as: where {q n+1 }, d n+1 , h n+1 and g n+1 are the unknowns at time t = (n + 1)Dt.

Introduction of a perturbation from the previous time
Instead of searching the unknowns {q n+1 }, d n+1 , h n+1 and g n+1 one seeks the perturbations {Dq}, DD, DH and DG from the solutions {q n }, d n , h n and g n at time t = nDt in the following manner.
The expressions (26) are introduced into equations (3), (22), and (25), we obtain a nonlinear problem: where ½K n T is the tangent matrix given by {FQ} is a quadratic form defined by: and {S n } is the second hand side expressed by: fS n g ¼ fF ext g þ ½Mfc n g À ½Cfv n g À ½K b fq n g À ðk nþ1 m i d n þ

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) and the parameter e of this transformation is considered as an expanding parameter [18][19][20]: ½K Ã fDW g þ eðð½K n T À ½K Ã ÞfDW g þ fFQgÞ ¼ efS n g ð32Þ By this way, the solutions DW(e) of the nonlinear equation (32) passe continuously from 0 for e = 0 to the solution of the nonlinear problem (27) for e = 1.

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: Pressure angle a 1 = a 2 = 20°D amping coefficients (1/s), (s) c 1 = 10 À2 , c 2 = 10 À6 Mesh frequencies (Hz) By introducing the expansions (33) into equation (32) and by equating like powers of e, we obtain a set of linear problems satisfied by the terms of the series (33) given by: 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 k. The step length is then computed by p ≥ 2 : defining a maximal value of the parameter e max [18][19][20] as follows: The solution of problem (25) is obtained while the parameter e max ≥ 1.

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 ½K n T defined at the initial conditions. If K Ã ¼ K n T , 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.

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 ½K n T evaluated at the initial time, a time step Dt = 3.7 Â 10 À6 s, a tolerance parameter k = 10 À4 and various values for the truncation order for our proposed algorithm, a tolerance parameter k = 10 À4 for Newton-Raphson algorithm and a regularization parameter h = 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: 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 h on the solution of the problem. In Figure 8, we represent the transmission error versus to time for h = 10 À2 , h = 10 À3 , h = 10 À4 , h = 10 À5 and h = 10 À6 . According to this test, we noticed that the solution of the problem is stabilized starting from the value h = 10 À3 (see Fig. 8).  In the following, we fix a truncation order p = 3, a tolerance k = 10 À4 for our algorithm and a tolerance k = 10 À4 for the Newton-Raphson method for a quality solution less than 10 À4 and a regularization parameter h = 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.

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. 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 h 1 = h 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 k = 10 À4 for the proposed algorithm and a tolerance parameter k = 10 À4 for the Newton-Raphson method.
In Figure 12a and b, we represent the time evolution of transmission errors d 1 and d 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. 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 k = 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 _ u 1 and _ u 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 _ u 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.

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 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 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.