Gear fatigue damage for a 500 kW wind turbine exposed to increasing turbulence using a flexible multibody model

This paper investigates gear tooth fatigue damage in a 500 kW wind turbine using FLEX5 and own multibody code. FLEX5 provides the physical wind (cid:28)eld, rotor and generator torque and the multibody code is used for obtaining gear tooth reaction forces in the planetary gearbox. Di(cid:27)erent turbulence levels are considered and the accumulated fatigue damage levels are compared. An example where the turbulence/fatigue sensitivity could be important, is in the middle of a big wind farm. Interior wind turbines in large wind farms will always operate in the wake of other wind turbines, causing increased turbulence and therefore increased fatigue damage levels. This article contributes to a better understanding of gear fatigue damage when turbulence is increased (e.g. in the center of large wind farms or at places where turbulence is pronounced).

The eciency of wind turbines is continously being improved in order to decrease the total cost of energy and therefore it is important to be able to in-vestigate the fatigue damage under dierent circumstances.One way of dealing with fatigue problems is mechanical: Changing the design, increasing safety factors, dimensions etc.
Another method is by developing or improving wind turbine controller mechanisms.Wind turbines are complicated machines with a lot of electrical equipment such as controllers for automatic adjustment of e.g.pitch/yaw regulation.In Molinas et al. (2010) an indirect torque control (ITC) technique has been investigated and the idea is that electromagnetic torque transients caused by grid faults and disturbances results in signicant gearbox fatigue.
As demonstrated in Schlechtingen et al. (2013) there is also a lot of development going on in the eld of condition monitoring which include recording large amounts of data using sensors for e.g. oil temperature, wind/rotor/generator speed, vibration/accelerometers, current and power output etc., which is then processed for early diagnosis of problems.
Fatigue loads for rotor and main gearbox bearings are calculated in Heege et al. (2007) using a method that couples non-linear Finite Element Method (FEM) and super element technique with a multibody approach (see e.g.Haug (1989); Géradin and Cardona (2001)).An approach to predicting wind turbine gearbox reliability for three generic gearbox congurations, based on assumptions of estimated failure rates, is given in Smolders et al. (2010) The most common cause of gearbox failure is surface contact fatigue (see Fernandes and McDuling (1997) which also presents examples and describes rolling contact fatigue, sliding-rolling contact fatigue and spalling).
Recently, researchers from National Renewable Energy Laboratory (NREL) initiated the Gearbox Reliability Collaborative (GRC) and guarantee privacy of commercially sensitive information, see Musial et al. (2007) in order to investigate gearbox problems and seek solutions that will lead to higher gearbox reliability.
A description of four stages of fatigue damage from initial crack to failure is given in Vasudevan et al. (2001), which also provides an overview of the topic from a historical perspective together with suggestions for the concept of fatigue damage control.
A typical procedure for fatigue load data analysis for wind turbine gearboxes is given in Niederstucke et al. (2003).In this article the commercial software package FLEX5 Øye (2001) will be used for obtaining wind turbine loads (forces and torque).The focus will be on investigating gear tooth fatigue, using a multibody dynamics program developed in Matlab.
The gear tooth stresses are found using either Comsol Multiphysics Multiphysics (1998-2012)  (2) a modied maximum likelihood method and (3) a graphical method.The Weibull cumulative distribution function Seguro and Lambert (2000) is given by where P (v ≥ 0) is the probability that the wind speed equals or exceeds zero while P (v < v i ) is the probability that the wind speed is below v i .Weibull parameters for this site is estimated to A = 6.operation of 20 years.
For the multibody program, a similar procedure as in Dong et al. (2012) with topology diagram shown in Figure 2 will be employed.Instead of using data from the NREL (GRC) 750 kW wind turbine Musial et al. (2007), data from a 500 kW planetary gearbox will be used.For aerodynamics FLEX5 instead of FAST is used.Our own multibody code instead of SIMPACK is used and rainow counting and Palmgren-Miner summation is made with Matlab.

Wake eects and turbulence in wind farms
The inow to wind turbines is always turbulent.This puts a special demand on the modelling of the incoming wind, as it is required to model not only the mean properties of the ow eld, but also the turbulence charac-

