Modeling of Wind Turbine Gearbox Mounting

In this paper three bushing models are evaluated to find a best practice in modeling the mounting of wind turbine gearboxes. Parameter identification on measurements has been used to determine the bushing parameters for dynamic simulation of a gearbox including main shaft. The stiffness of the main components of the gearbox has been calculated. The torsional stiffness of the main shaft, gearbox and the mounting of the gearbox are of same order of magnitude, and eigenfrequency analysis clearly reveals that the stiffness of the gearbox mounting is of importance when modeling full wind turbine drivetrains.


Introduction
Multibody dynamics has been used for load calculation for wind turbines for more than 20 years.Peeters (2006) lists 15 different computer codes for simulating wind turbines.Some of these codes are used for certification of wind turbines, and they have in common that they are intended for fatigue analysis.Most effort has usually gone into modeling the wind loads.The structural part is often limited to between 16 to 24 degrees of freedom (DOF) (Peeters, 2006, p59), and rarely attention has been given to details of the gearbox.The relatively simplified gearbox models are well suited for simulation of large time sequences required for certification (over 10 5 [s]).However, there is an increased demand for more detailed analysis, which includes load distribution in the drive train that is not accommodated by said codes.
An element that is often overlooked or not given detailed attention when modeling geared wind turbine drive trains is the interfaces between the gearbox, main shaft and the nacelle bedplate.In most wind turbines the configuration of this combination is statically indeterminate; hence flexible mounting elements are introduced.Rubber bushings are often used for connect-ing the gearbox to the bed plate through torque arms mounted on the gearbox.
The main torque path in a geared wind turbine drive train enters via the rotor and leaves via the bushings, i.e., the effective stiffness of this torque path may be seen as a number of series connected rotational springs that roughly may be divided into three that represent the main shaft, the gear meshes of the gearbox and the bushing suspension, respectively.
The literature is scarce on rubber mounts for wind turbine applications.Peeters (2006) uses one flexible element to connect the gearbox to the tower top -hence this stiffness must represent the stiffness of the yaw system, bed plate and the rubber bushings combined.
Extensive research has been done on rubber bushings for road vehicle and railroad suspension applications.Bushings have large influence on cornering capabilities and noise, vibrations and harshness (NVH) of road vehicles.The bushing models developed by Wineman et al. (1998), Lee and Kim (2002), Ledesma et al. (1996) and Svensson and Håkansson (2004) have been developed for road vehicle applications.The bushing model presented by Berg (1997) is developed for simulating the performance of train suspensions.The rubber bushings used in road vehicle suspensions are of similar design as those used in wind turbines, only smaller.
A commonly used bushing model is the Kelvin solid Lee and Kim (2002) which is represented by a spring in parallel with a viscous damper, this model is also known as the Voigt model (Lee and Kim, 2002).Ledesma describes this model as state-of-theart in multibody simulation (MBS) of vehicle suspensions (Ledesma et al., 1996).
Other references take a more general perspective on rubber.Pipkin and Rogers (1968) explains how creep force can be modeled for general purposes.Ledesma et al. (1996) has further developed Pipkins creep model into a bushing model for MBS.Ledesma approaches the bushing modeling by using the theory of nonlinear viscoelasticity.This model has been further extended by Wineman et al. (1998) with a view to fit the model to experimental data.
The approach used by Svensson and Håkansson (2004) is to combine a number of different flexible elements.They suggest using a nonlinear spring in parallel with a number of fluid elements and Elasto-plastic elements.Berg (1997) uses a similar approach, only he uses a linear spring, one friction force, and one viscous force.The flexible elements are not formulated in the same way although they are supposed to model the same phenomenon.
A fourth approach is to model the bushing by a transfer function as done by Lee and Kim (2002).Most bushing models are in one dimension, hence, the authors assume the stiffness is statically decoupled Ledesma et al. (1996), Svensson and Håkansson (2004), Wineman et al. (1998).The references that do not explicitly state this assumption are modeling bushings with one DOF and therefore they cannot take into account coupling effects.
The advantage of using measured data is that it describes the bushing in a useful way for MBS.The disadvantage is that it only describes the actual bushing that is measured e.g.experiments have to be repeated for each bushing (Wineman et al., 1998).Most of the models presented in the references are intended for implementation in commercial MBS codes: Lee and Kim (2002), Svensson and Håkansson (2004), Ledesma et al. (1996), Wineman et al. (1998) and Berg (1997).
In this paper the following models are investigated: • Linear spring and damper • Nonlinear spring and linear damper • The Hydro-dynamic bushing presented by Svensson and Håkansson (2004).
These models are incorporated into spatial multibody models and calibrated based on experimental results for a typical drive train with a gearbox having a single planetary stage and two parallel stages.The rotational stiffness of the power conducting part of the drive train is modeled by means of a spatial multibody model of the main shaft and the gearbox including flexibility of shafts, bearings, planet carrier and gear meshes.
The main contributions from this paper are twofold: • to put forward a best practice for modeling of bushings when simulating drive trains.
• to compare the torsional stiffness of the gearboxto-bedplate connection with that of the power conducting drive train.

