Static and dynamic models for spiral bevel gears

– This paper synthesizes a whole of work carried out on the spiral bevel gears from quasi-static and dynamic models viewpoint. A sophisticated quasi-static model makes it possible to calculate the tooth loads, the pressures, the instant mesh stiﬀness, and the deﬂections on the ﬂanks of spiral bevel gears. Based on these results, two three-dimensional lumped parameter dynamic models are presented. The mechanical system under consideration comprises: a spiral bevel pinion and gear connected by a time-varying non-linear mesh stiﬀness function, and mounted on two shafts simulated by Timoshenko’s beams supported by bearings. Two variants are considered which rely on diﬀerent contact stiﬀness simulations: (a) using an averaged mesh stiﬀness function acting at the centroid of the loaded areas on tooth ﬂanks and (b) a more local approach based on a discrete distribution of the local mesh stiﬀness, elements over the contact areas. A number of results are presented and commented which illustrate the interest of these dynamic models.


Introduction
The dynamic behaviour of gears has been extensively studied over last decades with two main objectives [1]: improve life time and reduce mesh noise.Because dynamic tooth loads can be significantly higher than those in quasistatic conditions, the life time of mechanical transmissions can be strongly affected.In this context, it is essential to predict the possible overloads on tooth contacts and the associated critical frequencies.From an acoustical point of view, gears are considered as potentially significant noise sources, mainly because of the vibration transfer from the a Corresponding author: michele.guingand@insa-lyon.fr meshing to the casing by the bearing-shaft assembly.In addition to experimental methods, theoretical analyses of the dynamic behaviour of the systems have been largely used.A number of dynamic models have been proposed for spiral bevel and hypoid gears [2][3][4].In particular, Lim and Cheng [4] have introduced a three-dimensional dynamic model of hypoid gears extended later to account for misalignments [5] and nonlinear time-varying mesh stiffness [6,7].In the continuation of these works, the influence analysis of four misalignment errors on the dynamic behaviour of spiral bevel gears has been performed by Peng and Lim [8].Li and Hu [2] have developed a 47 degree of freedom model for spiral bevel gears, based on their initial model of bevel gears [9].Finally, Gao et al. [3] have proposed a model dedicated to the study of shocks in spiral bevel gear systems.In this paper, two dynamic models are presented which make it possible to obtain instantaneous displacements and tooth loads over a broad speed range.In addition, the second model provides dynamic tooth loads and contact pressures on tooth flanks.The two models incorporate actual spiral bevel gear geometries and some of the results delivered by an accurate quasi-static contact model (ASLAN, software developed by LaMCoS of INSA of Lyon).The fundamentals of the quasi-static modeling are synthesized in the first section and it is shown how the instant mesh stiffness and kinematical error obtained using ASLAN are used as input data in the dynamic model.

Static model
The load sharing computation is decomposed as follows: definition of gear geometry by simulating the gear manufacturing process, simulation of no-load kinematics (taking into account misalignments), load calculations: contact pressure distribution, transmission error, stiffness. . .

Definition of the geometry
The proposed method enables the definition of spiral bevel gear geometry based on the Gleason type generation process.In this study, all the UMC (Universal Motion Concept) motions are not considered and a simplified case is analyzed.The method presented is based on Litvin works [10].

Tool characteristics
The envelopes of the tool cutters are considered as the union of a conical surface generating the active tooth surfaces and a toroidal surface cutting the tooth root shapes.
Once the tool characteristics are defined, the cutting process is simulated to obtain the pinion and gear geometry.The process is illustrated in Figure 2, where the different reference frameworks are defined as follows: -R m2 (x m2 , y m2 , z m2 ): linked to the machine, -R f (x f , y f , z f ): linked to the machine, -R s (x s , y s , z s ): linked to the tool, -R 2 (x 2 , y 2 , z 2 ): linked to the gear.

Gear positioning during the cutting process
Generation is kept limited to two generative movements which, in the cutting machine framework, correspond to the rotation of the cradle (tool support) and the rotation of the cut pinion.These two rotations are linked by a constant ratio: m P2 .Compared with the UMC model, the hypotheses are: the two gear flanks are manufactured simultaneously, the tool axis remains parallel to the principal machine tool axis (no tilt angle, no swivel), ratio m P2 is constant (no Modified Roll Motion), the pitch cone tip of the gear is on the principal axis of the cutting machine (no radial or axial offset).
From a general point of view, the equation of an envelope is obtained with the following form: where α o and β o are two parameters describing the location of a point on the tool whose envelope is sought, γ o is a tool positioning parameter and f e is an equation linking these parameters deduced from: e, m and o refer respectively to the reference axes linked to the envelope, cutting machine and tool.− → V and − → n correspond respectively to the velocity and normal vectors at any point.

