Investigation of the load distribution in planetary gearboxes of wind turbines using multibody-system simulation

The high costs for the development, erecting und operation of wind turbines are connected to very high expectations for a reliable and low-maintenance operation and require a precise knowledge of the loads and stresses to be expected. The transfer of knowledge from smaller wind turbines and possibly other concepts succeeds only to a limited extent. Rather comprehensive simulation approaches to determine wind loads, operational conditions and possible resonances are already used since many years. By means of simulation models, the natural frequencies can be determined and compared to possible excitations. The simulation of the operation of the wind turbine under different wind speeds allows the calculation of component loads as a basis for the further design process. The paper concentrates on the possibility of using multibody-system simulation models in the design process of gearboxes for wind turbines and the associated dynamic properties of the complete system by the example of the 15 MW reference wind turbine of the National Renewable Energy Laboratory (NREL). The comprehensive factors which influence the load distribution in the gearing of the first planetary gear stage require a detailed consideration of the elasticity of all relevant components. Based on the developed gearbox design and a detailed multibody-system simulation model, the influence of the level of detail of the model on the resulting natural frequencies and the occurring load distribution in the gearing of the planetary gear stages can be discussed. The present results show that findings on the required level of detail of simulation models cannot be applied to new turbines independently of the power class.


Introduction
Wind turbines are operating under complex operational conditions that are varying daily but also over the year, depending on the operation site and the position in a wind park. In comparison to the dimensioning of large drive systems for the energy production, which are operating under nearly constant operational conditions, the torque component is only one part of the loads that have be considered. Already under normal operational conditions, the shear of the wind field leads to changing thrust loads over the height of the rotor surface. In combination with the inclination angle of the drive train, the resulting load distribution over the rotor leads to a relevant bending moment around the horizontal axis perpendicular to the rotation axis. The pitching of the rotor blades out of the wind reduces the thrust to zero, so that the dead weight of the rotor has be supported by the drive train. The loads that occur are acting in the opposite direction to the described operating load. The analysis of the load situation around the vertical axis can identify the inflow angle as main influence factor. The wind direction is influenced by short-and long-term weather conditions, the location in the wind park and the control strategy of the wind turbine. The azimuth drives react in different sensibilities to a change of the inflow angle. This becomes obvious in wind parks, where in most cases even wind turbines of the same manufacturer are oriented in different directions.
The combination of the described loads and the further consideration of extreme load cases requires a detailed knowledge of the loads to design a sufficient load bearing supporting structure. Dependent on the design of the wind turbine, different concepts are in use to support the rotor loads. Due to limitations on space and weight, not all possible design measures can be realized. The resulting deformation of the drive train has an influence on the alignment between rotor and stator for direct driven wind turbines or on the introduced loads for wind turbines with gearboxes. The following analysis will focus on gear driven wind turbines where especially for the dimensioning of the first gear stage the accurate knowledge of the loads is very important.

Development of detailed wind turbine simulation models
The 5 MW baseline model of the National Renewable Energy Laboratory (NREL) already served as the basis for the first detailed multibody-system simulation model, which has been developed at the Chair of Machine Elements using the software SIMPACK since 2015. For a rotor diameter of 126 meters, a wind speed range from 3 to 25 m/s and a nominal rotor speed of about 12 rpm, different sections of the drive train and the supporting structure were developed. The distribution of the transmission ratio of the individual stages of the planetary-planetary-helical gearbox was optimized using an algorithm in order to achieve a high power density [10]. The constructive design of the gearbox housing, the planetary carriers, the machine carrier, the hub and the nacelle completed the entire model and provided the basis for investigations on the influence of the resilience on the overall system behavior [6,10].
With the publication of the new NREL baseline model for a power of 15 MW, the template for the development of a multibody-system simulation model for a new power class was available [1]. The turbine has a rotor with a diameter of 240 meters and is also designed for the wind speed range between 3 and 25 m/s. The nominal speed of the rotor for this turbine is around 7.5 rpm and the turbine is designed as a direct drive. Extensive information about the baseline model is available, which can be used to adapt the present 5 MW multibody-system model and provide the basis for the development of a new gearbox for this power class.
The rotor blades of the new simulation model are generated using the data set for the 15 MW baseline wind turbine for the two available SIMPACK-approaches "Basic" and "Advanced". In addition to the information on mass and stiffness, the geometries of the blade cross sections over the rotor blade length are required to generate the elastic beam model. The data sets for the rotor blade cross sections are taken from the database and transferred into an appropriate format. Compared to the effort required to create the rotor blades, the model of the wind turbine tower can be created based on the available cross-sectional areas for different tower heights. K Fig. 1 Changing load distribution over the planet carrier rotation An interface to the wind load simulation program Aero-Dyn in SIMPACK can be used to represent the wind loads and requires the adaptation of the corresponding input files. The input files for AeroDyn and for the rotor blade were rebuilt based on the information provided for the baseline wind turbine using AeroDyn version 15. Since the information on 50 sections along the length of the rotor blade is available for the new turbine, the information required for AeroDyn is also provided in a higher resolution than for the 5 MW wind turbine. The information on the structural and aerodynamic properties are transferred from OpenFast to SIMPACK. Due to a different number of nodes for the sections in both databases, it is necessary that data points have to be interpolated during the transfer.