Considered System
A gearbox with a main shaft is used as experimental and theoretical study in this paper.The gearbox and main shaft are coupled to a dynamometer in an experimental set-up that is shown in Figure 1.The motor of the dynamometer drives the main shaft through the dynamometer shaft and coupling.In the test set-up it is also possible to apply axial and radial loads on the main shaft that are independent of the torque load.In this work this type of loading is, in general, referred to as NTL (non-torque load).The radial loads are applied by two hydraulic cylinders that are positioned in the yz-plane in such a way that a radial force in the y-and z-direction can be applied without any nominal influence on the shaft torque.The gearbox, see Figure 2, is composed of three gearing stages; one planetary stage and two parallel.The four shafts are: main shaft (MS) which is attached to the planet carrier (PLC) using a hydraulic shrink fit, the low speed shaft (LS), the intermediate speed shaft and the high speed shaft (HS).The gearbox has three planets (PL) in the planetary stage.The gear ratio of the gearbox is 81.5.
The output shaft of the gearbox is loaded by a generator as when the gearbox is used in a wind turbine.During tests, the gearbox, main shaft and generator is mounted on a bed plate from a wind turbine.The main bearing is bolted to the bedplate forming a stiff connection compared to the bushings that fastens the gearbox to the bedplate, see Figure 3.
A right hand coordinate system with the x-axis pointing downwind is used, see Figure 1.The origin is at the axis of rotation of the main shaft and located in the center between the bushings.

Bushing Modeling
To investigate different bushing models relative to measurements a simulation model of the entire gearbox The mass of the dynamometer shaft and dynamometer coupling is significant.The dynamometer shaft is supported in the opposite end of the gearbox (indicated with red in Figure 1).To reduce the complexity of the model an equivalent mass of the dynamometer shaft and coupling has been obtained using parameter identification, where measurements of bending moment and rotational position of the main shaft was used to obtain the equivalent mass.By using the equivalent mass it is ensured that the calculated bending moment is equal to the measured.The equivalent mass is 4100[kg], and center of gravity is located 1750[mm] upwind relative to the main bearing, indicated in Figure 2.
The Adams/View model of the gearbox considers only the gearbox and the main shaft, shown in Figure 2. It has been chosen not to model the gears of the gearbox in order to speed up the parameter identification of the different bushing models.In all, the simplified gearbox model consists of three bodies: housing, main shaft and dynamometer coupling.The bed plate is included in the model as frame.The only flexibility that is introduced to the model is that of the bushings.The main shaft is constrained to ground using a rigid revolute joint.The main shaft and housing are coupled using kinematic joints.The high speed shaft is neglected in this model.Due to the high gearing ratio  The bushing models described here consider no coupling effects between deformation axes.Further it is assumed that the bushings have the same properties in the radial y-and z-axes.

Linear bushing
A model with linear damping and linear spring stiffness is used as the first choice.In that case we have the following relationship between force, deflection and deflection rate:

Nonlinear bushing
To account for nonlinear stiffness and deflection dependent damping a nonlinear model is also introduced: The nonlinear bushing model uses four parameters to describe the properties in axial direction and four in radial direction; in all eight constants.

Hydro-Dynamic Bushing
The Hydro-Dynamic Bushing is developed by Svensson and Håkansson (2004).The force may be written as the sum of an elastic force, a fluid force and an elastoplastic force: where the elastic force is given as The fluid force is computed from The states of the intermediate fluid inertia are obtained by time integrating the corresponding equation of motion: where u m is a coordinate giving the position of the mass m.Finally, the elasto-plastic force must be accumulated according to Where the superscript j indicates current time step and j − 1 indicates previous time step.The maximum (yield) force of the elasto-plastic coupling is F y .
The Hydro-Dynamic model uses 8 parameters to express the material properties for both axial and radial displacement; in all 16 constants.case has a duration of 655[s] and is characterized by a number of distinctly different intervals.Firstly, a period of no excitation and no motion is observed.Next, radial NTL is gradually introduced followed by starting the actual rotation of the drive train.Next, the rotational speed is gradually increased and, finally, some torque is added to the main shaft.Quantitatively, this may be expressed as: • From approximately 260[s] to 300[s] the speed drops to zero and rises again.This transient was not planned; however, it is included in the parameter identification.
• At around 535[s] the torque is increased (see Figure 6).
The gearbox proximity sensors are zeroed with the gearbox at rest with no torque, no NTL and no rotation.Hence, any deflection measurement is relative to those caused by the weight of the gearbox, main shaft and dynamometer coupling.The motion sensors are located as shown in Figure 5.A difference between s pz and s sz would indicate a twist in the gearbox about the x axis, θ x .The sensors for detecting motion along the x-and z-axis (s px , s sx , s pz and s sz ) are located close to the bushings at each side of the gearbox.The sensor s bx , s px and s sx measure the rotation about the y-axis while the motion along the y-axis is measured only by s py at the port side of the gearbox.The measured data from the motion sensors are shown in Figure 6.The six measured positions are transformable to rigid body motion of the gearbox expressed in the coordinate system shown in Figure 5.The transformation can be expressed as where The displacements corresponding to the measured data may then be computed as: These values that are measured in the shown coordinate system may be compared with those obtained for the gearbox housing from the computer model.The equations (8)... ( 10) are only valid for small rotations.The derived gearbox rigid body motion is plotted in Figure 7.

