INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING, VOL. 40, 357–381 (1997) THREE-DIMENSIONAL GEARS MODELLING IN MULTIBODY SYSTEMS ANALYSIS ALBERTO CARDONA∗ Computational Mechanics Laboratory - INTEC, Universidad Nacional del Litoral - Conicet, Guemes 3450-3000 Santa Fe-Argentina SUMMARY We present a formulation for describing gear pairs in the three-dimensional analysis of exible mechanisms. The set of holonomic and non-holonomic constraint equations that denes the behaviour of gears is developed. The formulation is capable of representing almost any kind of gears used in industry: spur gears, bevel gears, hypoid gears, etc. All reaction forces due to gear engagement are accounted for. Coulomb friction forces due to relative sliding motion are also considered. Gear trains and planetary systems can be easily modelled by superposing several gear pairs in the same system. KEY WORDS: non-linear dynamics; multibody dynamics; nite rotations; gears 1. INTRODUCTION Gears are devices of widespread use for the transmission of rotary motion from one shaft to another, which are present in almost every machine. A pair of gears with non-intersecting and non-parallel axes develops three-dimensional gearing relative motion. All kinds of gears used in industry can be seen as particular cases of three-dimensional gearing: e.g. spur gears, bevel gears, hypoid gears, worm gears, etc. The literature on gears is vast and many classic texts treat them thoroughly (see, for instance, Reference 1 and the references therein). Several studies on gear modelling have been published in recent years, most of them with the objective of improving the analysis capabilities to allow the design of smaller, lighter and quieter transmissions. Rama Mohana Rao and Muthuveerappan2 evaluated root stresses for dierent positions of the contact line when it moves from the root to the tip in helical gear teeth using three-dimensional nite elements. Kahraman3; 4 developed a model to simulate the dynamic behaviour of a single-stage planetary gear train with helical gears. The formulation allows to model planetary gear sets with any number of planets. The model was used to predict the eect of manufacturing and assembly errors, and of teeth separations and mesh stiness uctuations on the load sharing between pinions.4 Geradin and Robert5 developed a general program to analyse torsional vibrations of gear trains. The model included eects as teeth exibility, manufacturing errors (eccentricity), and stiness uctuation. The model of the gear train could be composed of parallel spur and helical gears, conical bevel gears and planetary systems. ∗ Professor, Universidad Nacional del Litoral, and research sta member of Consejo Nacional de Investigaciones Cientcas y Tecnicas de la Republica Argentina CCC 0029–5981/97/020357–25 ? 1997 by John Wiley & Sons, Ltd. Received 25 April 1995 Revised 28 February 1996 358 A. CARDONA Xiao and Yang6 used the concept of velocity screws and dual matrices to derive expressions to describe the kinematics of a pair of gears with non-parallel and non-intersecting axes. In this paper we develop the equations of motion for representing three-dimensional gearing, using the formulation of a general purpose code for exible multibody analysis.7; 8 The gear is implemented as a particular element of the code, adding a new functionality to it. The code takes into account full three-dimensional motion, elasticity of members, damping, friction dissipation and many other eects. Code elements comprise a beam, membrane, rigid body, spring, hinge, cam, screw, prismatic, and cylindrical joints within many others, allowing to make a highly realistic model of any machine. It can be linked to nite element packages for dynamic analysis, incorporating thus easily arbitrary nite element models of elastic subsystems—either xed or movable—to the numerical model of the mechanism. The joint is described by two physical nodes: one at the centre of each wheel. Each node can be connected in turn to the other members of the mechanism. The following simplifying hypotheses have been made in order to build the numerical model of the gear pair: 1. Bodies in contact are considered rigid (gear teeth exibility have been neglected). 2. Time-varying eects as uctuation of the number of tooth pairs in contact are neglected. 3. Gear backlash is neglected (zero teeth separation is assumed). Also, gear teeth spacing errors and misalignment are not considered by the model. Although these simplifying hypotheses can be eliminated, in this paper we have centred the study only on geometric aspects of three-dimensional gearing for conciseness. The formulation developed for the gear pair is general enough to represent almost any kind of gears used in industry: spur gears, helical gears, conical gears, internal gear and pinion, worm gears, hypoid gears, etc. All components of contact forces induced by the gear pair are taken into account: tangential, radial and thrust forces. These forces are appropriately transmitted to the element nodes and to the other members of the mechanism. Friction forces are taken into account whenever appropriate (e.g. worm gearing). They are a function of the transmitted force, and of a friction coecient. The friction coecient varies with the instantaneous value of sliding velocity between worm and gear. Gear trains and planetary systems can be easily represented by superposing several gear pairs in the mechanism model. The paper is organized in two parts. The rst part treats the case of gearing engagement between two wheels. The second part analyses the case in which one wheel has innite radius, i.e. rack /gear engagement. Several examples of application are shown, illustrating the generality and power of the approach. 2. FORMULATION OF THE CONSTRAINED DYNAMICS PROBLEM The equations of motion for a dynamic system subjected to holonomic constraints can be stated in the form8; 9 @L d @L − = Q + Q0 dt @ q̇ @q (1) k ( q̇; q; t) = 0 where an augmented Lagrangian approach is followed, L is the Lagrangian of the unconstrained system (it accounts, for instance, for the internal strain and kinetic energies of the elastic members and for the potential energy of the conservative loads); q and q̇ are the generalized displacements and velocities of the dynamic system, and Q are the non-conservative forces acting on the structure THREE-DIMENSIONAL GEARS MODELLING IN MULTIBODY SYSTEMS ANALYSIS 359 (they include, for instance, the friction and damping forces). The constraint forces Q0 , which can be thought of as the generalized forces that oblige the system to verify the imposed holonomic constraints, are Q0 = P i (ki − pi ) @i @q (2) where k and p are the scale and the penalty factors, respectively, and i is the Lagrange multiplier associated to the constraint i . The penalty term does not modify the stationary properties of the original Lagrangian function (i.e. the value of the solution is not altered) but adds some positive curvature in the range space of @i =@q, with a signicant improvement of convergence of the iterative scheme.8; 9 The values of k and p are selected to balance adequately the equations: usually, a mean value of the stiness of the mechanism is an appropriate choice for them. Equations (1) form a system of non-linear dierential-algebraic equations. The solution is advanced step-by-step using a time integration scheme (in particular, we used the Hilber–Hughes– Taylor method, see References 10 and 11 for details). At each time step, a system of non-linear algebraic equations has to be solved. In order to solve this system, we use the Newton method. The contribution of the holonomic constraints to the Jacobian matrix is obtained by dierentiating constraint forces and constraints with respect to the generalized displacements vector # "P # " 0 pBBT −kB i (ki − pi ) ∇qq i + (3) S= 0 −kBT 0 0 where B = [ @1 =@q @2 =@q · · · ]. The rst term on the right-hand side is formed only by rstorder derivatives of the constraints, while the second one comprises the second-order derivatives. Experience has shown that neglecting the second-order derivatives of constraints, does not alter signicantly the convergence rate of standard Newton iteration; therefore, these terms are usually ignored.8; 9 PART I: GEARED WHEELS PAIR 3. KINEMATICS OF THE JOINT 3.1. Basic denitions We will consider a joint formed by two geared wheels with centres A and B. We dene in this section some elements that will be used to express the kinematic relations at the joint. The position of each wheel centre in the inertial frame is given by xA and xB . Let us dene one triad of unit orthogonal vectors at each wheel centre at the reference conguration: (a) {\1 ; \ 2 ; \ 3 } is attached to the rst wheel, with origin at its centre A; (b) {^1 ; ^ 2 ; ^ 3 } is attached to the second wheel, with origin at node B. Both triads are dextrorsum and have their rst vector perpendicularly oriented to the wheel plane and their second vector oriented towards the contact point between wheels (see Figure 1). Explicit expressions for vectors \ 2 ; \ 3 ; ^ 2 ; ^ 3 obtained according these guidelines, are given afterwards in Section 3.2. Remark. The distance between centres at the current conguration as well as the relative orientation of both wheels, should be kept constant by some external means (e.g. a frame) to maintain the correct engagement. 360 A. CARDONA Figure 1. Gear pair kinematics The orientation of both material triads at the current conguration is obtained applying the rotation operator of each wheel to the vectors at the reference conguration: \0i = RA \ i ; ^0i = RB ^ i ; i = 1; 3 (4) RA ; RB are the rotation matrices at nodes A; B, which are related to the rotation vectors A ; B e A ); exp( e B ).8; 12 through the exponential forms exp( We dene a second pair of unit triads at the current conguration. These triads are attached to the supporting frame of the gear pair. They are built following the same rule we used to get the reference triads {\1 ; \ 2 ; \ 3 } and {^1 ; ^ 2 ; ^ 3 }, i.e., vectors \100 and ^100 are normal to each wheel: \100 = \10 = RA \1 ; ^100 = ^10 = RB ^1 (5) while vectors \002 ; ^002 point towards the contact point and vectors \003 ; ^003 complete a dextrorsum triad. We also dene a triad of unit vectors {W100 ; W002 ; W300 }, which are oriented along the teeth in contact at the current conguration: W100 is parallel to the tooth baseline; W002 points along the tooth vertical line, from the rst wheel to the second one; W300 is normal to the tooth midplane. The triad {W00i } is related to {\00i } through the angles of cone A and helix A , as shown in Figure 1: [ W100A W200A W300A ] = [ \100 \002 \003 ] YA (6) THREE-DIMENSIONAL GEARS MODELLING IN MULTIBODY SYSTEMS ANALYSIS 361 Figure 2. Sign conventions for cone and helix angles The explicit expression of the orthogonal matrix YA follows: cos A cos A YA = Y(A ; A ) = sin A cos A − sin A cos A sin A sin A sin A cos A − sin A 0 (7) cos A The sign conventions for the cone and helix angles are described in Figure 2. are here Remark. (1) The superscript A in equation (6) is used to emphasize that vectors W00A i computed in terms of kinematic variables at wheel A only. (2) Note that internal gears can be represented by selecting A = 180◦ . Similar relations hold between {W00i } and {^00i } through the angles B and B , with the only dierence that signs are changed to match the condition of gear engagement (note that the triad W00i has been dened as pointing from the rst wheel to the second one): [ W100B W00B 2 00 W00B 3 ] = [ ^1 ^002 1 ^003 ] Y(B ; B ) 0 | 0 {z YB 0 −1 0 0 0 −1 } (8) 362 A. CARDONA We get the relation between triads {\00i } and {^00i } at engagement by equating equations (6) and (8): [ ^100 ^002 ^003 ] = [ \100 \002 \003 ] YA YBT = [ \100 \002 \003 ] Z(A ; A ; B ; B ) (9) Remark. It can be easily shown that Z is in fact a function of the three independent parameters (A ; B ; A − B ). 3.2. Explicit expressions of vectors \ 2 ; \ 3 ; ^ 2 ; ^ 3 at the reference conguration First, we determine the position of the contact point xC . It can be computed in terms of kinematic variables at wheel A or at wheel B: xCA = xA + rA = xA + \002 rA xCB = xB + rB = xB + ^002 rB (10) where rA ; rB are the radius vectors from the wheels centres to the contact point rA = rA \002 ; rB = rB ^002 Using equation (9), we compute vector ^002 in terms of vectors \00i : P 00 \ i Z i;2 ^002 = (11) (12) i=1;3 Then, we can write the equation for the relative position of xB with respect to xA as indicated below: P 00 \ i Z i;2 rB (13) xAB = xB − xA = \002 rA − i=1;3 \003 Noting that of equations = \100 × \002 , and since \ i = \00i at the reference conguration, we obtain the system e1 ] \ 2 = −xAB − rB Z1;2 \1 [(rB Z2;2 − rA )1 + Z3;2 rB \ (14) We get the expression of \ 2 by solving the latter system: 2 −1 e1 ] {xAB + rB Z1;2 \1 } ) [(Z 2;2 rB − rA )1 − rB Z3;2 \ \ 2 = −((Z2;2 rB − rA )2 − rB2 Z3;2 = c1 xAB + c2 \1 + c3 \1 × xAB (15) where constants (c1 ; c2 ; c3 ) are dened as follows: c1 = (rA − Z2; 2 rB )=c c2 = (rA − Z2; 2 rB )Z1; 2 rB =c c = (rA − Z2; 2 rB )2 − Z3;2 2 rB2 (16) c3 = Z3; 2 rB =c The third vector \ 3 is then computed: \ 3 = \1 × \ 2 = c1 \1 × xAB − c3 [1 − \1 ⊗ \1 ] xAB (17) THREE-DIMENSIONAL GEARS MODELLING IN MULTIBODY SYSTEMS ANALYSIS 363 Vectors ^ 2 ; ^ 3 are afterwards computed using equation (9): [^1 ^ 2 ^ 3 ] = [\1 \ 2 \ 3 ]Z. Remark. It can be shown numerically that c 6= 0 for rA ¿rB . Therefore, we ask this condition be satised to assure the system of equations (14) is solvable. 3.3. Degrees of freedom of the joint The pair is formed by two rigid bodies (the two wheels) related by one kinematic constraint that links the rotations of both wheels about their respective normal axes. Fourteen kinematic variables are used to describe the joint, grouped into the generalized co-ordinates vector q: q = hxTA TA xTB TB A Bi (18) Recall xA ; xB give the position and A ; B give the rotation vector of each wheel centre. The angular displacements A ; B measure the relative rotation of each wheel in the local frame, resulting in the following relations between quoted and doubly quoted frames: \002 = \02 cos ^002 = ^02 cos − \03 sin A B − ^03 sin A; \003 = \02 sin B; ^003 = ^02 sin A B + \03 cos + ^03 cos A (19) B These relations can be written in matrix form [ \100 \200 \300 ] = [ \10 \02 \03 ] R(e1 ; A) [ ^100 ^200 ^300 ] = [ ^10 ^02 ^03 ] R(e1 ; B) with 1 0 0 R(e1 ; ) = 0 cos − sin 0 sin cos (20) (21) Since the dimension of the vector of kinematic variables q exceeds by three the number of physical degrees of freedom of the joint, we will have to impose a total of three constraints. One constraint describes the kinematic relation resulting from teeth contact, and two additional ones x the value of the relative angular displacements A ; B . This set of constraints is described in the next section. Three Lagrange multipliers—conjugated to these constraints—are added to the generalized co-ordinates q to form the vector of unknowns of the joint. 3.4. Variation of vectors \00i ; ^00i ; W00i By taking variations in equation (10), we can write xCA = xA + \002 rA (22) xCB = xB + ^002 rB The variation of vector \002 is computed using equation (19): \002 = \02 cos A − \03 sin = (RA A ) × \002 − A − \02 sin \003 A A A − \03 cos A A (23) 364 A. CARDONA For the rest of the doubly quoted vectors we get, in the same fashion, \003 = (RA A ) × \003 + \002 ^002 = (RB B ) × ^002 − A ^003 B ^003 = (RB B ) × ^003 + ^002 (24) B After replacing the variations of \002 ; ^002 in equation (22), we can compute the variations of position of the contact point computed in terms of kinematic variables at wheel A or at wheel B: xCA = xA + (RA A ) × \002 rA − \003 A rA xCB = xB + (RB B ) × ^002 rB − ^003 (25) B rB Equations (23) and (24) can be written in matrix form [ \100 \200 00 \300 ] = [R] A A ] [ \1 \002 \003 ] + A [ 0 − \003 \002 ] [ ^100 ^200 00 ^300 ] = [ R] B B ] [ ^1 ^002 ^003 ] + B[ 0 − ^003 ^002 ] (26) 00B Finally, using equations (6)–(8) we establish the expressions for the variations of vectors W00A i ; Wi : [ W100A W200A 00A W300A ] = [R] A A ][ W1 W00A 2 W00A 3 ] + A[ 0 −\003 \002 ] YA [ W100B W200B 00B W300B ] = [ R] B B ][ W1 W00B 2 W00B 3 ]+ −^003 ^002 ] YB B[ 0 (27) 4. COMPUTATION OF THE CONSTRAINT FORCES Three constraints are needed to determine fully the joint kinematic variables. The rst equation of constraint gives the kinematic relation between angular displacements of both wheels: 1 = (− A zA + B zB ) mn cos n = 0 2 (28) where mn is the normal module of the gear teeth, n the pressure angle in the normal plane, and zA and zB are the numbers of teeth at each wheel. It can be shown easily that by formulating the constraint in this way, the conjugated Lagrange multiplier—scaled by factor k—has the physical meaning of normal contact force: F = k1 . The second constraint represents the hoop contact constraint produced by engagement between teeth: 2 = (xCA − xCB ) · W00A 3 =0 (29) where xAC ; xBC give the position of the contact point computed on wheel A and B, respectively. The third constraint expresses that the triad W00i is unique, when computed in terms of kinematic variables of either wheel: 00B 3 = W00A 2 · W3 = 0 (30) THREE-DIMENSIONAL GEARS MODELLING IN MULTIBODY SYSTEMS ANALYSIS 365 The variation of 1 is simply written 0 xA 0 A 0 xB @1 1 = q · = · @q 0 B 1 − zA mn cos n 2 A 1 z m cos B B n n 2 (31) The variation of the second constraint can be written in the form 00A A B 2 = (xCA − xCB ) · W00A 3 + W 3 · (xC − xC ) (32) Using equations (25) and (27), we get xA W00A 3 T 00A B R (W × (x − x )) A A A 3 C 00A xB −W 3 @2 = · 2 = q · @q RTB (W00A B 3 × rB ) 00 00 B A YA 3;3 \ 2 − YA 2;3 \ 3 · (xA − xC ) 00 00A B rB ^ 3 · W 3 (33) Variations of the third constraint are next computed 00B 00B 00A 3 = W00A 2 · W 3 + W 3 · W 2 giving 3 = q · xA A x 0 00B RTA (W00A 2 × W3 ) (34) 0 @3 B = · 00B @q −RTB (W00A B 2 × W3 ) 00 00B \ · W −Y A A 2;2 3 3 00 00 00A B YB 3;3 ^ 2 − YB 2;3 ^ 3 · W 2 (35) Afterwards, by substituting equations (31), (33) and (35) into equation (2), we evaluate the vector of constraint forces Q0 of the joint: (k1 − p1 ) @1 @2 @3 0 Q = (36) (k2 − p2 ) = B{ k[ − p } @q @q @q (k3 − p3 ) The expression for the stiness contribution follows immediately from equation (3). 366 A. CARDONA 5. RADIAL COMPONENT OF THE CONTACT FORCE Contact between teeth is produced along the pressure line. The orientation of the pressure line varies with the sign of the transmitted torque, in such a way that the contact force can be seen as always trying to separate both wheels. The contact force lies in the plane {W200 ; W300 }, normally oriented to the teeth in contact. Its magnitude F is equal to k1 (the Lagrange multiplier conjugated to the rst constraint 1 times the scale factor k). The hoop and axial components—oriented along W003 —have been taken into account when formulating the holonomic constraint 2 . However, the radial component of the contact force is of non-holonomic nature and has to be added explicitly to the formulation as a non-conservative force. The radial component of force acting on wheel A at the contact point is F = −|F| sin n W002 = −k|1 | sin n W002 (37) while the opposite force −F acts upon wheel B at the same point. We remark that F can be expressed indistinctly in terms of variables at either wheel. We will note. F A = −k|1 | sin n W00A 2 (38) F B = −k|1 | sin n W00B 2 and use one or the other expression when appropriate. The radial contact force is conjugated to the variations of position of the contact point C, resulting in the following virtual work expression : W = xAC · FA − xBC · FB (39) = xA · FA + A · RTA (rA × FA ) − xB · FB − B · RTB (rB × FB ) where we have used the identities: \003 · FA = 0; ^003 · FB = 0. The vector of non-conservative forces Q0 of the gear pair is identied as formed by the conjugated components to the variations of the generalized co-ordinates of the joint xA FA T A R (r × F ) A A A B x −F B 0 (40) W = q · Q = · T B B −RB (rB × F ) 0 A B 0 The contribution to the tangent stiness is obtained by dierentiation of equation (40). We compute rst the increment of the contact force FA : FA = −k|1 | sin n W200A − sign(1 )k sin n W200A 1 T 00 00 = k|1 | sin n W̃00A 2 RA A + k|1 | sin n −YA 2;2 \ 3 + YA 3;2 \ 2 A −sign(1 )k sin n W200A 1 A = −F̃ RTA A − cA \003 A + (F A =1 )1 (41) 367 THREE-DIMENSIONAL GEARS MODELLING IN MULTIBODY SYSTEMS ANALYSIS with cA = k|1 | sin n YA 2;2 (42) Similarly, B F B = −F̃ RTB B − cB ^003 B + (FB =1 )1 (43) with cB = k|1 | sin n YB 2;2 (44) The increments of moments are computed next. First note that (RTA rA ) = −rA RTA \003 (RTA FA ) = A RTA (FA =1 )1 − cA RTA \003 (45) A Therefore, A (RTA (rA × FA )) = RTA r̃A (F A =1 )1 + RTA (rA F̃ − cA r̃A )\003 A (46) B (47) Also, B (RTB (rB × FB )) = RTB r̃B (FB =1 )1 + RTB (rB F̃ − cB r̃B )^003 Equations (41), (43), (46) and (47) give the contributions to the rows of the (non-symmetric) tangent stiness matrix S, such that ( ) " #( ) S q S q qq q · (48) 2W = Sq S where Sqq = 0 −F̃ RTA A 0 0 −cA \003 0 0 0 0 RTA (rA F̃ − cA r̃A )\003 0 0 0 −F̃ RTB B 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 A Sq 0 00 −cB ^ 3 B T 00 RB (rB F̃ − cB r̃B )^ 3 0 0 0 (FA =1 ) 0 0 0 0 0 RT r̃ (FA = ) 0 0 1 A A (FB = ) 0 0 1 = RT r̃ (FB = ) 0 0 1 B B 0 0 0 Sq = S = 0 (49) 368 A. CARDONA 6. FRICTION BETWEEN WHEELS Several gear designs have a non-zero mutual sliding speed along teeth in contact (e.g. worm gearing). In such cases, there appears a non-negligible friction force between wheels. Let us suppose the system is working with a contact force of modulus |F| = k|1 |; the friction force acting on wheel A at the contact point can be written Ffr = −sign(vrel )(vrel )k|1 |W100 (50) where vrel is the relative sliding speed, measured along the teeth longitudinal direction W100 : vrel = (ẋA + (RA A ) × rA − ẋB − (RB B ) × rB ) · W100 The contribution to the internal forces due to friction then reads xA 1 RT r̃ A A A x −1 B q · Qfr = · Ffr T B −RB r̃B A 0 B 0 (51) (52) We note that the internal friction forces depend both on the generalized displacements and on the generalized velocities. When computing the tangent matrices, we will neglect the contribution to the stiness (derivatives with respect to the generalized displacements) and compute only the derivative of the internal forces with respect to the generalized velocities, which result in the tangent damping matrix C. First note that Ffr = −R (vrel ) W100 k sign(1 ) 1 − R (vrel ) k |1 | W100 − R0 (vrel ) k |1 | W100 vrel (53) The increment of sliding speed reads vrel = W100 · (ẋA − rA × RA A − ẋB + rB × RB B ) ẋA A T T = h W100T W100T r̃A RA −W100T −W100T r̃B RB i ẋB B (54) By replacing the latter into equation (53), we get @Ffr q̇ = −R0 (vrel ) k |1 | W100 h W100T @q̇ W100T r̃TA RA −W100T ẋA A T −W100T r̃B RB i ẋB B (55) THREE-DIMENSIONAL GEARS MODELLING IN MULTIBODY SYSTEMS ANALYSIS 369 Then, the tangent damping matrix can be written in the form W100 T 00 r̃ R W A 1 A −W00 1 C = −R0 (vrel ) k |1 | −RTB r̃B W100 0 0 T W100 T 00 r̃ R W A 1 A −W00 1 −RTB r̃B W100 0 0 (56) PART II: GEAR-RACK PAIR 7. KINEMATICS OF THE JOINT 7.1. Basic denitions Let us describe briey the modications to be introduced in order to model a gear-rack system. The joint connects two nodes: node A, a material node on the rack, and node B, located at the wheel centre. The position of each node of the kinematic pair in the inertial frame is given by xA and xB . We dene a triad of unit orthogonal vectors at each node at the reference conguration: {\1 ; \ 2 ; \ 3 } is attached to the rack, with origin at node A; {^1 ; ^ 2 ; ^ 3 } is attached to the centre of the wheel. Vector \ 3 is oriented along the rack, and vector \ 2 points along the tooth vertical line from the rack to the wheel. Triad {^ i } is oriented following the same convention as for the wheel pair (see Figure 3). The relative orientation between rack and wheel should be kept constant by an external means to maintain the correct engagement of both members (e.g. a frame). The orientation of both material triads at the current conguration is obtained using the rotation of each member: \0i = RA \ i ; ^0i = RB ^ i ; i = 1; 3 (57) We dene a second triad {^100 ; ^002 ; ^003 } at the current conguration for wheel B, following the same conventions as in the two-wheel case. We also dene a triad of vectors {W100 ; W200 ; W300 }, oriented along the teeth in contact at the current conguration (W100 parallel to the tooth baseline, W200 points along the tooth vertical line, from the rack to the wheel, and W300 is normal to the tooth midplane). The triad {W00i } is related to {\0i } through the helix angle A : [W100A W200A W300A ] = [\10 \02 \03 ] YA (58) with YA = Y(0; A ) and where the explicit expression of the orthogonal matrix Y(A ; A ) was given in equation (7). The cone angle for the rack is 0, since vector \02 is assumed coincident with W002 . The superscript A is used to emphasize that these vectors are computed in terms of kinematic variables at the rack A only. 370 A. CARDONA Figure 3. Kinematics of the gear-rack pair Similar relations hold between {W00i } and {^00i } through the angles B and B , with the only dierence that, again, signs are changed to match the condition of gear engagement: 1 0 0 (59) [W100B W200B W300B ] = [^100 ^002 ^003 ] Y(B ; B ) 0 0 −1 0 {z YB | 0 −1 } By using equations (58) and (59), we get the following relation between triads {\0i } and {^00i } at engagement: [ ^100 ^002 ^003 ] = [ \10 \02 \03 ] YA YTB = [ \10 \02 \03 ] Z(A ; B ; B ) (60) 7.2. Degrees of freedom of the joint The vector of generalized co-ordinates of the gear-rack pair is given next: q = h xTA TA xTB TB uA B i (61) where uA measures the displacement of the rack with respect to the wheel, and rotation of the wheel B with respect to the pair frame. B gives the 7.3. Variation of vectors \0i ; ^00i ; W00i First variations of vectors \0i and ^00i can be written in the following form: [ \10 \20 0 \30 ] = [ R] A A ][ \1 \02 \03 ] [ ^100 ^200 00 ^300 ] = [R] B B ][ ^1 ^002 ^003 ] + B[ 0 −^003 ^002 ] (62) 371 THREE-DIMENSIONAL GEARS MODELLING IN MULTIBODY SYSTEMS ANALYSIS Using equations (58) and (59) it is directly seen that [ W100A W200A 00A W300A ] = [R] A A ][ W1 W00A 2 W00A 3 ] [ W100B W200B 00B W300B ] = [R] B B ][ W1 W00B 2 W00B 3 ]+ B [0 (63) − ^300 ^200 ] YB The position of the contact point xC can be computed in terms of kinematic variables at either rack A or wheel B: xAC = xA + \02 h + \03 uA = xA + rA (64) xBC = xB + ^002 rB = xB + rB where h is the (constant) rack thickness and uA is the displacement of the rack with respect to the wheel along axis \03 , and where rA ; rB are the vectors of relative position of the contact point with respect to each member node: rA = \02 h + \03 uA ; rB = rB ^002 (65) By taking variations in equation (64), we get xAC = xA + (RA A ) × rA + \03 uA xBC = xB + (RB B ) × rB − ^003 B (66) rB 8. COMPUTATION OF THE CONSTRAINT FORCES The rst constraint gives the kinematic relation between internal variables uA and B zB mn cos n = 0 1 = uA cos A + 2 B: (67) The second constraint represents contact between teeth: 2 = (xAC − xBC ) · W00A 3 =0 (68) where xAC ; xBC give the position of the contact point computed on the rack A and on the wheel B, respectively. The third constraint expresses that the triad W00i is unique, when computed in terms of kinematic variables from either side: 00B 3 = W00A 2 · W3 = 0 Constraints variations 1 ; 2 ; 3 are detailed next: xA 0 0 A 0 xB @1 1 = q · = · @q 0 B cos u cos A A n 1 B z m cos B n n 2 (69) 372 A. CARDONA xA W00A 3 T 00A B R (W × (x − x )) A A A 3 C 00A xB −W 3 @2 = 2 = q · · @q B RTB (W00A 3 × rB ) u Y 00 B A A 3;3 \ 2 · (xA − xC ) 00 A B rB ^ 3 · W 3 (70) xA 0 T 00A 00B RA (W 2 × W 3 ) A xB 0 @3 = · 3 = q · 00B @q −RTB (W00A B 2 × W3 ) u 0 A 00 00 00A B YB 3;3 ^ 2 − YB 2;3 ^ 3 · W 2 9. NORMAL COMPONENT OF THE CONTACT FORCE Engagement between rack and wheel produces a contact force, whose magnitude F = k1 is given by the Lagrange multiplier conjugated to the rst constraint times the scale factor k and its orientation follows that of the pressure line. The orientation of the pressure line depends on the sign of the transmitted torque, such that the contact force always tries to separate both members. The normal component of contact force is, therefore, of non-holonomic nature. The force acting on the rack A at the contact point results F = −k|1 | sin n W200 (71) while the opposite force −F acts upon wheel B at the same point. F can be expressed in terms of variables of either side, and we note FA = −k|1 | sin n W200A FB = −k|1 | sin n W200B (72) We will use one or the other expression when appropriate. The virtual work expression of the non-holonomic constraint forces follows: W = xAC · FA − xBC · FB = xA · FA + A · RTA (rA × FA ) − xB · FB − B · RTB (rB × FB ) (73) THREE-DIMENSIONAL GEARS MODELLING IN MULTIBODY SYSTEMS ANALYSIS 373 where we have used the identities: \03 · FA = 0; ^003 · FB = 0. In vector form FA xA T A A RA (rA × F ) B xB −F · W = q · Q = −RTB (rB × FB ) B u 0 A B 0 (74) where Q is the contribution to the vector of internal forces of the joint due to the normal component of the contact force. The contribution to tangent stiness is obtained by dierentiation of equation (74). Derivatives of FB and of RTB (rB × FB ) are given in equations (43) and (47). The increment of the contact force FA reads FA = −k|1 | sin n W200A − sign(1 )k sin n W200A 1 T 00A = k|1 | sin n W̃00A 2 RA A − sign(1 )k sin n W2 1 = −F̃ A RTA A (75) A + (F =1 )1 The increment of moment at node A can be written in the following form: (RTA (rA × F A )) = RTA r̃ A (FA =1 )1 (76) Therefore, the contribution to the (non-symmetric) tangent stiness matrix S of terms associated to the normal component of the contact force in the rack-gear pair can be written in the form ( W= 2 q ) " · Sqq Sq Sq S #( q ) (77) where Sqq 0 0 0 = 0 0 0 −F̃ RTA A 0 0 0 0 0 0 0 −F̃ RTB 0 0 0 0 0 0 0 0 0 0 0 0 0 −cB ^003 B 0 RTB (rB F̃ − cB r̃B )^003 0 0 0 0 0 B 374 A. CARDONA Sq (FA =1 ) 0 0 RT r̃ (FA = ) 0 0 1 A A (FB = ) 0 0 1 = RT r̃ (FB = ) 0 0 1 B B 0 0 0 0 0 (78) 0 Sq = S = 0 Expressions for friction terms between rack and wheel are entirely similar to those developed for the geared wheels pair. 10. EXAMPLES The contributions of the gear pair element are assembled to the rest of the mechanical system by following a standard procedure, leading to a system of dierential /algebraic equations (DAEs). This system is in turn time integrated using the Hilber–Hughes–Taylor algorithm (the whole procedure is described in References 10 and 11). We present here four examples of application of the techniques for modelling gear pair systems. In the three rst cases, we have neglected dynamic eects (i.e. we have assumed zero mass) making kinematic analyses in which the displacement at some degrees of freedom of the system are imposed to determine a new conguration. The resulting system of equations is solved using the Newton method. The fourth case corresponds to a dynamics simulation: the analysis of a spur gears transmission including exibility of axes and bearing. 10.1. Conical straight bevel pair The rst example consists of a conical straight bevel gear pair, in which one wheel is xed and the other one is mounted on a rigid bar which turns around the y-axis (Figures 4 and 5). The centre of the rst wheel is xed at the co-ordinates origin. The wheel is lying in the x–z plane. It has radius rA = 5, normal modulus mn = 0·2, pressure angle n = 20◦ , teeth number zA = 50, cone angle A = − tan−1 (5) = −78·7◦ and helix angle A = 0. At the initial conguration, the second wheel has its centre at (0; 1; 5) and lies in a plane parallel to the x–y plane. Its radius is rB = 1, the teeth number zB = 10, cone angle B = − tan−1 (0·2) = −11·3◦ and helix angle B = 0. It is hinged to a rigid bar which goes from (0; 1; 5) to the point of co-ordinates (0; 1; 0). At this point, the bar is connected through a hinge with axis along the y-axis to the foundation. The system has one degree of freedom. A new conguration is found by imposing an incremental angle of 1 rad to the rigid bar. Also, a torque Mext = 1000 is imposed at the articulation between the bar and the second wheel. Figure 4 displays the computed conguration together with the resulting actions of the rst wheel upon the foundation. The computed normal contact force between gears was Fgear = 1·064 × 103 . 10.2. Three-wheeled system The second example consists of three equal spur gears mounted on a rigid bar (Figures 6 and 7). The system is lying in the x–y plane. The bar is articulated to the foundation at the co-ordinates THREE-DIMENSIONAL GEARS MODELLING IN MULTIBODY SYSTEMS ANALYSIS 375 Figure 4. Conical bevel pair Figure 5. Schematic view of conical bevel pair (program output) origin. A string (kstring = 1000; l0 = 50) connects the extreme of the bar to the foundation at point (30; −50; 0) (not shown in Figure 6). The three wheels have normal modulus mn = 0·2, pitch diameter d = 5, pressure angle n = 20◦ and teeth number z = 50 (clearly, the cone and helix angles are = = 0). The rst wheel— centred at the co-ordinates origin—is xed to the foundation, while the rigid bar is articulated at the same point. The bar is aligned along the x-axis at the initial conguration. The system is put in motion by imposing a relative angular displacement rel = 1 rad at the third wheel, thus inducing a rotation of 1 rad of the bar about the origin as shown in Figure 6. At the same gure we display the values of actions of the rst wheel on the foundation for the new conguration, and also the value of force exerted by the string. The computed value of contact force between gears was Fgear = 1·156 × 105 . 376 A. CARDONA Figure 6. Three-wheeled system Figure 7. Schematic view of three-wheeled system (program output) 10.3. Double rack system This third example consists of a conical straight bevel wheel xed to the foundation at point (−3; 0; 3), which connects two racks oriented parallel to the y-axis. One rack passes through point (0; 0; 6) with teeth normal to the y–z plane and oriented towards the negative x semiaxis; the second rack passes through point (−6; 0; 0) with teeth normal to the x–y plane oriented towards the positive z semiaxis (Figure 8). √ ◦ The cone angle of the wheel is B = −45◦ , its radius is rB = 6= 2, the pressure angle √ n = 20 , the teeth number is zB = 30, helix angle B = 0, and normal modulus mn = 0·4= 2. The helix angles of both racks is 0. THREE-DIMENSIONAL GEARS MODELLING IN MULTIBODY SYSTEMS ANALYSIS 377 Figure 8. Double rack system Figure 9. Spur gears transmission Both racks are restricted to move parallel to the y-axis. Therefore, the system has only one kinematic degree of freedom. A new conguration is found by imposing an incremental displacement uimp = 5 to one rack. Also, an axial force Fext = 1000 is imposed to the second rack. Figure 8 displays the computed conguration together with the resulting actions of the connecting wheel upon the foundation. 10.4. Dynamic analysis of a spur gear transmission The last example concerns the analysis of a transmission composed by two axes linked through a pair of equal spur gears (see Figure 9). Most data in this example correspond to an experimental uven.14 However, the analysis set-up made by Kubo et al.,13 and was used also as test-case by Ozg we are presenting here does not account for the exibility of teeth; therefore, we are not able to 378 A. CARDONA Figure 10. Applied torque at wheel 2 compare results with those of the cited references. An extension of our formulation to include teeth exibility, clearance and loaded error transmission is described in Reference 15. The system is driven at a constant speed of 12 000 rpm at wheel 1, and transmits a torque T2 applied at wheel 2, whose time variation law is dened as follows (see Figure 10): 945·8; t ¡ 0·001 T2 = (79) 945·8 + 100 sin(104 t); t ¿ 0·001 Gear properties are: normal modulus mn = 0 · 15748, pitch diameter d = 3 · 937, pressure angle n = 20◦ and teeth number zA = zB = 25 (cone and helix angles are trivially = = 0). Their mass and rotary inertia are: mA = mB = 5·36 × 10−3 and IA = IB = 0·0102. Torsional stiness of shafts are: kt1 = 1·7 × 104 , kt2 = 3·0 × 104 , and damping ct1 = 0·07603, ct2 = 0·1237. Disks at the extremes of axes have rotary inertias I1 = 0·051 and I2 = 0·0102. The axes are mounted on exible supports, with stiness values k1 = k2 = 1·477 × 107 and damping c1 = c2 = 50. The system response is almost stationary. To compute these results, the system equations were time integrated up to a time instant in which all initial transients have been damped out. In order to reach quickly the stationary phase, we have performed the following computation steps: rst, we evaluated the static deformation under the initial constant value of torque T2 (= 945·8); secondly, we made a kinematic analysis in which the system was put into motion at 12 000 rpm, by ignoring all inertial forces and starting from the deformed conguration computed at the rst computation step; nally, a full transient analysis was performed, starting from the conguration and velocities computed at the previous phase of analysis. Figure 11 displays the time evolution of the torque transmitted by axes 1 and 2. We can see that the frequency content of torque at axis 2 is much higher than that of axis 1. The oscillations evidenced in axis 2 are produced by coupling between the response of the system supports and the perturbation in the excitation of frequency × 104 Hz. We can see that the frequency of response at axis 1 is much lower than that of axis 2, and is determined by the axis lowest torsion vibration eigenfrequency. THREE-DIMENSIONAL GEARS MODELLING IN MULTIBODY SYSTEMS ANALYSIS 379 Figure 11. Transmitted torque at axes 1(a) and 2(b) The reaction forces at the axes supports are displayed in Figure 12. Finally, the time evolution of the total mesh force at the spur gears pair is shown in Figure 13. We can appreciate some high-frequency oscillations at t = 0 and at t = 0·001 which are purely numerical and are rapidly damped out by the algorithm. 11. CONCLUDING REMARKS A general methodology to formulate gear pair connections in three-dimensional kinematic and dynamic analysis of mechanisms has been proposed. The joint was formulated using a set of holonomic and non-holonomic constraints. Holonomic constraints were imposed by an augmented Lagrangian technique, while non-holonomic radial and friction forces were explicitly added to the Euler–Lagrange equations of the dynamic system. 380 A. CARDONA Figure 12. Horizontal (a) and vertical (b) components of the reaction forces at the axes supports The joint discussed here forms part of the vast library of joints of the multiple purpose program for mechanism analysis MECANO,7 and adds new functionalities to the code. The capability of modelling three-dimensional gears and to connect them with the rest of (exible and rigid) joints and bodies of the code results in a very valuable tool for the design of mechanisms and machines. Four examples of application have been shown. The kinematics of mechanisms including conical straight bevel pairs, three spur gear wheels mounted on a moving rigid bar, and a system composed by two racks connected by a conical straight bevel gear has been analysed. Finally, the dynamic analysis of a spur gears transmission has been made, illustrating the generality of the formulation. The gear pair joint developed in this work is rigid (i.e. it consists of an ideal connection between rigid gears). Some further developments to include eects as teeth exibility, stiness uctuation THREE-DIMENSIONAL GEARS MODELLING IN MULTIBODY SYSTEMS ANALYSIS 381 Figure 13. Time evolution of the total mesh force and gear backlash were reported in Reference 15. Future research will be directed to improve the capability of modelling teeth exibility and dierent eects of interest in gear dynamics. ACKNOWLEDGEMENTS This work received nancial support from Conicet through grant PID-BID 238 and from Fundacion Antorchas through grant A-13015/1-16. REFERENCES 1. J. E. Shigley and J. J. Uicker, Jr., Theory of Machines and Mechanisms, McGraw-Hill, New York, 1980. 2. C. Rama Mohana Rao and G. Muthuveerappan, ‘Finite element modeling and stress analysis of helical gear teeth’, J. Comput. Struct., 49, 1095–1106 (1993). 3. A. Kahraman, ‘Load sharing characteristics of planetary transmissions’, Mech. Mach. Theory, 29, 1151–1165, (1994). 4. A. Kahraman, ‘Planetary gear train dynamics’, J. Mech. Design, 116, 713–720 (1994). 5. M. Geradin and G. Robert, ‘Vibrations de Torsion et de Flexion d’un Train d’Engrenages’, Rapport VF-42, LTAS, Universite de Liege, Belgium, 1981. 6. D. Z. Xiao and A. T. Yang, ‘Kinematics of three dimensional gearing’, Mech. Mach. Theory, 24, 245–255 (1989). 7. SAMCEF—Module d’Analyse de Mecanismes MECANO (Manuel d’Utilisation), 1994. 8. A. Cardona, ‘An integrated approach to mechanism analysis’, Ph.D. Thesis, Universite de Liege, Belgium, 1989. Collection des Publications de la Faculte des Sciences Appliquees 127. 9. A. Cardona, M. Geradin and D. B. Doan, ‘Rigid and exible joint modeling in multibody dynamics using nite elements’, Comput. Methods Appl. Mech. Eng., 89, 395–418 (1991). 10. A. Cardona and M. Geradin, ‘Time integration of the equations of motion in mechanism analysis’, J. Comput. Struct., 33, 801–820 (1989). 11. A. Cardona and M. Geradin, ‘Numerical integration of second order dierential-algebraic systems in exible mechanism dynamics’, in M. F. O. Seabra Pereira and J. A. C. Ambrosio (eds.), Computer-Aided Analysis of Rigid and Flexible Mechanism Dynamics, NATO ASI Series, Vol. 268, Kluwer, Dordrecht, 1994, pp. 501–529. 12. A. Cardona and M. Geradin, ‘A beam nite element nonlinear theory with nite rotations’, Int. j. numer. methods eng., 26, 2403–2438 (1988). 13. A. Kubo, K. Yamada, T. Aida and S. Sato, ‘Research on ultra speed gear devices (reports 1–3)’, Trans. Japan Soc. Mech. Eng., 38, 2692–2715 (1972). 14. H. N. Ozguven, ‘A non-linear mathematical model for dynamic analysis of spur gears including shaft and bearing dynamics’, J. Sound Vib., 145, 239–260 (1991). 15. A. Cardona, ‘Flexible three dimensional gear modelling in mechanism analysis’, Europ. J. Finite Elements, 4, 663–691 (1995).

1/--страниц