Fatigue damage
Employing SN or Wöhler curves is a suitable approach that is generally accepted for fatigue lifetime estimation.The idea of the chosen method is to use a rainow counting algorithm on the gear tooth stress time-series.
The result is a measure of cycle amplitudes or ranges which can be counted in bins and collected in a rainow matrix, which is a function of mean values and ranges.See e.g.Madsen (1990) for a description of the rainow counting algorithm, which together with the linear Palmgren-Miner rule, based on Palmgren (1924) and Miner (1945) forms the following important relationships Equation (2) states the accumulated fatigue damage, also known as the linear damage hypothesis.Generally the material/component will fail when the sum D >= 1.At stress level number i, the number of counted cycles n i is a result of the rainow counting algorithm evaluated on the time series.The algorithm returns the range (or amplitude) of all stresses for the i th case, i.e. (n i , s i ).Capital N i is the maximum number of cycles to failure and also shown in Accurate information about the SN-curve for huge wind turbine components like the gears described in the present article, is assumed to be known.The gear root bending stresses are obtained using Finite Element Modelling (FEM) using Comsol Multiphysics Multiphysics (1998-2012).Surface stresses are considered as well and calculated using Hertzian contact stresses Norton (2000).

Planetary gearbox model
The following section describes planetary gearbox details (tooth data and facewidth etc.) followed by a description of the implemented multibody code for obtaining gear tooth forces and stresses.
The number of teeth in the rst (planetary) gearbox stage is 20, 35 and 91 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.Minimum contact facewidth for sun/planet and planet/ring is 210 mm.
It is assumed that bearings and wind turbine model components are rigid such that the gear tooth forces depend on axial rotor and generator torque, regardless of transverse forces and torque (a 2D model).A real wind turbine and gearbox which has been running since 1995 is considered.It is therefore expected that the calculated accumulated fatigue, D in Equation ( 2), for 20 years of operation should be acceptable.

Multibody program code
Multibody dynamics software is used to couple or link mechanical rigid or elastic bodies with each other.
General introductions to the eld are given in e.g.Nikravesh (1988); Shabana (1989); Géradin and Cardona (2001).Equations of motion are expressed using the following n dierential equations coupled with m algebraic equations 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.For the model here the primary input is the load.At one end the rotor torque is applied and on the other end of the transmission system the torque from the generator is applied.The model load input is the output generated by the FLEX5 program.Constraint equations in the Φ(q, t)-vector are expressed at acceleration level, to form equations of motion from the DAE-system (4) into an ODE in the form given by Nikravesh (1988):  (2014a), where it is compared to a method using a constraint formulation instead.
The force element algorithm requires knowledge of the total gear tooth force, i.e. the sum of the gear tooth forces on either one or two teeth.Fatigue is a local phenomenon in single points so knowledge of the force, i.e.
magnitude and direction and where it is acting on a single tooth is required.The main specic points of interest are either at the gear root or at the tooth surface contact point.The dierence between obtaining and storing the force acting on a single tooth instead of the total gear force, can be regarded as a matter of calculating a mesh stiness ratio and then internally store the mesh stiness of a single tooth.This method is described in details in Jørgensen et al. (2014a).The force acting on a single tooth is found by using Hooke's law, by multiplying the tooth penetration depth and the stiness from the multibody code.

Details of constraint functions
A planetary gearbox with 3 planets as depicted in Figure 5 consists of 3 planets, 1 carrier, 1 sun gear and 1 ring gear or in total, 6 bodies.All the bodies can rotate except the ring gear which is xed.Six bodies which can rotate (θ) and translate (in the x− and y−directions) result in 18 generalized coordinates, which will also be the size of the q-vector.The following is a list of 6 bodies: 1. Sun gear.
In the mathematical denition of the constraints, r, describes the global position of the origin of the local coordinate system attached to a rigid body.A transformation matrix, A, transforms a vector from local coordinates to global coordinates.A prime is added to a vector to indicate that the vector is dened in local coordinates, while a vector symbol with no prime attached, represents a vector in the global/inertial coordinate system.The following lists the 13 constraint equations (see Figure 5): 2 constraints x the center of the sun gear to the ground (revolute joint): 2 constraints x the center of the carrier to the ground (revolute joint): 2 constraints x one of the endpoints of the carrier to the center of planet gear 1 using a revolute joint: 2 constraints x the ring gear to the ground.

