Comparison and Implementation of a Rigid and a Flexible Multibody Planetary Gearbox Model

We propose algorithms for developing (1) a rigid (constrained) and (2) a (cid:29)exible planetary gearbox model. The two methods are compared against each other and advantages/disadvantages of each method are discussed. The rigid model (1) has gear tooth reaction forces expressed by Lagrange multipliers. The (cid:29)exible approach (2) is being compared with the gear tooth forces from the rigid approach, (cid:28)rst without damping and second the in(cid:29)uence of damping is examined. Variable sti(cid:27)ness as a function of base circle arc length is implemented in the (cid:29)exible approach such that it handles the realistic switch between one and two gear teeth in mesh. The (cid:28)nal results are from modelling the planetary gearbox in a 500 kW wind turbine which we also described in Jłrgensen et al. (2013).


Introduction
Of increasing interest are wind turbine gearboxes, because new wind turbines are large and gearboxes are one of the most expensive and critical components.
Many researchers/engineers consider not only the gearbox itself but also the system it is part of.By analysing the whole drivetrain, it is possible to simulate the dynamics of bearings, gear wheels, shafts and generators, i.e. stators and rotors.In Peeters et al. (2006) three types of multibody wind turbine drivetrain models have been implemented using the multibody software package DADS.It is found that adding bearing exibility is important for an accurate prediction of the eigenfrequencies.
The disadvantage of modelling a large system, i.e. the whole drivetrain or even the dynamics of the whole windturbine, is that the complexity increases.Modelling large systems should not be an issue, when us-ing most commercial multibody software tools.The commercial multibody software package Adams/View (MSC software) is used in Haastrup et al. (2011) to investigate best practices for gearbox bushings in wind turbine drive-trains.Other people combine their own gear programs into commercial tools.This approach is made in Rasmussen et al. (2012); Hansen et al. (2011) where an aeroelastic tool is used to model a bigger system, i.e. a whole wind turbine including exibility of the blades and tower with turbulent inow.
In this article we present algorithms that are implemented in Matlab.Focus is on the planetary gearbox without modelling the system it is part of, i.e. how the gearbox is connected to motors/generators etc.The models are presented with results using simplied and realistic loads.For the latter, we look at a planetary gearbox of a 500 kW wind turbine and use loads calculated by the aeroelastic and industry-accepted tool FLEX5 Øye (2001).For additional details of the experimental validation of this particular 500 kW gearbox, including the rigid multibody approach, see Jørgensen et al. (2013).

Multibody modelling
Introductions to the eld of planar and spatial multibody dynamics are given in Nikravesh (1988); Shabana (1989); Haug (1989); Géradin and Cardona (2001); Amirouche (2004).A mechanical system is dened in Nikravesh (1988) as a collection of bodies or links, which relate one body to the other.Bodies in these systems are connected by joints and/or force elements like actuators, springs and dampers.Mechanics is dened as the branch of analysis which involves the study of motion, time and forces on objects.The classical approach to solving the equations of motion for a constrained mechanical system is to solve equations for an index 3 Dierential Algebraic Equation (DAE) system (see e.g.Hairer and Wanner (1996)).The n dierential equations coupled with m algebraic equations is written Where M ∈ R n×n is the mass/inertia matrix, q ∈ R n is a vector of cartesian coordinates, Φ q ∈ R m×n is the Jacobian of the Φ(q, t) ∈ R m kinematic constraint equations, λ ∈ R m are the the Lagrange multipliers and g ∈ R n are the external body forces applied in the global reference system, aecting the origin of the local coordinate system of corresponding bodies.The local coordinate system is located in the center of mass, because this is a straightforward way to ensure a time-independent constant inertia/mass matrix.The vector Φ(q, t) is a symbol of general rheonomous constraints if at least one of the constraints is explicitly time-dependent.Geometrical constraints are used to x the gears of the planetary gearbox, hence Φ(q) = 0, which are scleronomous or time-independent.An example of a gear constraint for external gears is Φ gear ≡ ρ i θi + ρ j θj = 0 (2) where the radii for gear wheel pairs i and j are given by ρ and θ denotes the angular velocity of both gear wheels Nikravesh (1988).However, this constraint does not provide any realistic reaction forces because there's no coupling to translational degrees-of-freedom (DOF).
Both methods described in this article provide a realistic coupling between rotational (for transforming the torque between bodies) and translational DOFs.
Various techniques for solving (1) were developed and are described by many authors in e.g.Nikravesh (1988); Ascher and Petzold (1998); Jalon et al. (1995);Gear (1971).Constraints are expressed at acceleration level, to form equations of motion from the DAEsystem (1) into an ODE in the form given by Nikravesh (1988): By a proper selection of initial position, velocity (q, q), constraints (Φ) and forces and moments on the right-hand-side (RHS), the equation system (3) can be solved for the accelerations q.A time-integration scheme such as e.g.RungeKutta can then yield the velocity-and position-vector of coordinates (q, q) for all timesteps.