Requirements for the design of the gearing
The rotor-side introduced torque loads must be transferred by the gearing of the planetary gear stage. Compared to the dimensioning of shear torque-loaded gearboxes, the definition of gearing modification must consider next to the torque-dependent misalignment of the gears, for example by the twist of the planet carrier, the misalignment of the planet carrier due to external loads. Such a misalignment leads to an inclination of the carrier and the planets against the ring gear and the sun gear. Under a constant inclination of the planet bolts, the load in the contact between the gears increases on the end, where the axis distance is minimal. That leads to a maximum load in the contact between planet-ring and planet-sun on the opposite ends of the gearing. The rotation of the planet carrier with the planets to a position after 180 degrees shifts the position of the maximum load to the opposite end. This change of the load distribution leads to the challenge, to define a modification of the tooth flank, which achieves an optimal load distribution factor over the complete carrier rotation, Fig. 1, [13].
In addition to considering a change in load distribution over the rotation of the planet carrier, the effects of different load cases and the frequency of occurrence of these load cases must be taken into account. Commercial purchasable load distribution calculation programs can be used to determine a required modification to the gearing. However, the significant influence of the flexibility of the surrounding structure can only be taken into account indirectly by applying a calculated misalignment of the gears to each other [8,11,12]. The determination of the misalignment requires the setup and calculation of different load cases and engagement positions by means of a finite element model. Due to the interactions between the rotor, main shaft, mainframe, gearbox housing, torque supports, planet carrier, and bearing stiffnesses, the required models become very complex, and the calculation times are long. The complete search for an optimal gear modification with a finite element model increases the effort in all previously mentioned points. Using the multibody-system simulation, and taking into account all relevant elasticities and the available force elements,  enables an evaluation of the load distribution in the tooth meshes with a modeling effort comparable to the pure FE approach but with shorter calculation times [14].