Gear tooth exibility and applied forces
A exible 2D multibody program has been made, which unlike a rigid model has gear tooth penetration as illustrated in Figure 7a.The details of this code is described in Jørgensen et al. (2014a) and shortly summarized here.A transformation matrix is dened as and the penetration depth l p along the line of action for gear bodies i and j is where A d is a driver/carrier transformation matrix (using θ d shown in Figure 5), A s/p is a transformation matrix from the direction shown using θ d (Figure 5) to one of the three planets, e.g.0 • , 120 • or 240 • .The unit tooth surface normal vector v n is a function of the pressure angle (typically α = 20 • or α = 25 • ), illustrated in Figure 6 and dened as Gear tooth penetration is illustrated in Figure 7a.
The time-dierentiated penetration depth lp is associ- ated with damping and is calculated using Figure 6: Unit tooth surface normal vector in the penetration direction as a function of gear tooth pressure angle α.
With stiness coecient k and damping coecient d, gear tooth stiness and damping forces are calculated using F k = kl p and F d = d lp respectively.Reaction forces are opposite of each other for bodies i and j and the force is applied to the g-vector on the right hand side (RHS) of Equation ( 5).Additionally, these forces contribute to additional torque.A radial unit vector from center of sun to center of planet gear v r is introduced.This radial unit vector v r can be rotated 90 and the tangential gear tooth force is the spring and damping force sum projected into the tangential direction using the dot product F k+d Hence, the torque is M = ρ w F k+d t where ρ w is working circle radius.The total force and torque vector on the RHS of Equation ( 5) become The frictional force in the gear tooth contact is neglected in the multibody simulation due to the relative small size.The friction is the source of loss in the gear box, and also has an inuence on the contact stress in the gear tooth.The friction is therefore included in the contact stress evaluation.
2.4 Gear tooth stiness, base circle arc length and gear tooth forces      7b) and corresponding FEM-stiness.
package Multiphysics (1998Multiphysics ( -2012) ) and expressed as a function of an average angle or base circle arc length.
Figure 7a shows gear tooth penetration (strongly exaggerated for illustrative purposes) and the average of s i and s j , is made equal to s 1 in Figure 7b, which is the contact point of two completely rigid gears in mesh.
The FEM-model is used to calculate the stiness at dierent contact points (Figure 9).In other words, the FEM contact point is located at a distance of s 1 from r b1 (or at a distance of s 2 from r b2 ) and it is necessary A ctive single gear tooth force is required for fatigue analysis.The reason is that it is necessary to only look at one gear tooth at the time, even though sometimes the load is shared among two gear teeth.
The problem is extensively described in Pedersen and Jørgensen ( 2014), but shortly explained the single gear tooth force is calculated as F t k 1 /k t , where F t is the total gear tooth force (sum of spring and damping forces), k 1 is the single gear tooth stiness and k t is the total gear tooth stiness.
Figure 8 shows 0.2 seconds of simulation time and it shows whether 1 or 2 teeth are in mesh at a given time (high stiness is when two gear teeth are in mesh and low stiness of around 3.5•10 9 N/m is when only one gear tooth is in mesh).The gure shows that most of the time 2 gear teeth are in mesh, except in the interval [0.145; 0.17] seconds.
The method used here for simulating the exibility of the gear tooth contact has the advantage that it is rather simple, leading to acceptable computer simulation times.At the same time the highly non-linear nature of the contact is included.The frictional force due to the sliding between the gear teeth is neglected in the multibody simulation, but later included in the estimation of the contact stress.

Gear tooth stresses in relation to fatigue analysis
The gear tooth stiness properties is the same as in 2. Gear tooth surface failure: Hertzian stress assumption is applied and pure rolling contact of cylinder-on-cylinder is assumed.The maximum stress a little below the surface is calculated.
The details of the two approaches are elaborated below and the 8 cases illustrated in Table 1 are considered using increasing turbulence intensity (TI).