Methods
Both rigid and exible multibody dynamics are specifically concerned with analysis of bodies undergoing large displacements (both rigid body motion and elastic deformations).For general formulations of exible multibody systems see Shabana (1997); Géradin and Cardona (2001) and references therein.
The following illustrates how to calculate the length between points from two dierent objects/bodies, using two dierent local coordinate systems i and j: where r is a vector in the global system to the local coordinate center, A is a transformation matrix from local to global coordinates and s is a vector in the local coordinate system.The primed vectors are all in the local coordinate system, e.g.ω and vectors without a prime, e.g.ω, r or s are dened in the global coordinate system.

Connecting bodies, rigidly and exibly
With a kinematical description of a joint, i.e. using a constraint formulation, no exibility exists in the joint and the joint acts as a rigid connection.But if the constraint (Figure 1a) is removed and replaced by a force pair (Figure 1b) then the eect of how e.g.springs and dampers aect the two connecting bodies can be programmed.
The mathematically most simple joint between bodies i and j is a spherical joint, which is illustrated in Figure 1a and described as Φ s : r i + A i s i − r j − A j s j = 0.The constraint equa- tions ensures that two points measured from body i and j coordinate systems coincide because l = 0. From optimization theory, Farkas lemma guarantees the existence of Lagrange multipliers λ Bae and Haug (1988) which can also be interpreted as ctive forces or moments.The spherical joint can be transformed into a exible connection by removing the joint and replacing it with a force element.
x y ) Joint removed and added a force element (spring, damper or actuator).Bodies i and j are pulled towards each other by the forces F i and F j , acting at a distance of s i and s j .
Figure 1: Illustration of dierence between a rigid and a exible connection for bodies i and j.
When a rigid multibody model is changed into a model where some of the constraints are formulated by use of force elements (springs/dampers/actuators), the idea is to remove constraints from Φ and instead formulate equations that explicitly describe the forces or moments acting on the aected bodies.This reduces the number of rows in Φ q and λ.Replacing constraints with force element contributions means that the g-vector in (3) is modied.For instance body i could have an external force vector of where f x and f y are the sum of horizontal and vertical external body forces, while n is the sum of moments.A force element as shown in Figure 1b can be inserted by changing the external force vector, e.g.
could be illustrated by either inserting a translational actuator, a spring or a damper.The forces F i and F j are moved to the center of the local reference of frame together with the torque which is calculated as: Extending the formulation to include element/body exibility is described in e.g.Pedersen and Pedersen (1998); Shabana (1989).Another description of exibility eects in multibody systems is given in e.g.Huston and Wang (1994), which describe approaches for connecting rigid bodies for elastic straight and taper segments subject to extension, torsion and bending.