Design of the gearbox
The first assembly of multibody-system models of drive trains of wind turbines was done in the first years of the 21st century [4,5,7,9]. The method as well as the used computers were further developed since this time. The main challenges to come to a detailed description of the load distribution for the gearing are the missing possibility to model detailed rotor-side loads and the elasticity of the gearing over many years. Since both open points are implemented in the multibody-system simulation software SIMPACK, the drive train for the NREL baseline models get assembled and further developed [6]. On basis of load assumptions taken out of the NREL baseline model, the design of the gear stages for the planetary-helical gearbox occurs, the calculation and the approval of the bearing lifetime are done, and the construction of the planet carriers takes place.
The design of the gearings occurs on basis of an overview to all constructive possible tooth combinations for all stages and the determination of the resulting gear ratios. Instead of a series connection of individual gear stages, whereby each stage transmits the full input power, a power split in the first planetary gear stage allows the transfer of the input power to the ring gear of the middle stage and the planet carrier of the third planetary gear stage. The second planetary gear stage is designed as a stationary gear stage, the planet carrier cannot rotate. In the third planetary gear stage, the previously divided power flows are recombined by the ring gear and the planet carrier. The helical gear stage ensures the axle offset and the achievement of the required total transmission ratio. To ensure that a gearbox can nevertheless be designed with the lowest possible weight, in a first step possible combinations of numbers of tooth are selected that lead to the required total gear ratio. Based on a rough determination of the tooth root strength and tooth flank safety, all required gear parameters can be calculated. In the diagram (Fig. 2), the standardized masses determined for a gear stage are plotted against the stationary gear ratio for various tooth combinations. Due to the different combinations of teeth that can result for the same gear ratio, the variations in the resulting mass are large.
The available information can be used to derive correlations between the stationary gear ratios and the mass of the assembly before the actual start of the design process. This leads to the conclusion that the stationary gear ratios of the first two stages should be minimized when defining the gear parameters (Fig. 2). In addition, the helical gear stage should contribute the largest share to the total gear ratio.
Next to the findings on reducing the mass, the arrangement of the individual stages results in further criteria that must be taken into account in the design process. Firstly, the ring gear diameter of the second and third stage should not be larger than the diameter of the first stage, since the following ring gears are to be designed as driven gears. The schema and the CAD model (Fig. 3) show that, in addition to the hollow shaft for supplying and controlling the pitch drives, the shaft of the first sun must also be guided through the sun of the second stage. The objective of a design that is as compact as possible also results in high demands on the constructive realization regarding mountability. These boundary conditions lead to an iterative process for determining the number of teeth and the total gear ratio. In this stage of the design process, commercial software is used to determine the safeties of the individual gear stages. In addition to verifying the available design space, the availability of bearings is a major challenge in the design process.

Assembly of the simulation model
The resulting first design of the gearbox components is used to develop a multibody-system model of the gearbox with In the models for the first, second and third stage the planets are supported by the carrier using force elements to describe the bearing characteristics. The carriers are supported in a rigid dummy body for the gearbox housing. The model for the fourth stage describes the dynamic behavior of the gear stage with the wheel shaft and pinion shaft, modelled by elastic beam elements, which are supported in a dummy body for the gearbox housing und connected by a force element to model the contact conditions in the helical gear stage [10]. The main model combines all substructures of the gearbox, defines the connection between the torque transferring components of the substructures by force elements and merges the dummy bodies to one body to represent the gearbox housing. This body has also six degrees of freedom, is supported by force elements on the mainframe and the mainframe is connected to the tower top. The gearbox housing and the mainframe will be replaced by flexible structures in the further development of the model.
The available model represents the NREL baseline model, extended by the described model of the gearbox, by an elastic beam element to represent the properties of the main shaft and by a generator model, connected by a coupling to the gearbox (Fig. 4), [10]. The complete model of the wind turbine can be used to simulate different load cases and allows the determination of the torque, the bending moments, and forces at the hub, but also at each modelled component in the gearbox. The resulting loads, simulated for different load cases and classified according to an assumed load case frequency, allows the definition of an approach for the design loads.