Gear tooth root bending failure
The gear tooth stresses are assumed to be linearly dependent on the magnitude of the gear tooth force obtained from the multibody program (from Figure 8).
In other words, the gear tooth stresses are assumed to be the normal single gear tooth force multiplied by the polynomial function describing the relationship between stresses and base circle arc length (see Figure 9).
Furthermore it is assumed that the maximum stress appears at the same point in the tooth root although this is not completely true, i.e. the point of maximum stress moves slightly.

Gear tooth surface stresses -pitting
Contact surface stresses are calculated as if two cylinders with radii equal to the base circle arc lengths are in contact.The cylinder length is the minimum gear facewidth, resulting in an ellipsoidal-prism pressure distribution as described in Norton (2000).Gear teeth are not only exerted to pure rolling but also to signicant sliding which changes the stress state.This quickly becomes a complicated problem and therefore the complexity is signicantly reduced by estimating the maximum von Mises stress using a simple average between a static and a dynamic loading coecient from Norton (2000).This can be justied because the averaging stress coecient is deemed not important to calculate exactly, since it is the same for all simulation A cylindrical geometry constant B is dened as where ρ 1 and ρ 2 are the base circle arc lengths of the two gears in mesh (also illustrated in Figure 10, i.e. corresponding to the values of s 1 and s 2 ).The contact-patch half-width a is illustrated in Figure 10 and can be calculated as The estimated Von Mises stress is taken as an average of a static and a dynamic coecient (from Norton (2000)) and it is estimated to σ = p max 0.57+0.73 2 = 0.65 p max .Because the relative fatigue damage is important, it is deemed that it is not that important whether 0.6, 0.65 or 0.7 is used but 0.65 can be justied, at least to illustrate the proposed method.The maximum von Mises stress is found subsurface Norton (2000) leading to fatigue micro-cracks and development of pitting.The von Mises stress is used together with a rainow-counting program that has been veried against a similar tool that is part of FLEX 5. • Section 3.2 describes some SN (Wöhler) curve considerations, i.e. material fatigue characteristics.

Adjusting time-factors due to dierent gear teeth in contact
Simulation time is 30 seconds (computation time is 10-15 hours per run).For fatigue-calculations the simulation time must be upscaled to the real number of operational hours given by the Weibull distribution for this particular site (shown in Table 1).Table 2 shows a simple (unadjusted/naive) time-factor calculation, i.e. to consider is that every full rotation of a tooth on the sun gear, is met by teeth from the three planet gears (see Figure 5).Therefore the fatigue from the calculated stresses for teeth in contact must be added together.Instead of having three 30-second simulations this situation corresponds to making a fatigue calculation on a longer 90-second simulation.The same idea for the sun gear can also be used for a full ring gear cycle, which sees 3 planet planet/ring gear connections.
One rotation of the planet gear however corresponds to only a 60-second simulation because one gear tooth will see another tooth from the ring and from the sun gear at all times.Tooth-adjusted time-factors are summarized in Table 3.

SN-curve details
Realistic SN-curves for these particular gear teeth are dicult to obtain and therefore the endurance limit as well as the slope must be estimated.It is known from Norton (2000) that the maximum von Mises stress is slightly below the surface which is therefore expected Table 3: Adjusted time scaling for teeth on dierent gears.Last three rows is the fatigue time-factor for a sun, planet or a ring gear tooth.Because the stresses (and gear tooth forces) are comparable, a fatigue calculation does not need to be done on all individual planets.
to decrease the endurance limit.SN-parameters (see Figure 4) will be estimated based on Hirsch et al. (1987); Boyer (1986) using the two points in Table 4.
The following expression is obtained (see Figure 4) where s 0 = 28.9GPa and m = 3.4, while σ e is the endurance limit in MPa.A large uncertainty is involved when choosing material and SN parameters, however for illustrative purposes the endurance limits in Table 5 are chosen.The endurance limit works as a threshold which eectively can neglect fatigue problems if raised just slightly to e.g.850-1000 MPa.Table 5 shows that the sun and planet gears have the same, relatively high, endurance limit.
Table 4: SN curve parameters using two points (P 1 and P 2 can be seen in Figure 4).
Sun gear 800 Planet gear 800