Implementation details of using force pairs
If a spring is inserted between the two bodies, the corresponding joint constraint should obviously be removed.The individual extra RHS force contribution to the g-vector in (3) becomes the scalar value of ∆F =k(l − l 0 ), where l is the deformed spring length cf.Equation (4), l 0 is the undeformed spring length and k is the spring stiness.Body i and j would sense the opposite force of each other as shown in Figure 1.Revolute joint connections between bodies can be transformed into exible damper/spring connections using the approach described above.The damping force is ∆F = d l where d is the scalar damping coecient and the derivative of the length calculated from Equation ( 4) is where: B = − sin θ − cos θ cos θ − sin θ

Constraints used for both models
Consider Figure 2 which shows the two fundamentally dierent approaches in this article for constructing a multibody planetary gearbox model with realistic reaction forces.The rigid approach works using Lagrange multipliers or constraints and the exible approach works by using linear springs and dampers with details described in section 2.3.Linear in the sense that the gear tooth penetration depth is measured along a line acting normal to the tooth surface.The gear tooth normal, radial and tangential forces are shown in Figure 3a, which also shows the path of contact as the tangent between two base circles (the green line).When two gear teeth mesh the distance from gear centers to anywhere on the contact line, varies along the yellow line.The gear loads depend on the transmitted torque and base radius.
Both illustrations in Figure 3 show that as the contact point (red dot in the illustrations) moves, the radii to the contact point continuously change.Despite that the radii of the contact point change, a fundamental property of involute gears is that the normal gear tooth force is always the torque divided by the base gear ra- The eect of vibrations, size of teeth, material errors, lubrication, number of teeth in mesh at the same time etc Klit et al. (2009) have no inuence in this model.

Description of rigid gear constraint
The problem when modelling a planetary gearbox with the simple gear constraint (2) is that only the torque is transferred through the rotational coordinates.There is lack of gear tooth forces, thus no physical reaction    F n , F r and F t , pressure angle α etc.
force contribution exists in λ, only torque.By considering an idealized gear, the constraint ( 2) is changed into the following non-holonomic gear constraint where v r is a unit vector between the center of both gears (see Figure 3b), A n is a transformation matrix (see eq. ( 4)) that depends on direction of rotation (CW/CCW) such that the matrix product (A n v r ) becomes a unit vector in the tooth surface normal direction as shown in Figure 3b.The gear constraint ( 6) is a dot product between a projection direction vector and a velocity vector of the red tooth contact point shown in Figure 3.This constraint equation gives a reaction force contribution because it includes all translational and rotational coordinates of involved bodies.It also has a physical meaning in the tooth surface normal direction.The constraint equation ( 6) is written completely out for 2D: where the unknown translational coordinates are x, y and the rotational coordinate is ω while vr is a rotated tangential unit vector.If the tangential It is easy to see that with DOF for ( ẋ, ẏ, ω), all DOF's are involved, which is important because it means that the internal reaction forces are treated correctly (see e.g.Jørgensen et al. (2013)).

Equations of motion
Equation ( 3) is solved for the accelerations q in each timestep and it is seen that Φ q q = γ.The Jacobian submatrix Φ q holds coecients for the acceleration terms.It is therefore necessary to collect coecients for q into the Jacobian (Φ q ).By using v = ω × v, Equation ( 7) dierentiated becomes Equation ( 8).
From this equation, LHS coecients to the Jacobian matrix and RHS contributions to γ can be found.

Flexible model
All the rigid model constraints are kept in the exible version of the multibody code except for the gear constraints that are removed.Body centers still need to be centered and xed to each other correctly.This means that e.g. the planetary carrier is still xed to planet gear wheel centers and sun/ring gear wheel centers are