Loaded analysis
A number of contributions on the calculation of tooth load sharing can be found in the literature.Sentoku [11] and Bruyere [12] considered bevel gears, when Gosselin [13], Icard [14], and Simon [15,16] studied spiral bevel and hypoid gears.In the majority of these works, approximate elastic models are employed which rely on curve-fitted experimental results and/or simplified tooth representations (for example, tooth structural deflections are deduced from a cantilever beam theory).The theory of Hertz is used to calculate the contact characteristics such as pressure distribution, contact area and deflections between the mating teeth.The main advantage of these models is that the computational times are limited and extensive parameter analyses can be performed; on the other hand, the accuracy of the simulations can be questionable in some cases.Alternative methods such as the Finite Prisms Method (Olakorede [17]) or the Finite Band Method (Gosselin [18]) have been proposed which are still effective from a computational time viewpoint but more precise in terms of deflections under load, although they are limited to standard cases.In this paper, the load sharing between the teeth is determined by combining the compatibility conditions in terms of displacements and the static moment balance.The procedure makes it possible to obtain the instantaneous contact pressures, the transmission error under load as well as the mesh stiffness and the root stresses.The compatibility conditions take into account, at the same time, the total bending and contact deflections which are both characterized by coefficients of influence.Structural and contact effects are separated (Icard [14]) the structural deflections of the parts are calculated by using standard Finite Element Model whereas the local contact compliance is computed via a local approach based on a discrete formulation of the results of Boussinesq [19] for elastic half-spaces.

Assumptions on the contact area
Since the pinion and the gear have the same elastic properties, they are quasi-identical solids in contact.For most of the practical applications, the contacts are fully lubricated and the coupling between the normal force and tangential displacements and between the tangential stresses and normal displacement can be neglected so that the zone of contact under load depends on normal loads only.
The potential contact area under load is supposed to be located in a plane parallel to the tangent plane at the unloaded contact point which is discretized in N rectangular cells of constant size on which the pressure distribution is considered as constant.

Algorithm of the load sharing computation
On the potential contact zone (taking into account several tooth pairs), load sharing is determined such that the compatibility conditions in terms of displacement (no inter-penetration of the parts) are satisfied, i.e.
outside the contact area (4) where U 1i and U 2i are the normal displacements of bodies 1 and 2 at point i, ei i is the initial gap at point i, α is the global normal approach of the contacting surfaces and p i the pressure at point i.The distance between the two bodies at point i after load is denoted y i .
Equations ( 3) and ( 4) can be synthesized as: Assuming that the pressure on a small rectangular cell of area S i is constant (S i is attached to the potential point of contact i and lies in the tangent plan), the global torque equilibrium leads to: where N is the number of rectangular surfaces considered in the model, p i is the pressure on the element i of surface S i and R i the associated lever arm.
The elasticity of the contacting solids is accounted for by using influence coefficients such that the displacement at point i reads: where C ij are the influence coefficients (effect of point j on point i).
The problem can then be formulated as follows: By combining all these equations and introducing: the pressure field can be calculated as: Based on this formulation, the system can be solved iteratively by using a Fixed Point Method which ensures a uniform, fast and accurate convergence.Parameters U 1i , U 2i , ei i , ef i and α are represented in Figure 4. Since the surfaces in contact are supposed to be planar and all parallel, convergence is achieved when the final gaps ef i at all the contact points are identical and equal to the normal approach α.

Definition of the influence coefficients
The matrix of influence coefficients takes into account both bending and contact compliances.However, these two effects are supposed to be independent so that three elemental matrices can be separated: one for contact deformations (C s ij ) and two for pinion (C Pf ij ) and gear (C Rf ij ) bending deflections respectively such that: A -Influence of bending coefficients The coefficients of influence for bending are calculated using a standard Finite Element Model because analytical formulations are not possible for complex rims or webs, particular assemblies, etc.However, full FE simulations including the modeling of rims, shafts, casings, etc., can be time consuming and/or require significant computational resources not compatible with parameter analyses for example.In order to significantly reduce computational times, the coefficients of influence are not determined for each kinematical position, neither for each potential contact points, but using the results at a limited number of points and for one position as described in [20] by Teixeira Alves.

B -Influence coefficients for contact
The pinion and gear teeth in contact are approximated by two half -elastic spaces and the displacements under normal forces are deduced by using the potential functions of Boussinesq [19].A general expression is derived which depends on the elastic constants of the two halfelastic spaces, on the dimensions of discrete cells, and on Pinion Gear Tangent plane the coordinates of the different points of contact in the tangent plan.

