Modeling and Parameter Identification of Deflections in Planetary Stage of Wind Turbine Gearbox

The main focus of this paper is the experimental and numerical investigation of a 750[kW] wind turbine gearbox. A detailed model of the gearbox with main shaft has been created using MSC.Adams. Special focus has been put on modeling the planet carrier (PLC) in the gearbox. For this purpose experimental data from a drive train test set up has been analyzed using parameter identification to quantify misalignments. Based on the measurements a combination of main shaft misalignment and planet carrier deflection has been identified. A purely numerical model has been developed and it shows good accordance with the experimental data.


Introduction
A common topology of wind turbine drive trains includes a gearbox with the main purpose of increasing the speed of the rotor to a value suited to the generator, power electronics and grid at hand.The steady increase in power rating of wind turbines together with the constraint on the blade tip speed has led to higher demands for gear ratio and, subsequently, to a higher complexity in the gearbox layout and in the load distribution.This trend has further emphasized the gearbox as a crucial element in any reliability assessment of a wind turbine.
Since model based prediction of performance, in general, and reliability, in particular, holds many advantages the obvious conclusion from the above is to employ models that take into account more degrees of freedom and more phenomena when doing time domain simulation of drive trains.This should, however, be done carefully since certification in the wind turbine industry is based on simulation of large time sequences.This work focuses on the investigation of the importance of including gearbox component flexibility with special emphasis on the planet carrier (PLC) of the planetary stage normally found in a wind turbine gearbox.In that respect the work can be viewed as a continuation of Haastrup et al. (2011) that focused on the distribution of the torsional flexibility of the drive train.Also, the work can be viewed as an integral part of the Gearbox Reliability Collaborative (GRC) initiated and managed by the National Renewable Energy Laboratory (NREL) where a 775 [kW] gearbox has been subjected to extensive testing.Part of the instrumentation involves proximity sensors used to measure the axial and radial motion and deflection of specific points of the PLC.The aim of this paper is to utilize these measurements to set up and calibrate a model that properly describes the behavior of the PLC with a minimum of complexity and computational costs.
Many different approaches have been suggested on how to model planetary gears.In the following papers researchers are investigating vibration phenomena using analytical approaches: Al-shyyab and Kahraman (2007), Kahraman (1994), Lin and Parker (1999).Others reach the same objective using finite element Parker et al. (2000a), Parker et al. (2000b) and Rigaud et al. (2006).Common for all of these papers is that they present methods for prediction of noise and/or vibration of gearboxes.Before computers became sufficiently fast to solve large finite element analyses (FEA) and multi body simulations (MBS) the same objectives were reached using analytical models.A survey is given in Özgüven and Houser (1988).
To model wind turbine gearboxes a gear model is required and some gear models for use in MBS has been suggested.Lundvall et al. (2004) proposes a gear model that imposes flexibility to the rigid body motion.The flexibility is obtained using FEM and tribology theory.
Another method proposed by Blankenship and Singh (1995) implements full 3D forces and moments relevant for a helical gear pair.Blankenship uses a combined stiffness to account for both gear wheel, tooth and contact stiffness.A third gear model suggested by Ebrahimi and Eberhard (2006) uses a lumped mass approach for modeling the flexibility of gears.In the present paper commercial available modeling tools developed by MSC.Software TM , see Table 1, is used in conjunction with measured data to set up guidelines for modeling planetary gearboxes in wind turbine drive trains.The gear model has a contact algorithm for detecting contact between mating gears taking profile modifications and contact stiffness into account.
In the present paper a gear model developed by MSC.Software in Germany will be used.The model has a contact algorithm for detecting contact between mating gears taking profile modifications and contact stiffness into account.Ultimately the goal is to predict fatigue loads on gear, bearings and other components and noise emission.However, none of these models are experimentally verified.