Modelling of flexible components
The characteristics of the planetary gear stages are also influenced by the geometry and the resulting stiffness of the planet carrier, whereby a stiffer carrier design increases the weight. In the first steps a sufficient safety factor against fatigue damages is used as design criteria and is the basis for an interactive process between the design in CAD, the stress determination using finite-element models and the determined design loads [3]. For the later integration in the multibody-system model, the last state of the resulting finite-element model gets modal reduced. This model represents the characteristic of the flexible structure by 100 mode shapes and allows the load introduction at specified load introduction points. These points must be defined in the finite-element model for the planet carrier bearings, for the joint connection points of the planets and the position of the bearings for the planets. In reality the load of a bearing gets transferred between the rings by the roller elements, which are moving in the zone of the introduced load. The load leads to a deformation of the component or the gearbox housing in the loaded area. To describe this characteristic in the flexible structure under consideration of the required linearization of the modal reduced model only a limited number of possible solutions is available. In the multibody-system model forces and joints with reference to a flexible body must be linked to discrete points called marker. For forces and joints a node on the axis of rotation must be defined in the finite-element model. The transfer of loads from this node to the surrounding flexible structure occurs by elements called multi-point-constraints. Currently different approaches are in use.
A first procedure connects the center point to nodes on the outer or inner diameter of the bearing seat in two planes which are positioned at both sides of the bearing. The resulting double-cone-shaped structure can act as a rigid connection which causes a stiffening of the structure (NAS-TRAN: RBE 2 elements). Using the second approach does not cause a stiffening. The displacement of the center point is dependent on the displacement of the connected points on the structure surface but no stiff connection between Hence the elements transferring pushing and pulling forces so that not only the housing stiffness of the loaded bearing seat area is considered (RBE 3). Elements which transferring only loads in one direction cannot be linearized by the modal reduction. But a behavior which is comparable to the described nonlinear elements can be modelled, if the center point is connected by a certain number of MPCs distributed over the circumference. They are only connected to nodes, which are in a specific angle range of the bearing seat. The transfer of pushing loads and the disconnection between load introduction and supporting nodes under pulling conditions can be described by a nonlinear bearing characteristic in the multibody-system model. At least the approach consists of a certain number of multi-point constraints (MPCs), characteristic curves for the bearing stiffness and force elements. The specified nodes, the number of considered mode shapes and additional nodes which are used to represent the deformation of the elastic structure with a higher level of detail are used to create the modal reduced finite-element model (Fig. 5).
The main shaft, the gearbox, the generator, and all electrical devices are supported by the mainframe. In many applications, a welded design is used for this purpose. The design criteria are driven by the required high stiffness, low weight and sufficient durability. The developed construction based on common designs for mainframes, the positioning of the sheet metal, the thickness, definitions of radii and cutouts is done under welding aspects and to ensure sufficient safety factors [2]. The assembly of the finite-element model is performed by usage of mid surface elements, the main bearing support and the torque arm supports are modelled by volume elements and the load introduction points are represented by the described multi-point constraints (Fig. 6).

Consideration of flexible gearings
The first state of the multibody-system model gets enlarged by the implementation of the assembled and modal reduced finite-element models. All connecting points for the bear-  Table 1 Natural frequency for different levels of detail of a wind turbine model  of the ring gear changes due the load introduction by the planets whereby in rotation direction of the carrier the ring gets elongated, behind the planet the ring gets compressed so that the maximum curvature occurs behind the planet and the deformation additionally influences the load distribution. To represent these effects in the multibody-system model also the gears can be integrated as modal reduced finite-element model (25 mode shapes are used) and using the force element FE-225 in SIMPACK (Fig. 7). Next to the described load introduction points on the axis of rotation, on each tooth, distributed over the flank width MPCs must be defined. The modelling of the elements and connection to nodes on the pitch diameter allows the description of the elasticity of the tooth by the modal reduced finite-element model. If nodes in the tooth root are used, the flexibility of the tooth gets calculated by analytical approaches in the multibody-system model. Over the flank width a minimum of three load introduction points should be specified. A higher number of MPCs increases the accuracy. The sun gear is implemented with the outer and inner geometry of the wheel, the ring gear is currently connected to the rigid structure of the gearbox housing.