Ring gear 550
Table 5: Endurance limits for SN curve.
Adjusting the SN curve and endurance limits makes a great dierence in the results.If all parameters are the same (e.g.endurance limit of 550 MPa), results clearly illustrate that the sun gear is most exposed to fatigue.A plausible explanation is that three planet gear teeth are acting on the sun and the sun gear has only 20 teeth.Results also clearly show that the ring gear is the least subjected to fatigue.Although it is constantly exposed to contact from three planet teeth, it has a total of 91 teeth and hence the lowest timefactor.

Results
Results are divided into a section dealing with root bending failure and surface stresses (pitting).

Root bending failure
The maximum gear tooth stresses are expected using a mean wind speed of 20 m/s.However the results in Figure 11 using TI=10% indicate that the maximum stresses at the gear root is approximately 60 MPa, which is clearly below the endurance (or fatigue) limit.
Even though the maximum stress would be higher for increasing turbulence intensities, it is not considered a problem because the endurance limit is much higher.
It can be concluded that for these particular gear teeth, no fatigue damage will happen at the root.

Surface stresses / pitting
The maximum von Mises stress using Ū = 20 m/s and TI=10% is found to be 6-700 MPa which will contribute to fatigue damage.Figure 12 is an example of the planet/ring stress level for a 30-second simulation.We select a specic location on a gear tooth (i.e.choose a corresponding s-value, see Figure 7) and draw a vertical line using the results from Figure 12b and process the von Mises stress using a rainow-counting algorithm and the SN-curve data from Table 4.         From Figure 13 it can be concluded that the average stresses are lowest for low values of s (the two red curves).This is a region with high stiness due to the fact that two teeth are in mesh at the same time.Additionally, a black curve is added to both graphs in Figure 13.This line indicates the contact point of the two working circles, i.e. the pitch point where only rolling appears (everywhere else sliding takes place).

Average von Mises stress
At this point only a single tooth is in mesh and that explains the higher stresses.Figure 12 shows that the stresses are not the same everywhere, i.e. it is ambiguous which value of s yields the maximum average stresses for both sun/planet and planet/ring results.