Unloaded gaps
In order to deal with the conditions of compatibility (Eqs.( 3) and ( 4)), the initial separations ei i between the mating surfaces must be known for each tooth pair potentially in contact.Projecting the meshing points (for example point I in Fig. 5) in the normal direction with respect to the tangent plane and denoting P the projection on the pinion surface and R that on the gear surface, the initial separation of gap at that point is expressed as the distance RP .
Once all the initial separations for all the possible tooth pairs in contact are known, the computation of the load sharing can be undertaken.As specified earlier in the paper, the fixed point method is used to solve the equations of displacement compatibility which leads to pressure distributions, transmission errors and the global mesh stiffness.

Mounting the gear and pinion
Assembly errors can be introduced which are linked to the change of basis described below (Figs. 6 and 7).
− → r a and − → r b represent the coordinates of a point in the Ra (O a , x a , y a , z a ) pinion reference system and Rb (O b , x b , y b , z b ) gear reference system respectively.P , G and E are the axial error for the pinion, for the gear, and the offset.
is the shaft angle.

Examples
In this section, two kinds of results are presented that can be obtained with the numerical model developed in this part: load sharing and instantaneous contact pressures.Other results, such as contact patterns under load, transmission errors or mesh stiffness are not shown here.All these results cannot be compared with those find in the literature since the simulated cutting process is simplified.

A -Load sharing
The first result deals with load sharing (Fig. 8).The number of teeth in contact is obtained for each kinematic position, as well as the global load supported by each tooth pair in contact.The total force is also shown.

B -Instantaneous contact pressures
Figure 9 presents, for a given kinematic position, the contacts on each loaded tooth flank (the study is achieved for 5 teeth potentially in contact), for both pinion and gear teeth.The pressure on the mean contact line is also indicated for all the teeth.

Dynamic models
Based on the results of the quasi-static model presented in Section 2, two dynamic models are developed and presented in this section.

Model
The model of spiral bevel gear is presented in Figure 10.The system of coordinates R(O, X, Y, Z) is rooted at O placed at the apex of the pitch cone of the pinion.The origins for the local coordinate systems attached to the pinion and the gear are their respective centers.The pinion-gear pair is discretized into 6 nodes denoted E, 1, 2, 3, 4, S (Fig. 10) respectively with six degrees of freedom each: w S , ϕ S , ψ S for bending displacements and rotations and finally, θ E , θ 1 , θ 2 , θ 3 , θ 4 , θ S for torsion.
The pinion and the gear are assimilated to rigid bodies and a simple Wrinckler foundation model is used for the tooth contacts [21].By so doing, the deflections at any potential point of contact Mi can be expressed in terms of the degrees of freedom at the centers of the pinion and the gear [22].The local mesh forces are derived by multiplying the local stiffness and the mesh deflection.Lagrange's equations lead to a non-linear time-varying mesh stiffness matrix and a forcing term vector generated by the initial separations between the mating flanks.The shafts are simulated by two-node finite elements, including shear secondary effects, whose mass and stiffness matrices are conventional [23].

Equations of motion
After assembling all the elementary mass and stiffness matrices, the equations of motion of the system are derived from Lagrange's equations: The inertial kinetic energy associated with the spiral bevel gear is given by Equation ( 13) and it generates a constant mass matrix [M E ] and a time-dependent forcing term [F 2 ] due to the inertial effects when rotational speeds are not constant (produced by the no-loaded transmission error).
where: I 1 , I 4 : moments of inertia of the pinion and the gear J 1 , J 4 : polar moments of inertia of the pinion and the gear.

Global model
In a first version of the dynamic model (referred to as model 1), potential energy is calculated by assuming that mesh stiffness can be reduced to a unique stiffness element in the normal direction and located at the centroid of the contact area.The deflection at this point is equal to the normal approach with respect to rigid-body positions and can be expressed by projecting the contributions of the degrees-of-freedom in the normal direction as: from which the mesh strain energy can be expressed under the form: where: projection (or structural) vector in the normal direction, j: location in the meshing, B m : centre of the contact area, n 1 , n 2 : direction of the action lines of the meshing loads.
The global mesh stiffness is given by the numerical code named ASLAN and it is estimated as: where F is the global mesh force and α is the normal approach.