Identification of PLC motion
The gearbox considered in this paper has three gear stages: one planetary and two parallel.The naming of the main components of the gearbox is shown in Figure 1.The PLC motion is measured by six proximity sensors.Four of them are measuring axial displacements (along the x-axis) and two are measuring radial displacement (perpendicular to the x-axis).All the displacements are measured relative to the housing of the gearbox, i.e., any displacement is a combination of rigid body motion of the PLC relative to the housing as well as deflection of both the PLC and the housing.
In order to produce valid experimental data, two faces on the PLC have been machined, as shown in Figure 2, and they are used as measuring reference.All sensors are mounted on the inside of the main flange of the gearbox.They are displaced a radius R s = 412[mm] from the x-axis, see Figure 2, and, further, rotated different angles around the x-axis with 0 • corresponding to the z-axis.The four sensors that measure axial displacement are located at θ x = [47 • , 137 • , 227 • , 317 • ] and the two sensors that measure radial displacement are located at θ x = [40  Identification of the PLC motion requires a model that is capable of showing the same behavior as can be observed in the measurements.From the observations two models are proposed: The first model concerns only rigid body motion (M1) while the second model also takes into account the deformation of the PLC (M2).
Two coordinate systems have been employed to describe the assumed rigid body motion of the PLC.Firstly, there is a fixed coordinate system (x H , y H , z H ) that is shown in Figure 3 and which is aligned with the global coordinate system (x, y, z) shown in Figure 2. Secondly, a coordinate system fixed to the PLC (x P LC , y P LC , z P LC ) is defined.This coordinate system is rotated around the x H -axis of the fixed coordinate system with an angle θ x (t), see Figure 3.The rigid body model assumes that the motion of the PLC can be described by a constant displacement, d x , along the x H -axis, along with a constant tilt φ H around a fixed tilt axis that is obtained by rotating the z H -axis an angle θ H around the x H -axis. Similarly, a constant tilt φ P LC around a rotating tilt axis that is obtained by rotating the z P LC -axis an angle θ P LC around the x P LC = x H -axis.In all this gives five parameters to be identified: d x , θ H , φ H , θ P LC , and φ P LC that describe the rigid body motion of the PLC.In the following equations θ x (t) denotes the rotation of the PLC defined as the rotation of the PLC-coordinate system relative to the H-coordinate system.The simulated sensor signals: s 47 m , s 137 m , s 227 m and s 317 m are obtained from the model.They are the motions in the x H -direction and may be computed for small tilt angles as: Where R s = 412[mm] is the radius from the x-axis to the sensor locations, and is the angle from the z H -axis to the sensor, v, location.The angle from the z P LC -axis to sensor v is: The relative angles defined in eq.( 2) and eq.( 3) are illustrated for a single sensor in Figure 3. Referring to eq.( 1) the location on the PLC having zero axial displacement due to misalignment of the housing relative to the main shaft is [R s , θ H ] in polar coordinates; the blue line in Figure 3 shows this location.The largest displacement caused by housing misalignment is given by R s • φ H .To obtain the displacement at a given sensor this has to be multiplied by sine to the angle from θ H to the sensor v, calculated in eq.( 2).This displacement is constant in the housing coordinate system and thereby at the sensors.
Similarly, the PLC is mounted with some misalignment relative to the axis of rotation.Again, the point where the misalignment yields zero displacement can be expressed in polar coordinates: [R S , θ P LC ], shown by the orange line in Figure 3. θ P LC is constant in the coordinate system rotating along with the PLC, therefore the rotation of the MS needs to be accounted for when calculating the angle from θ P LC to sensor v; this angle is given by eq.(3).In Figure 4 the different meanings of the two tilt angles are shown exaggerated.

Deformation model
For the purpose of determining the deformation pattern of the PLC a FE model is created and a unit  drive torque of T = 1[kNm] is applied.The Planets are mounted in the PLC using planet pins which are also modeled.See Figure 6 for the planet pins.The tangential force caused by the drive torque is applied at each planet pin in the FE model.Since helical gears are used in the gearbox, a tilt moment is applied on each planet.The definition of helix angle and pressure angle is explained in Figure 5b.The force normal to the tooth and the axial force is respectively: and the tilting moment is

Planet pin
Figure 6: Boundary conditions applied to the FEM.
where d = 400.4[mm] is the pitch diameter of the planet gears.The basic gear data for the planetary gear is given in Table 3.The upwind end of the PLC, where it is connected to the main shaft, is constrained to zero displacement and the force F P LC is applied at each planet pin, see Figure 6.Also, the tilt moment M t is applied.The forces and moments are applied through massless rigid body elements.
The planet pins are mounted in the PLC using a clearance fit, which preferably should be modeled using contact elements.However, detailed knowledge of the actual tolerances are not available and, instead, it is chosen to model two extremes: a stiff model where the pins are structurally connected to the PLC and a flexible model where the pins are absent from the structure and simply transmit forces to the PLC as if it was simply supported.It is believed that a mix of the corresponding deflection patterns would be capable of representing the actual deflections.The deformation patterns are presented in Figure 7 and 8 where the axial deflection of the rim is plotted.
To investigate the influence of the tangential load relative to the tilt moment, analyses have been con- ducted with the loads applied separately and together; the conclusion is that the tilt moment is insignificant.The deflection of the PLC is modeled by scaling the curves of Figure 8.A model of the deformation is proposed: where The parameters ρ s and ρ f define the scaling of the mode shapes shown in Figure 8.A model of the PLC deflection has been created where the data labeled δ s (θ) in Figure 8 scaled by ρ s and the data denoted δ f (θ) scaled by ρ f eq.( 7) is used.
In the following two different models are investigated.Model (M1) simply employs the rigid body motion as described in eq.( 1) and therefore has the following five parameters to determine: d x , θ P LC , φ P LC , θ H , and φ H . Model 2 (M2) includes the flexibility of the PLC using eq.( 5) and therefore has two additional parameters to be determined: ρ s and ρ f .

Parameter identification
There is no exact solution for finding the five parameters used in eq.( 1) therefore they are determined by optimization.The result of the parameter identification is shown in Figure 9 for the test, LC2, where 50% torque is applied to the gearbox.For each load case the RMS values of the difference between the measured and simulated displacements has been calculated and is shown for both M1 and M2 in Table 4.In Figure 10 and 11 the parameters are plotted.There is a sub plot for each of the parameters in the rigid body model.For M2 the deflection parameters are shown, see Figure 11.
When the PLC is moderately loaded, less than 50% of nominal torque, the contact between the PLC and pins does only transfer force and no bending moment, but under larger loads the pins are stiffening the PLC and the deformation pattern changes to have six peaks instead of three.This behavior can be verified by looking at the value of ρ 2 which is almost zero in the load case with 50% torque (LC2).In the load cases with higher loads ρ 2 has a value different from zero which means that more peaks are present.In general, it is concluded that the proposed models describe the motion of the PLC quite well.It is also obvious, that the overall motion is a combination of rigid body motion and a superimposed deflection.Clearly, the shape of the M2 model shows a better resemblance with the measured data than that of the M1 model, however, this is only partly reflected by the actual RMSdeviations as shown in Table 4.

Numerical simulation model
In this section a time domain simulation model of the gearbox is introduced.The main purpose of this model   is to compare the measurements with the results of the parameter identification and, at the same time, evaluate the influence from the modeling on the estimated life of the planet gears.Since emphasis is on gear box modeling the misalignment between the main shaft and PLC identified in the previous section is imposed on this model.The modeling is done in the commercial software Adams/View and the description of the model is divided into four sub sections: the bodies, the gear meshes, the bearings and the boundary conditions.

Bodies
All of the shafts in the gearbox can be considered as flexible bodies.However, to reduce the computational time the LS, IMS and HS, see Figure 1, are considered to be rigid.Previous studies have shown that introduction of flexible shafts only have local impact on motion and load distribution.Therefore, the exclusion of the flexibility of these shafts will not affect the motion of the PLC (Haastrup et al., 2010).Similarly, the housing can also be considered flexible and since the sensors are mounted on the housing its deflection will influence the measurements.However, initial studies have shown that the deformation of the housing is one order of magnitude smaller than that of the PLC and therefore the housing deflections are neglected.Flexible bodies are modeled using an assumed modes method that is available in MSC.Adams.This method relies on the assumption that the deformation of a body under load can be described by a finite number of mode shapes.The mode shapes are calculated using finite element software, in this case MSC.Nastran.The mode shapes are composed of a selected number of free body modes and fixed boundary modes (one static mode for each DOF).All of these modes are combined using Component Mode Synthesis, the Craig-Bampton method (MSC.Software, 2010).Each fixed boundary mode is associated with one boundary DOF.When a flexible body is imported into Adams/View the boundary DOF's become interface nodes, which can be used for applying forces and constraints.Before the MS and PLC were meshed using MSC.Patran, all small details such as fillets, chamfers and grooves were removed to allow the mesher to create a coarser mesh.It is assumed that the stiffness is not significantly changed by removing these details.
The MS and PLC are modeled as separate flexible bodies that are joined together using a kinematic constraint to allow for the modeling of the misalignment described in section 2.1.The planet pins are clearance fitted to the PLC and this can be modeled in several ways.In this study the interface between the planet pins and the PLC has been modeled using massless rigid elements where one node (pilot node) defines the DOFs for the other nodes, see Figure 13.The effect of using the pilot node approach is shown in Figure 12, it can be seen that the pins are allowed to pivot relative to the PLC at the interfaces.
In Figure 13 the pilot nodes on the PLC are shown.There are three planet pins, one for each planet.Each planet pin is connected to the PLC in both ends, which gives six interfaces.An interface consists of two pi- lot nodes, one connected to the planet pin and the other to the PLC.The pilot nodes are connected to the structure (PLC or pin) using so-called RBE3 elements, which form a rigid connection to the structure.The RBE3 elements use the surfaces of the interface where the clearance fit is.The two pilot nodes of an interface are kinematically connected locking all translational DOF's.Some bending stiffness has been added to the interface using torsional springs with a stiffness of 10 6 [Nm/rad] to model the six-peak behavior of the stiff model, see Figure 7 and 8.The torsional stiffness of the interfaces is selected to be equal to the bending stiffness.

Gear model
The main requirement to the gear model is that it is able to model the tangential and radial forces acting between the gear pairs.For this purpose Adams Gear Generator is selected, which calculates a contact force based on involute curves.Adams Gear Generator en-ables modeling of spur gears, hence the effects of the helix angle, β = 7.5 • , must be neglected.The axial forces caused by the helix angle are opposite for the planet-sun and planet-ring meshes, thereby creating a moment on the planet.An FE model has shown that this moment has negligible influence on the axial deformation of the PLC.
Other possibilities exist to model gears, for example Adams/Gear AT which can model helical gears.It uses a 3D contact algorithm to calculate contact forces but it is several order of magnitude more time consuming than the Adams Gear Generator model.Since it does not provide any significant extra information regarding the PLC behavior it is discarded.

Bearing model
Bearings can be modeled to various levels of details.The simplest way is to use kinematic joints; however this method does not take the bearing flexibility into account.Models exist that take into account the contact stiffness of each rolling element in the bearing.These models suffer from being less computationally efficient as kinematic joints and they require detailed information about the geometry of the bearing which is often difficult to obtain.The level in between is to use a three or six component linear stiffness element.It is only possible to use this model if the bearing stiffness can be obtained from some other analysis tool or by the bearing manufacturer.Since the fatigue life of the planet bearings is to be evaluated a detailed bearing model (Kabus et al., 2011) is used at the planets while the rest of the bearings are modeled using linear springs and dampers.
In Adams/View the linear stiffness model has three linear stiffness components: k x , k y and k z and three rotational linear stiffness components: k yz , k xz and k xy .Since the axis of rotation is the z-axis the stiffness k xy has to be zero, all other components can have a value.For each bearing in the gearbox linear stiffness values have been provided by GRC, which have been linearized for loadings corresponding to nominal working load.
To simplify the Adams/View model some of the bearings have been merged into one bushing.For example, each of the shafts supporting the HS and IMS are supported by a pair of tapered roller bearings which are modeled by one linear stiffness element per set.

Boundary conditions
The bedplate, which the gearbox is mounted to, acts as the ground in this model and is therefore modeled as rigid.In order to model the rigid body motion of the MS relative to the housing as investigated in the parameter identification it would be necessary to include the flexibilities of both the bedplate and the main bearing.That is, however, outside the scope of this paper and the main bearing is constrained to ground using a spherical joint which removes all translational DOF's.
Two rubber bushings are used to mount the gearbox to the bedplate; these are modeled using a three component linear spring and damper model.The stiffness and damping values are calculated by Haastrup et al. (2011).
In the tests a drive torque was applied to the main shaft by a dynamometer.This drive torque was mea-sured and used as input to the model.In the model this torque is applied at the flange of the hub, which is also the case for the dynamometer testing.The dynamometer is also capable of applying radial and axial loads, see Table 2 for overview of load cases.In the experiment and in wind turbines in general, the rotational speed of the shafts is determined by the generator controller and/or the torque on the main shaft.The main shaft torque is controlled by the dynamometer or pitch system.In the model the rotational speed of the shaft is controlled on the output shaft.The torque applied on the high speed shaft is where ω s (t) is the desired rotational speed which has been measured in the experiments, ω is the actual simulated rotational speed and K = 10 3 [Nms].Eight simulations have been conducted which are summarized in Table 5.In this table the RMS values showing the simulations deviation from the measurements are also listed.The misalignment between the main shaft and the PLC has been implemented by introducing the θ and φ angles obtained in the parameter identification.
From Figure 14 it is seen that the same tendencies as observed in the parameter identification are repeated, namely that the rigid body motion of the main shaft relative to the gear housing is dominant and that the deflection of the PLC can be seen superimposed on this.The RMS deviations clearly show that the inclusion of the main shaft misalignment is most important.The inclusion of the flexibility of the PLC has the same effect as seen earlier in the paper, namely an improved resemblance between modeled and measured deflections but with limited impact on the RMS deviations.In general, the proposed gearbox model is capable of reflecting the actual behavior of the PLC; however, the MS-PLC misalignment must be forced upon the model from the measurements.In order to get this misalignment from a model it would be necessary to include more components from the drive train, such as the main shaft, the bed plate and the main bearing.

Gear life reduction due to misalignment
By introducing a more advanced gear model (Adams/Gear AT) into the Adams/View model described in section 3 the contact stress can be obtained for the gear mesh.The contact algorithm within Adams/Gear AT uses a fine surface mesh (MSC.Germany, 2011).Adams/Gear AT is able to change the contact pattern when the gears are misaligned.In the simulations conducted here crowning and tip relief were modeled.Simulations have been conducted both with and without MS-PLC misalignment in order to evaluate the influence.For the planetsun an increase in contact pressure of 1.5% is found.The contact pressure is considerably higher for the planet-sun mesh than for the planet-ring mesh due to different geometry.An increase in contact pressure will reduce the surface lifetime.The magnitude of lifetime reduction may, for example be estimated based on the AGMA Surface-Fatigue strength life factor C L (AGMA, 1989).Based on the relation between C L and number of load cycles N, the following relation can be derived where S f c,1 and S f c,2 are the fatigue strengths equivalent to N 1 and N 2 load cycles respectively.By inserting the change in contact pressure the reduction in lifetime can be expressed as It is noted that an increase in contact stress by 1.5% decreases the expected life time by 23%.This result indicates that misalignment of gears has a significant influence on expected life time and thus needs to be taken into account.

Conclusions
In the presented work best practice for modeling of gearboxes in wind turbine drive trains has been investigated.Emphasis has been on the planetary stage and the axial displacement of the planet carrier relative to the gear housing.Experimental work combined with numerical studies have shown that the misalignment between the main shaft and the planet carrier is the dominant source for these displacements and, subsequently, on the expected life of the planet gears.Adding flexibility of the planet carrier clearly improves the resemblance between the measured and modeled deflections, however, this is not reflected in the RMS values of the deviations between measured and modeled.The experimental and numerical results presented in the paper clearly indicate that although the flexibility of the planet carrier influences the load distribution within the gear box, it is more important to have a model that reflects the planet carrier misalignment.

Figure 1 :
Figure 1: Naming of the main components.

Figure 2 :
Figure 2: CAD drawing showing the location of the proximity sensors.

Figure 3 :
Figure 3: Definition coordinate systems and angles.

Figure 4 :
Figure 4: Illustration of the PLC and housing misalignment.F t

Figure 7 :Figure 8 :
Figure 7: Axial deformation as function of y P LC -z P LC coordinates.The deflection pattern of the stiff model is δ s and the deflection pattern of the flexible model is δ f .

Figure 9 :Figure 10 :
Figure 9: Measurements plotted together with results from the parameter identification.

Figure 12 :
Figure 12: Seventh mode shape of the combined MS and PLC.The planet pins are able to pivot.

Figure 14 :
Figure 14: Comparison between measured and simulated deflection of the PLC (LC2).

Table 1 :
Software used in this paper

Table 2 :
Haastrup et al. (2011)In two of the load cases radial load was also applied by means of a non-torque radial loading device, seeHaastrup et al. (2011).The load cases are listed in Table2.Load cases.

Table 3 :
Basic gear data.Length in [mm] and angle in [deg].