Accumulated rainow counting results
Figure 14 shows an example of processing the 30 seconds sun/planet von Mises stress calculation results through the RFC-code.Figure 14a shows the stress range for each of the 6-20 m/s mean wind speed bins.
It can be seen that the change from 6-8 and 8-10 m/s is huge in comparison with the change from e.g.16-18 or 18-20 m/s. Figure 14b is the cumulative number of full cycles, using the time factor for the sun gear (TF s from

Cumulative fatigue damage due to increased turbulence
The IEC 61400-1 standard about design requirements for wind turbines specically operate with 3 classes for turbulence intensity: 12%, 14% and 16%.Turbulence intensity depends on e.g.surface roughness i.e. how many disturbances are present from e.g.
trees/buildings etc and also hub height.In other words, the turbulence intensity is generally lower for oshore wind turbines (undisturbed) when compared to that of landbased wind turbines.With increasing altitude, the ow is generally also less disturbed leading to less average turbulence.It is found that the investigated 500 kW wind turbine gearbox has no fatigue-problems related to gear tooth root bending.The Hertzian contact stress assumption implies calculating the maximum von Mises stresses slightly under the gear tooth surface.Fatigue damage occurs under the assumption of using an endurance limit of 800 MPa for the sun and the planet gear wheels.
The interpretation of this result is that fatigue just below the gear tooth surface can lead to gear tooth surface pitting, which eventually increases the risk of gearbox failure.Furthermore it can be concluded that: • The average von Mises stress attens out in the high mean wind speed region, i.e. in the interval of mean wind speeds approximately between 15-20 m/s, the stresses are not much higher than at 15 m/s.This behavior resembles that of a typical power-curve for a wind turbine and it also resembles a similar result published in Jørgensen et al. (2014b).
• The stresses at specic locations have been calculated using interpolation and it can clearly be seen that the von Mises stress on a single gear tooth in the region where two gear teeth are in mesh (low base circle arc lengths), are lower than in the region where only a single tooth is in mesh.
Specically it can be concluded that the stresses at the working circle/pitch point and for higher base circle arc lengths are higher, than at points specied by low base circle arc lengths.The gear tooth stresses on a specic place on the tooth are lower in the region where two gear teeth share the load, when compared to the case where only a single gear tooth must take the whole load.
• When estimating the eect of turbulence on gear tooth surface pitting, a marginal increase in turbulence in the high turbulence region relatively contributes much more to accumulated fatigue damage, than the same increase in low turbulence regions.
A practical perspective of these results suggests that for low turbulence intensities, wake eects (turbulence generated from other wind turbines in a wind farm) are not deemed important as a cause of gearbox failure due to pitting caused by subsurface cracks.
It is common to assume that the turbulence intensity for oshore wind turbines is lower than that of similar landbased wind turbines (because these are located near obstacles such as trees, buildings, uneven terrain etc which increase the turbulence).Our results suggest particular attention to pitting as a cause of gearbox failure, should be taken when making wind turbine layouts of wind farms in regions with high turbulence.
emergency stops caused by drivetrain problems are very expensive for the wind turbine industry.Drivetrain components are critical for operation, they're typically large, heavy and expensive to repair or replace.Costs for operation and maintenance (O&M) increase by the time wind turbines are non-producing.Component lifetime reliability is important.A key parameter to designing drivetrain components is calculation (estimation) of accumulated fatigue damage.Related problems during the lifetime of a wind turbine are issues with corrosion, wear and leading edge blade erosion, which can reduce the performance.

Figure 1 :
Figure 1 demonstrates the concept of this article, which is to calculate gear tooth fatigue based on turbulent wind input.Accumulated fatigue damage for a wind turbine producing power using wind speed from 6-20 m/s and a bin size of 2 m/s will be investigated (the 6 m/s bin contains wind speeds from 5-7 m/s etc.).Based on Jørgensen et al. (2014b) the turbulence intensity (TI) is set to around 10-20% throughout the whole interval.For this particular site the Weibull parameters based on 5 years of wind data from Hansen and Larsen

Figure 3 :
photo of the ow eld inside the Horns Rev wind farm.The photo, which was taken an early morning, shows the turbulence generation due to the mixing of wakes visualized through condensation of water drops.The rst turbines in the wind farm are seen to only subject to the turbulence of the incoming wind, whereas the remaining turbines are subject both to the ambient turbulence and the additional turbulence generated by the wakes of the surrounding turbines.This increases the overall turbulence level and causes the turbulence to become anisotropic.In the interior of the farm, the turbulence settles at a constant level and, due to mix-
5) Force elements (springs and dampers) are used to model the connecting exibility between rigid gears, instead of using (rigid) gear constraints.A full description of the applied force element algorithms as illustrated by Figure 2 is presented in Jørgensen et al.
additionally prevents the ring gear from rotating.Φ 9 : θ 3 = 0 (10) 2 constraints x one of the endpoints of the carrier to the center of planet gear 2 using a revolute joint.Φ 10−11 : r 4 + A 4 s 4,P 2 − r 5 = 0 (11) 2 constraints x one of the endpoints of the carrier to the center of planet gear 3 using a revolute joint.Φ 12−13 : r 4 + A 4 s 4,P 3 − r 6 = 0 s 4,P 2 and s 4,P 3 can be found using a transformation matrix where the dierence is that planet gears two and three are 120• and 240• apart away from planet gear number 1.The dimension l c is illustrated in Figure5.There are no gear constraints.Instead the forces that make the gear wheels rotate are evaluated by modifying the RHS of the equations of motion (5)

Figure 7
Figure 7 illustrates some important concepts in relation to gear geometry and the base circle arc lengths s 1 and s 2 .With exible connected bodies both gear wheels can penetrate each other.The gear tooth stiness has been calculated using a Finite Element Method (FEM) Base circle arc length, s1 and s2 with line of action seen from the point of view of gear 1. Base circle radii are denoted r b and working pressure angle is αw.Green arrows indicate 2 teeth in contact (high stiness) and the red arrow indicates 1 tooth in contact with the other gear.

Figure 7 :
Figure 7: Tooth penetration and estimation of equivalent stiness from exible multibody model.The sum of distances s 1 + s 2 = s c is a constant, see additional details in Pedersen and Jørgensen (2014).