Analysis in the frequency and time domain
The resulting detailed multibody-system model of the wind turbine considers the complete flexibility of the drive train and allows a detailed analysis of design loads and the influence of each single increase of the level of detail. Further it is theoretical possible to analyze and optimize the load distribution in the gear stages in detail.
The comparison of the torsional natural frequencies in dependency on the model level of detail shows the influence of the different components on the results. The torsional mode shapes of the drive train are mainly influenced by the degrees of freedom of the gearbox housing and the flexibility of the shafts and rotor blades (Table 1), [13].
Further investigations in the frequency domain focus on the excitation behavior of the gearbox considering imbalances, misalignments as well as tooth meshing frequencies and serve to localize possible critical operational ranges. The Campbell diagram cannot be used to identify critical speed ranges because the natural frequencies are very close together. The simulation of a slow run-up processes with the simulation model with the highest level of detail provides more precise information about possible resonances. Exemplary shows Fig. 8 the axial velocity of the first planet carrier over the generator speed and the frequency. The maximum of the axial velocity is reached at 1000 rpm and approximately 120 Hz as well as at 1200 rpm and 150 Hz, caused by the gear meshing frequency of the third planetary gear stage in the first order.
The first investigations in the time domain serve to validate the load assumptions for nominal bearing loads used for the design of the gearbox. The diagram (Fig. 9) shows a comparison of analytically determined bearing loads and the corresponding simulation results. On the generator side, good agreements can be proven, and the assumptions for the planet bearings can also be made reliably based on the torque. Only the assumptions for the planet carrier bearings differ significantly from the assumptions, whereby these values must be checked again by the load calculation with wind field.
The influence of the used modelling approach to represent the gearing can be seen for a more detailed analysis of the resulting contact forces over the width of the gearing and over the angle of rotation of the planet carrier. In Fig. 10 and 11 the results for the contact between sun and planet are shown in the left diagram, for the contact between planet and ring gear in the right diagram. Figure 10 shows the results for flexible modelled planet carrier, Fig. 11 the load distribution for the rigid modelled planet carrier. In the resulting plots the influence of the flexibility of the planet carrier can be especially seen at both ends of the gearing. The loads for the contact between sun and planet are higher and the load distribution gets shifted to the non-rotor end side (width defined with zero) if a rigid instead of flexible planet carrier is used. For the contact between planet and ring a comparable increase of the load can be seen.
Based on the simulation model, the effects of the modelling approach and of the rotor loads can be compared. The values for tooth flank modification (1125 µm) and crowning (100 µm) specified for shear torsional loading remain unchanged for the following investigations. The resulting face load distribution factors for the different analysis are shown in Fig. 12. The above described changes of the load distribution for a simulation model loaded under torque with flexible or rigid modelled planet carrier can be seen in the left two columns. The load distribution factor increases for the contact between sun and planet from 1.73 to 2.28 and for the contact between planet and ring gear from 1.63 to 2.21. It should be mentioned at this point that a further optimization of the load distribution is possible for the corresponding load cases and that the values here are intended to illustrate the basic progress which is influenced by the modeling approach. The resulting changes prove the large influence of the planet carrier stiffness. The comparison to the resulting load distribution factors for the model with flexible planet carrier and rigid modeled gearing shows the minor influence of the flexibility of the teeth.
In the fourth to sixth column the simulation models loaded with torque and a bending moment, and the results are again shown for different levels of detail of the modelling. It becomes obvious that under combined loads the modeling approach has less of an impact on the already high value of the load distribution factor between planet and ring gear. Also, the influence of the flexible modelled planet carrier is of lower importance. The dead weight of the rotor gets supported by the wind loads which reduces the inclination of the planet carrier in the gear contacts. The resulting load distribution factors for the rigid modeling of the gearings and the planet carrier show a significant increase for the contact between sun and planet which can be traced back to the applied tooth flank modification that take the twisting of the flexible sun gear into consideration.

Conclusion
The presented modelling approach and the results for analysis in the frequency and time domain show various possibilities to evaluate the influence of different parameters on the system properties. Based on the results for determining the torsional natural frequencies, it can be seen that almost every step taken to improve the level of detail of the model is associated with a reduction in the natural frequencies. Only the influence of the elastically modeled gearing on the resulting torsional natural frequencies can be assessed as low. In order to analyze the load distribution in a tooth contact of the first planetary gear stage, it is necessary to take into account all elasticities of the torque-carrying components and the support against the mainframe. The influence of the elastically modeled planet carrier appears especially significant in the case of sheer torque loading. The investigations with all load components put this finding into perspective and show that especially the load distribution factors for the tooth meshing between the sun and the planet require a high degree of detail in the simulation model. In the next step, the shown investigations must be extended to a multitude of load cases with simultaneous optimization of the gear modification. The determination of an optimum mod-ification for the entire operating range could be carried out using the frequency of occurrence of a load case and the utilization of the gearing resulting from the load distribution. The minimization of the utilization by adapting the gear modification would be conceivable as a target value. The influence of the so far unconsidered elasticities of the gearbox housing and mainframe must be examined in further investigations and the results have to be validated with detailed finite element investigations.
Funding Open Access funding enabled and organized by Projekt DEAL.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4. 0/.