Description of force element algorithm in the exible model
For every planet in a planetary gearbox (e.g.see Fig- ure 2), g s/p contain the sun/planet gear tooth forces and moments and g p/r holds the planet/ring gear tooth forces and moments.Every pair of sun/planet and planet/ring forces and moments consists of spring and damper forces/moments such that g = g 0 + g k + g d where g 0 is the initial body forces and torque, g k is the forces and torque from springs and g d is the force and torque from dampers (in the global reference of frame).
The total vector of external forces and moments for all bodies is modied such that g = g 0 + k g s/p + g p/r for k number of planets.This approach makes it easy to add/remove planets in a systematic way because the algorithm is generic and modular.
From Figure 4 it is seen that for a 1D spring and damper, forces are proportional to either a tangential distance interpreted as the penetration depth measured from the equilibrium position or to an angular veloc-ity dierence (the time derivative dierence in rotation speeds).These equations have to be changed a bit, in order to calculate the forces in higher dimensions, i.e. a transformation matrix is added.A 2D model is sufcient to illustrate the concept as this can easily be implemented in a 3D model.First the tangential tooth penetration depth is found, then the time derivative dierence is found and nally the spring and damping coecients are multiplied and the force contributions are added together.

Calculation of forces and moments
The driver/carrier angle θ d is the angle from the global x-axis to the local x-axis of the coordinate system for the rst planet.This angle has a corresponding rotation matrix A d .Similarly, A s/p (φ) is a transformation matrix for the angle from the center of the sun to the center of the current planet in mesh.Example: For the planetary gearbox in Figure 2, every sun/planet angle (insert the angles in eq. ( 4) to obtain each of the transformation matrices).There will be a loop that sums up all force contributions for all involved planets.The gear tooth penetration depth l p is now given by the projection where the unit tooth surface normal vector is: which is illustrated together with the pressure angle in Figure 5.In equation ( 9), θ 0 is the initial angle at t = 0 seconds.Normally θ 0 = 0 • and if this is not the case it means that the gear will start rotated in a position that is not the equilibrium position.Additionally, ρ b is the base circle radius for the involved body.The penetration depth is associated with the stiness properties (constant/linear/non-linear etc) of the gear(s).
The velocity dierence is associated with the damping properties of the gear(s) and can be written The absolute forces in cartesian coordinates are proportional to the spring/damping coecients so F k = kl p and F d = d lp .Reaction forces are opposite of each other for bodies i and j.
is the same vector rotated 90 • CCW and the tangential direction is expressed as a function of the radial unit vector: v t = vr .Therefore the scalar tangential gear tooth force is the dot product: and the absolute torque is M = ρ w F t where ρ w is the working circle radii.For gear wheels i and j the change on the The equations ( 12) are calculated twice as many times as there are planets, i.e. both for the sun/planet and for the planet/ring connections and all the force contributions are added together.

Stiness model with variable (multiple) teeth in mesh
The above assumptions rely on a simplication, where the multibody code knows nothing of the current number of gears in mesh since the stiness is given as a con- The base circle arc length s 1 is also shown in Figure 6.
In the present paper we use the same method as described in Pedersen and Jørgensen (2014) We also need to relate the change in angle ∆θ 1 to the change in base circle arc length ∆s 1 of the contact point, using the simple linear relation An oset-adjustment ∆θ a is the relationship between the mean angle and the stiness lookup-values.
The nal relationship between gear angle and base circle arc length is where s 1min is the minimum value of the arc length for the specic set of gears.

Without prole shift
For the given conguration used here the oset of the angles is the same for all planets.
The osetadjustment for the sun/planet is: The oset-adjustment for the planet/ring is:

Including prole shift
Inclusion of prole shift, x 1 , means that all planets interact with the sun and the ring at dierent places, i.e. the discontinuity in the stiness does not happen at the same time for all three planets (see Figure 12a).
An individual oset for all gear contacts is: • Sun/planet:

Gearbox with 1-4 planets
The two types of models described above are implemented and compared against each other.The number of bodies is equal to the number of planets plus the planetary carrier and sun and ring gears.shows that the model works as expected, with a cyclic harmonic motion.