Figure 8 :
Figure8: Top graph: The total gear tooth force is shown in blue, while the gear tooth force from a single gear (due to lower stiness) is shown in red.When only one pair of teeth is in mesh, the red and blue curves are coinciding.The multibody program uses the blue curve, but for the fatigue analysis, the red (single-tooth contact) curve will be used.Lower graph: The base circle arc length of the inner gear (s 1 on Figure7b) and corresponding FEM-stiness.
Figure7bshows that the base circle arc length (s 1 or s 2 ) represents the current rotation of both gear wheels as well as the gear tooth penetration depth of each gear tooth into the other gear(s) in mesh is a function of stiness and tooth force.The sum of s 1 and s 2 is constant, but both distances change as a realistic dynamic model of the gear teeth requires this.By means of the gear tooth penetration depth, the average base circle arc length and the gear tooth stiness, the total gear tooth force for all teeth is calculated and plotted in Figure8.This gure shows both total gear tooth Jørgensen et al. (2014a) and the damping coecient is estimated to d = 40 • 10 3 Ns/m.The following two stress situations have been examined 1. Gear tooth bending failure: A FEM-program has been used for obtaining stresses in the situation where a unit force is applied to dierent contact points (dened as the base circle arc length).
21) using material constants m 1 = m 2 = 1−ν 2 E where ν = 0.3 is the Poisson ratio and E = 210 GPa is the Young's modulus.The gear tooth force calculated by the multibody program is denoted F .The maximum pressure/stress is

Figure 9 :
Figure 9: Maximum von Mises stress at gear tooth for a unit force (1 N) applied at a given base circle arc length distance s.A second order polynomial t approximates the FEM-results well, for the base circle arc length intervals of interest.
when using the rainow-counting algorithm, the number of cycles must be multiplied by this number to get an idea of the accumulated fatigue damage over the lifetime.The time-factor is calculated as the operational hours times 3600 seconds, divided by the simulation time in seconds.It is assumed that everything within those 30 seconds of simulation time is representative for the whole lifetime.Due to the high stiness of the gear teeth the timestep is rather small and the simulation time limit depends on computation time and memory.The sun has 20 teeth, the planets have 35 teeth and the ring 91 teeth.In other words on average those 30 seconds of simulation time and corresponding von Mises stresses, correspond to 600 simulation seconds for one specic sun gear, 1050 simulation seconds for one specic planet gear and 2730 seconds for one specic ring gear.The time-factor will simply be reduced by the number of gear teeth.Another important thing

Figure
Figure 7a illustrates how the base circle arc lengths is measured and Table6indicates important base circle arc lengths (the middle value for s is the interface/switch-point between a single and two teeth
(a) 3 transparent plates are inserted, left region shows twoteeth contact and right region shows single-tooth contact.(b)2D plot of (a) including both regions (black vertical lines).Mean wind speed is 20 m/s, TI=10%.

Figure 12 :
Figure 12: Estimated von Mises stress for 30 second simulation time, as a function of base circle arc length s.For this to be plotted properly, interpolated von Mises stress values are used.

Figure 12
Figure 12 shows an example of the stresses calculated for a simulation with indication of 3 dierent base circle arc length locations.The multibody program only calculates the gear tooth forces as a function of the dynamics of the whole system and that system does not

Figure 13
Figure 13 resembles a behavior similar to those found using a rigid multibody model in Jørgensen et al. (2014b).It assumed that multiple exible planets on one hand decrease the gear tooth force because additional planets share the load.On the other hand, the gear tooth force is probably increased because not all springs and dampers are aligned equally (especially not when addendum modication is used).If one of the planets tries to rotate in one direction, one of the other planets could try to rotate in another direction due to the spring/damper force.With high stiness these spring and damper forces are not negligible and cause the three planets to act against each other while trying to rotate the sun and ring gears at the same Half cycle count, 6-20 m/s mean wind speed.

Figure 15 :
Figure 15: Accumulated fatigue on sun, planet and ring gears for lifetime operation (20 years).

Table 1 :
62 and k = 2.19 us-Topology diagram of the Matlab multibody gearbox model (as illustrated in Dong et al. (2012)).
simulations are rather short compared to the timescales in Table1.Therefore it is important to ensure the statistical behavior of the simulation input and output as simulation results will be extrapolated to a lifetime Hours of operation at dierent wind speeds (binned).Total operation is 102013 hours.

Table 2 :
Time scaling for 30 seconds simulation time (no adjustment).

Table 6
Table 6 but not completely the same.This also