Parameter Identification
The strategy of the parameter identification is to define an objective function that reflects deviation between measured and simulated results and minimize that by means of an optimization scheme that uses the parameters to be identified as variables.This corresponds to design optimization where the design is the model parameters and the design performance is the correlation between simulated and measured gearbox motion.
The optimization is performed using the complex method.Its main advantage, in this context, is that it does not require any gradient that otherwise would have had to be determined numerically because of the interaction between Adams/View and Matlab.Potentially, any non-gradient method may be time consuming because of a high number of iterations.This has, however, not been an issue in the current work and therefore alternative methods have not been investigated.The complex method generates a population of designs randomly distributed within some specified limits.A design is a vector consisting of the design variables of a bushing model.Hereafter the worst design is mirrored in the mean of the rest of the designs.This is repeated until the difference between the worst and best design has reached a specified tolerance (Box, 1965).The population size used here is, in general, two times the number of model parameters.These relatively small populations has yielded both consistent and satisfactory results with small computational costs.
The parameter identification is illustrated by means of a flow chart diagrams in Figure 8.The complex optimization is programmed in Matlab whereas the time domain simulation for the evaluation is carried out in Adams/View.The Adams simulation is started with the values from the current design that the complex algorithm needs to evaluate.In the Adams simulation the measured torque of the main shaft and the measured NTL are used as input to the simulation model.The output from the simulation is the displacement of the housing.The only parameter that is changed from one simulation to the next is the design variables that describe the bushing properties.
The output from the simulation is given as timeseries collected into one matrix: 11) where n is the number of time samples.In the pre-processing of the parameter identification ( 8) is used to convert the measured data to rigid body motion of the housing that is comparable to those obtained in Adams.The objective function value is Where ∆d i = d m i − d a i and W is a diagonal weighing matrix.In order to ensure that the sample time of the measured data matches the simulated, the measured data are interpolated.The dominant weighing term has been the one associated with θ x since this is in series with the main degree of freedom of the drive train.
The sensors on the gearbox are zeroed at standstill, which means the gearbox is subjected to gravity when the sensors are zeroed.This is accounted for by adding the initial offset from the measurements, d m ini , obtained from the first sample, (13).The virtual sensors are zeroed without gravity applied in the Adams model; here gravity applied is when the simulation is started.The effect of gravity is calculated as an average of the first 20 samples of the simulation as given by ( 14).Since 520[s] is simulated, the first and last part of the measurements are omitted, and 2000 samples are requested from a simulation the first 20 samples corresponds to 5.2 [s].
The parameter identification is conducted on three Adams models; one for each of the bushing models.In addition to the bushing parameters a parasitic torque that originates from the NTL is introduced as a parameter.The parasitic torque is added to the main shaft torque and is calculated in this way: Where F N T L,z is the vertical load produced by F N T L,p and F N T L,s .This is based on reports from observations drawn from the experiments that the radial NTL were not applied ideal, i.e., without introducing an additional torque.The parameter identification has justified this observation.

Results
The bushing models are compared on a number of parameters.The first parameter is relatively straight forward, since the complexity may be represented by the number of independent parameters that must be determined.Therefore it is desirable to have as few constants as necessary.The computational time should be  According to the results presented in Table 2 the linear bushing is best in terms of number of parameters and computational time whereas the nonlinear bushing is better on accuracy.Obviously, decisive conclusions on accuracy can only be drawn based on several measurements, however, within the scope of this work it is safe to conclude that the relatively simple linear and non-linear models seem a better choice than the more complex hydro-dynamic model overall wind turbine drive train modeling due to their ability to capture the important dynamic characteristics with relatively few parameters.
The parameters obtained by parameter identification are listed in Table 3. Results of the parameter identifications are shown in Figure 9 for θ x .Since the model produces six sensor signals that are comparable to the measured data presented in Figure 7 these are plotted together in Figure 10.In general, it seems that the ability to capture the torsional degree of freedom is quite good.There are obviously some offset errors on, particularly, the z-deflection and the y-rotation.