Gear tooth forces -rigid vs exibly-connected model (1-4 planets)
For the rigid gear constraint, the reaction forces are already in the gear tooth surface normal direction.For a dynamic system the term −Φ T q λ yields the reaction forces but only the gear tooth reaction forces are interesting.The exibly connected method does not have Lagrange multipliers for the gear tooth reaction forces and therefore it is necessary to implement some book-keeping that tracks the penetration depth and penetration velocity if damping is included.
Figures 9a 9d show the gear tooth forces, for producing the same results as provided in Table 1c.The gear tooth force from the rigid model is compared with the exible gear tooth forces, given by the expression The eect of using only a single planet is easy to understand (see Figure 9a).It can easily be seen that after 2 seconds the carrier torque is removed and left is only a simple harmonic motion which without damping will continue forever (this innite harmonic behavior will become even more clear when the stiness is increased by a factor 10 in the following section).By using numerical integration in Matlab from 05 seconds, it can be seen that the area under the rigid/exible sun/planet curves deviates around 5%.The area under the planet/ring curves deviate approx 1%.In other words, the antiderivative of the gear tooth force curve is approximately the same in either case.
When using 2 or more planets (gures 9b-9d), it is more dicult to assess the inuence of how the rigid gear constraint behaves in comparison with the exible gear connection.In the rigid case, the rst planet behaves dierently from the other planets while in the exible case all planets behave equally.This is true as long as no prole shift is used and as long as the initial rotation error is inexistent.Planets 2-4 have 1 gear constraint on the sun gear as well as on the ring gear, causing a single gear tooth reaction force two places.
The exible model yields gear forces everywhere but they are the same when no prole shift is used and with the same perfectly aligned initial conditions.For this reason two subplots are shown for the simulation with 2-4 planets because the subplot for planet 2 is the same for all following planets.
When 2 or more planets are used, it is easy to realize that the exible gear connection is highly recommended in favor of the rigid gear constraint.Physically the exible method is more correct when comparing with the rigid method.However, the latter is easier to implement, it is faster if the stiness is high and realistic and it does not require any book-keeping at all.

The eect of damping
The stiness is now increased by a factor 10 so k = 5000 N/m and the number of planets is chosen to be 3 which   is deemed to be a typical case.Dierent damping coecients are added but prole shift is not used and perfect initial rotation angles as described in the previous section are used again.1c shows that as the stiness is increased the relative dierence between the two approaches decreases and in addition the amount of damping has little inuence on the relative movement of the carrier angle after 5 seconds.From basic theory of vibrations it is known that for a 1-DOF spring-mass system, the natural frequency is

Table 10a in comparison with Table
If the mass is constant and the stiness increases by a factor 10 the natural frequency is expected to increase by a factor of √ 10 = 3.16 which also seems to be the case.We can estimate the half period length with k = 500 N/m from Figure 9a-9d and it is approximately 1.75 seconds.From gures 10b-10d the half period length is approx 0.55 seconds.The ratio between the two is very close to 3.16, so the theoretical decrease in period length is conrmed as we increase the stiness by a factor of 10.
After 2 seconds, Figure 10b shows the innite harmonic motion of the gear teeth penetrating each other.
As long as no damping is included, the gear tooth forces are oscillating to an amplitude given by the penetration depth at t = 2 seconds, which is the moment where all external forces and moments are removed.Only the inertia of the system keeps the planetary carrier rotating and since no damping is involved, the energy of the system is preserved which can be seen by the regular periodic oscillations.

Figure 10d illustrates the equilibrium condition
where the stiness and damping forces locks up the gear teeth so after a while, the gears will end up moving in a steady and similar way as the rigid model.
With damping gear teeth can rotate while penetrating each other (otherwise the spring forces would be 0).
By analyzing the forces in the program, it can be seen that at exactly 2 seconds the gear teeth penetration depth is not large, but due to inertia, the spring forces are at a maximum around 2.2 seconds.The damping forces are almost constant but decrease slightly from 2.2 to 5 seconds.The result is that the sun/planet and planet/ring gear tooth forces end up at around ±2.3 N while the planetary gearbox rotates as steady as in the rigid model with small (or no) acceleration jumps.
Equilibrium with no external forces mean that if the sun/planet force is 2.3 N then the planet/ring force should be -2.3N and if this is not the case, the planet will rotate in a non-ideal way which is seen as a kind of rotation, which the rigid model would never do.This behavior of using damping coecients was fully expected.The bigger the damping coecients is, the more the exibly-connected planetary gearbox will behave like a rigid model and this is something that must be taken into account for real practical examples like those described in the following section.Figure 11 shows some realistic input loads that to some extent have been veried by using the aerodynamic FLEX5 software tool Øye (2001) and comparing with experimental data Jørgensen et al. (2013).
Because the stiness for this model is much higher than before, the timestep has to be much smaller, which increases computation time and consequently the total simulation time is rather short.It must however be long enough to be able to remove the rst part of the simulation results due to turbulent (wind) inow, giving the model time to nd something that can be regarded quasi-steady.In other words, it takes a while before inertia forces and accelerations are in a semi-equilibrium state.