Discretized model
For the second version of the model (model 2), a discrete distribution of the local stiffness elements on contact areas is taken into account and it leads to the following formulation: where A is the contact area and i the points in the contact area.
In this case, it is to be noted that an additional forcing term F 1 appears which is caused by the initial separations between the tooth flanks (as opposed to what is obtained using model 1).The elemental stiffness elements k i are approximated by dividing the elemental load at point Mi by the local deflection at the same point as: where ef i and ei i are the final and initial gaps at point i delivered by ASLAN with the method presented in Section 2 (see Fig. 4 for their definitions).

Energy dissipation
In what follows, energy dissipation is accounted for by introducing a global Rayleigh viscous damping matrix [C] (a linear combination of the global mass and the averaged stiffness matrices).
The equations of motion resulting of Equations ( 14), (15) and/or (17) point to a parametrically excited nonlinear differential linear system of the form: where q represents the global generalized displacements vector.
The equations are solved by combining a time-step integration scheme which is combined with a unilateral contact algorithm as described in [22].

Dynamical loads
Considering model 1, Figure 11 shows an example of dynamic tooth load time variation (the abscissa represents the number of mesh periods) for a pinion speed of 17 200 rpm.The numerical transients, observed for the first mesh periods, are due to the initial conditions (here, the static solution with averaged mesh stiffness) which can be substantially different from the actual dynamic solution.However, it is observed that steady state conditions are obtained very rapidly.The corresponding spectrum (calculated by FFT once steady motion is established) exhibits the classic discrete peaks associated with the mesh frequency and its harmonics.

Pressure distribution
By using model 2, the spatial distribution of mesh forces on the flanks of the teeth can be estimated (Fig. 12).At low speed (1 rpm on the pinion), the results are expressed in terms of contact pressure and compare very well with the results given by ASLAN in quasistatic conditions.It is to be noted in this example, that the pressure peaks are located in the central part of the tooth flank.The dynamic contact patterns at higher speed (17 600 rpm) are represented in Figure 13 where significant fluctuations of the pressure field can be observed.Compared with the classic quasi-static analyses in the literature, the proposed model makes it possible to have access to dynamic pressures (and stresses) thus analyzing more precisely the actual behaviour of high-speed spiral bevel gears.

Dynamics factors
The preceding results can be generalized over a broad range of speeds by introducing the dynamic tooth load factor defined as: F d , F s : dynamic and static loads.The maximum dynamic tooth load results of the integration of the dynamic pressures, and care was taken to make sure that the maximum value at every speed was sought after the system had reached steady state conditions.The results in Figure 14 reveal a structure of dynamic response close to that obtained for spur and helical gears.For both models 1 and 2, a major tooth critical speed emerges around 17 000-18 000 rpm on the pinion whereas secondary response peaks excited by the higher harmonics of the mesh frequency are present in the lower speed region.Some differences can be reported between the results of models 1 and 2: the amplifications are not the same and the critical frequencies are slightly shifted.In both cases, significant dynamic effects are observed (maximum DF between 1.5 and 1.6) which certainly reduce the reliability of the system and must be integrated in any realistic strength analysis.

Conclusion
The definition of the geometry of spiral bevel gears has been presented along with a static model aimed at calculating tooth load distributions and transmission errors.Two lumped parameter dynamic models with several degrees of freedom have been then introduced, which use accurate descriptions of the pinion and gear geometries but employ a different mesh stiffness model.The first model is similar to the classic approaches in the literature, with a single time-varying mesh stiffness connecting the pinion and the gear.The second model relies on a local   description of the instantaneous contact conditions via a distribution of mesh stiffness elements (Wrinckler's elastic foundation) over the potential contact area between the pinion and the gear.One of the advantages is that dynamic pressure distributions on the flanks can be estimated.Whatever the model used, the dynamic curve response presents similar trends with i) significant tooth force amplifications at the major critical speed and ii) several secondary resonances generated by the harmonics of the mesh excitations.
In order to validate the numerical results presented in this article, a test bench is currently realized.

Fig. 1 .
Fig. 1.Envelope definition of the cutter for the numerical model.

Figure 3
Figure 3 corresponds to a pinion and gear set obtained with the developed numerical model, based on the simulation of a simplified Gleason type generation process.

Fig. 5 .
Fig. 5. Projections of the meshing points on the pinion and the gear.

Fig. 9 .
Fig. 9. Contact lines and instantaneous pressures for a given kinematic position.

Fig. 13 .
Fig. 13.Dynamic contact pattern at high speed in 3 dimensions (A), and projected on the tooth flank (B).

Fig. 14 .
Fig. 14.Dynamic coefficients by the two models versus pinion rotational speed.

Table 1 .
Geometrical data of the spiral-bevel gear.

Table 2 .
Geometrical data of the shafts.