Torsional Stiffness
The gearbox has three stages that all employ helical gears.The gear modeling tool used features a 3D contact algorithm and is capable of modeling helical gears.
The basic gear data is assembled in Table 4.
The gearbox has a planet carrier which is mounted to the main shaft using a press fitting, therefore, in the flexible model the planet carrier and the main shaft have been joined using bonded contact.All shafts have been meshed using tetrahedron elements; the mesh is relatively coarse because only stiffness is of interest.The bearings are modeled using linear spring and  The housing is considered rigid; its movement is constrained to ground preventing it from moving.The spherical main bearing is modeled by a kinematic constraint only allowing rotation about all 3 axes.The output shaft is locked against rotation about the xaxis.A torque is gradually applied at the hub flange of the main shaft.During simulation the rotation of the PLC was recorded at the upwind bearing of the PLC.A force-rotation plot can be seen in Figure 11.The stiffness is obtained by linear regression.
The stiffness of the gearbox without bushings and main shaft is 27[MNm/rad] while the torsional stiffness of the main shaft is 50[MNm/rad].When comparing with the torsional stiffness of the bushings of 11.4[MNm/rad] which means all of the torsional stiffnesses are in the same order of magnitude.

Eigenfrequency
The main implication of the torsional stiffnesses identified in this paper is best expressed by means of the lowest eigenfrequency of the drive train.This has been carried out on a simplified and linearized model in the commercial simulation software SimulationX.In Figure 12 the main inertias, torsional stiffnesses and kinematic constraints are shown.
The stiffnesses of the meshes have been transferred to the MS-PC axis via the simulation results of section 5.2 and the kinematic constraints may be derived from the gear data in Table 4.
Two different first eigenfrequencies are found depending on whether the torsional stiffness of the bush- .11[rad/s],respectively.Clearly, there is a substantial relative difference.The importance of this difference in dynamic characteristics will affect both the transfer of power in the drive train, in general, and the transmission of torque from the gearbox to the bedplate in particular.

Conclusions
One of the main purposes of this work has been to investigate a best practice for modeling gearbox bushings used in wind turbine drive trains.Comparison between three different models, a linear, a non-linear and a hydro-dynamic bushing has been performed by means of parameter identification from experimental results.All three models are capable of displaying the main dynamic characteristics of the gearbox bushing in an adequate way, and the relatively simple linear and non-linear models are easily competitive with the more complex hydro-dynamic model on accuracy.
A second purpose was to estimate the importance of including the bushing model in an overall drive train model.This has led to the investigation of the three main torsional flexibilities within the drive train: main shaft, gearbox and bushings.Clearly, the torsional stiffness of the bushings should be included.In the gearbox examined throughout this paper it was shown that the lowest eigenfrequency of the drive train would increase with approximately 50% if the bushing flexibility was omitted.

Figure 2 :
Figure 1: The dynamometer set-up Hub flange for connection to dynamometer coupling Main shaft (MS) Planet carrier (PLC) Sun shaft Low speed shaft (LS) IMS

Figure 4 :
Figure 4: Graphical representation of hydro dynamic model

Figure 5 :
Figure 5: Locations of motion sensors

Figure 6 :
Figure 6: Measured data from the investigated load case

Figure 7 :
Figure 7: Data used in the parameter identification

Figure 8 :
Figure 8: Flow charts of the parameter identification

Figure 10 :
Figure 10: Comparison of results obtained with nonlinear bushings

Figure 11 :
Figure 11: Torque as function of main shaft rotation

Figure 12 :
Figure 12: Diagram of SimulationX model of the gearbox.Inertia, stiffnesses and gear ratios are added

Table 1 :
Bushing displacements have been measured with the dynamometer subjected to a number of different excitations that have been assembled in a single load case.The measured variables are listed in Table1.The load Measured signals, see also Figure5

Table 2 :
Benchmarking table The computational time is the time it takes to simulate the full time-series(520[s]).On the other hand, the model should produce accurate results.The accuracy is measured by the objective function value calculated using the constants found by the parameter identification.If the gearbox model is able to represent the real gearbox and the bushing model is suitable, then the objective function value should approach zero.The values presented in Table2are normalized with respect to the accuracy of the hydro-dynamic bushing.

Table 3 :
Parameters obtained using parameter identification.All units are SI-units

Table 4 :
Basic gearbox data