Results from the 500 kW planetary gearbox model
In order to save computation time the initial velocity of the planetary carrier, connected directly to the rotor, is set to 2.86 rad/s or 27.3 RPM which can be seen from FLEX5 simulation results.The torque on the carrier and sun is then interpolated using Figure 11 and in the beginning of each timestep the g-vector (1) on the RHS of the equations of motion is updated accordingly, causing some acceleration uctuations due to turbulence.
Table 2 is the mass matrix for our realistic 500 kW wind turbine.As the gears rotate, the base circle arc length changes and because prole shift is used, the switch between 1 and 2 gears in mesh will happen at (a) Rotor/main shaft torque in a real 500 kW wind turbine in a wind eld with a mean speed of 10 m/s.   Figure 14 shows a comparison of the rigid vs. the exible approach for this particular planetary gearbox using two dierent damping coecients.It is interesting to note that the mean values are almost identical and approximately 1/3 of the rigid gear tooth forces.This corresponds fully to the fact that Figure 2 only shows 1 sun/planet gear constraint with the rigid model (because no more constraints can be added).By experimenting it can be seen that higher damping de-creases the amplitude of the gear tooth forces, because then gear teeth are not allowed to penetrate each other to the same extent.
Finally it becomes interesting to look at the torque from accelerations and moment of inertia of e.g. the sun, shown in Figure 15.It can be seen that the standard deviation of the torque becomes small the closer we get to the rigid case, which is equivalent to increasing the damping.In other words, damping prevents gear tooth penetration.

Conclusions
Two fundamentally dierent approaches for modelling multibody planetary gearboxes have been described and implemented in Matlab.The rst approach is a completely rigid model, where bodies are connected through a rigid gear constraint that allows reaction forces easily to be extracted, from the vector of Lagrangian multipliers.
The second approach also uses rigid bodies, but the gear constraints are exchanged with springs and dampers.Comparing the two approaches shows that gear tooth reaction forces must be in good agreement between each other because when running a simple example with 1-4 planets and no damping, the total angular movement of the carrier seems to be the same for all cases.The deviation of the position of thenal carrier angle after 5 seconds is found to be less than 1 %.These results could probably also be obtained by using a simpler rigid gear constraint such as

Desc. F [kN]
Rigid p 1 /r  Figure 14: Gear tooth forces using rigid/exibly connected bodies for modelling a real planetary gearbox using prole shift and realistic stiness/damping parameters.Final carrier angle is 409.622• by using the rigid approach, meaning that the dierence is negligible.
φ i θi = −ρ j θj from Nikravesh (1988), but then it would not be straightforward to obtain the gear tooth reaction forces, which is the scope of the present paper.
As we increase the stiness, the natural frequency behaves well accordingly to the theory of harmonic motion for springs.By adding more and more damping, the gears in the planetary gearbox clearly rotate more and more like they would do, if implemented using a rigid gear constraint.Overall the results of the two methods are described as: • The rigid gear constraint is generally very fast, but experiments with 1-4 planets show that the gear tooth forces are only a rough estimate of the real planetary gear tooth forces.A major drawback with this method is that only one of the planets can have a gear connection to both the sun and the ring, meaning that not all planets will feel the same gear tooth forces when operating under the same ideal initial conditions.Also fatigue calculations are inaccurate because the force variations are too small.
• The exibly-connected method is more realistic and should therefore be used for accurate calculations.However, a major drawback is that by using realistic stiness and damping parameters the time-step decreases signicantly and therefore the simulation time becomes longer than by using the rigid approach.The method is suitable for fatigue as well as other detailed calculations.
For the exible method we describe a method which implements gear tooth stiness depending on the rotation angles of the two gears in mesh.This method is also capable of taking prole shift into account, for making more realistic simulations.Intuitively it makes sense, that the mesh stiness is higher when 2 teeth are in mesh compared to when only a single tooth is in mesh.Simulations show that with no prole shift and perfect initial conditions, such as initial gear body rotation angles, all planets experience the shift between 1 and 2 gears in mesh at the same time.But with prole shift, the tooth clash when switching between 1 and 2 teeth in mesh, is shifted such that it happens at dierent times for dierent planets.
Finally, we present simulations for a model of a real planetary gearbox in a 500 kW wind turbine.For the planetary gearbox with 3 planets in the 500 kW wind turbine, the gear tooth forces using the exible model are close to 1/3 of those by using the rigid model.We demonstrate that the standard deviation of the moment from accelerations and moment of inertia becomes smaller the higher damping we use.Two fundamentally dierent approaches are theoretically vali-dated and usable for planetary gearbox modelling with any number of planets, using multibody dynamics.
Appendix (3D model) The skew-symmetric matrix of a vector is a 3×3-matrix and for ω it is ω.A useful method of calculating the cross product between two vectors is to multiply the skew-symmetric matrix of the rst vector with the second vector, e.g.ωs = ω × s.The skew-symmetric matrix for e.g. the vector of local angular velocities ω The 3D version of Equation 7 is: Φ3D : (A n v r ) T ri + ρ i (A i ω i ω i ) × v r + (A i ω i ) × v r + (A i ω i ) × vr −(A n v r ) T rj − ρ j (A j ω j ω j ) × v r + (A j ω j ) × v r + (A j ω j ) × vr Equation ( 26) dierentiated becomes Equation 27and a maximum 1 skewmatrix above mathematical symbols is used for cross products.From this equation, LHS coecients to the Jacobian matrix and RHS contributions to γ can be found.
Rigid planetary gearbox model with 4 gear constraints and 1 rotational DOF.Planet 1 transfers torque to/from the sun gear.Planets 2 and 3 can only have 1 gear constraint (either to sun or to ring gear) to prevent an over-constrained system.planetary gearbox model all planets are active.For illustrative purposes, the spring and damper is shown in the tangential direction.This is merely a projection of the real spring tooth forces which work in the tooth normal direction.
b) Radial direction unit vector (seen from the left gear) vr and two rotation-direction dependent unit vectors Anvr (CW/CCW).

Figure 3 :
Figure 3: Illustration of gear tooth force components

Figure 4 :
Figure 4: Illustration of spring k g and damper d g between two gears for a simplied spring/damper system.Gear tooth force is perpendicular to teeth surface.

Figure 5 :
Figure 5: Surface normal vector in the penetration direction as a function of the pressure angle α.
stant.Improved realism and accuracy of the exiblyconnected model can be obtained by letting the stiness be a function of the rotation angles of the gears.The stiness coecients of gear teeth are generally non-linear and a function of the point of contact and the size of the contact load, see e.g.Pedersen and Jørgensen (2014).The meshing stiness is the combined stiness of two gears in contact.The number of gear teeth pairs in contact is also changing, usually between one or two pairs in contact.This leads to discontinuities in the meshing stiness.The point of contact on a gear tooth can be given as a function of the base circle arc length as shown inPedersen and Jørgensen (2014).
for extracting the stiness as a function of the base circle arc length s 1 .The FEM package Comsol MultiphysicsMultiphysics (1998Multiphysics ( -2012) )  is used for obtaining the stiness coecients.Other methods for including tooth stiness are found inEbrahimi and Eberhard (2006);Ziegler et al. (2006).

Figure
Figure7aand 7b show respectively the sun/planet and planet/ring stiness.As soon as two teeth are in mesh at the same time, the stiness almost increases to the double.The gures can also be used to calculate the contact ratio, i.e. the average number of teeth in mesh at the same time.It can thus directly be seen that the planet/ring contact ratio is higher than the sun/planet contact ratio.Introducing this kind of stiness into the multibody program causes small shock impulses which can decay fast or slow, depending on the amount of damping added.In order to obtain stiness expressions we need to relate the rotation of the involved gears to the dened stiness variations.The gear rotations are dened by the variables θ 1 and θ 2 .Since these two rotations are unconstrained we dene a mean value.The stiness is dened relative to the base circle arc length of gear number one, i.e. a mean value θ1 is dened as the mean value of the actual θ 1 and the value θ * 1 (the latter being the value that θ 1 would have relative to the θ 2 value if a gear constraint was used instead):

Figure 6 :
Figure 6: Base circle arc length, s 1 and s 2 with involute curve shown on gear 1.The load from gear 2 acts along the line of action.Base circle radii are denoted r b and working pressure angle is α w .Angle from center to contact points are denoted β.

Figure 7 :Figure 8 :
Figure 7: Stiness as a function of the base circle arc length s (cyclic dependence).

Figure 8
Figure 8 is shown with 3 planets although there is space for 4 planets.Input torque is put on the planetary carrier i.e. the driver.The carrier angle θ d is measured from the global x-axis to the local x-axis of the carrier.This local carrier x-axis is aligned in the direction of the center of the rst planet.The rst planet always starts at 0 • , then the next planet rigid/exible comparison.The rigid gear constraint provides a good average for gear tooth force.rigid/exible comparison.As shown in Figure 2, only planet 1 has both gear constraints at sun/planet and planet/ring.rigid/exible comparison.It can be seen that the rigid gear constraint is a rough inaccurate estimate of the exible gear tooth forces.

Figure 10 :
Figure 10: Gear tooth forces for the rigid vs. exibly-connected articial planetary multibody gearbox model for dierent damping coecients.

3. 2
Figure 7a (sun/planet) and Figure 7b (planet/ring) is used in the following, for modelling a planetary gearbox with 20, 35 and 91 teeth for respectively the sun, planet and ring gears.The prole shift is x s = 0.582, x p = 0.419 and x r = −0.840,respectively and for input, realistic loads from Jørgensen et al. (2013) are used.The rotor torque is used on the planetary carrier and the generator torque is used on the sun gear, i.e. their values are inserted to the RHS of Equation (3).
variation of one planet due to switch between 1 and 2 teeth in mesh, based on (a).

Figure 12 :Figure 13 :
Figure 12: 6-body multibody model and demonstration of realistic stiness model using prole shift.If the base circle arc length is below a threshold value, 2 teeth are in mesh and otherwise only a single teeth is in mesh.

Figure 15 :
Figure 15: Inertia term for the sun gear (M sun = I zz θ) using dierent damping coecients, for the realistic planetary gearbox in a 500 kW wind turbine.

Figure 2
Figure2shows two versions of a 6-body multibody model, i.e. consisting of sun, ring, 3 planetary gears and a planetary carrier.With 6 coordinates per body (3D) there are 36 DOF.All 6 bodies are xed at their mass center to either ground or to the planetary carrier.In 3D it require 30 constraint equations but as the ring gear cannot rotate and with 4 gear constraints, it means that there is one rotational DOF left.When torque is applied to either body, everything will rotate appropriately.The 3D version of Equation2.1.1 becomes: is easy to see that with DOF for ( ẋ, ẏ, ż, ω ξ , ω η , ω ζ ), all DOF's are involved which again is important for treating the reaction forces correctly (see e.g.Jørgensen et al. (2013)).