Simulation methodology for the identification of critical operating conditions of planetary journal bearings in wind turbines

The usage of journal bearings as planetary bearings in wind turbines instead of roller bearings has become more common in recent years. Their usage is advantageous, due to smaller installation space needed compared to roller bearings allowing for higher power densities of wind turbine drive trains. However, this technology presents a challenge since there is currently no standardized approach for the design of planetary journal bearings regarding wear. Due to varying wind speeds and dynamic operating events a large variation of loads has to be considered in the design process of a planetary journal bearing for wind turbines. Some of these loads are considered potentially critical to the journal bearing in terms of wear. Identifying these critical load areas early in the design phase supports a reliable bearing design and wind turbine operation. This paper introduces a method to identify critical operating conditions for planetary journal bearings using a simulation tool chain, which couples a multi body simulation (MBS) model of a wind turbine with an elasto-hydrodynamic (EHD) model of the planetary journal bearing. Based on the EHD results critical operating conditions are determined for the planetary bearing. Furthermore, methods are implemented to reduce the number of required EHD simulations for analysing the bearing design. The combination of the identification of critical operating conditions, while reducing the computational effort leads to a simulation methodology, which enables a faster bearing design assessment considering the wide variation of wind turbine operating conditions. The applicability of this method is demonstrated by a simplified use case. Firstly, this paper introduces the MBS model and the parameter space that describes possible combinations of bearing loads such as forces, moments and rotational speed. Due to the number of combinations and the EHD computing effort, the identified parameter space is secondly sampled statistically to reduce the simulation effort. A risk map is derived from the EHD results, to easily indicate potentially critical operating conditions for the planetary journal bearing.


Introduction
Power production from wind energy is important for the renewable energy supply. More than 20% of the electric power consumption in Germany is fed by energy generated by wind turbines and this trend is increasing [1]. To further reduce the levelized cost of energy (LCOE) even more compact wind turbines with higher outputs are desired. This leads to higher torques, which need to be transferred in the planetary gear stages of the gearbox, while not increasing its installation space or weight [2]. In order to achieve this, the usage of journal bearings has become more common, since they require less radial space than the conventionally used roller bearings. Furthermore, journal bearings are less susceptible to typical failure mechanisms that occur in roller bearings such as fatigue or white etching cracks [3,4]. However, the usage of planetary journal bearings is challenging regarding a wear-safe operation and design. Firstly, the planetary bearings need to operate with higher specific loads at lower sliding speeds then are common for most journal bearings. Secondly, the tilting moments caused by the helix angle of the planetary gears leads to a risk of high edge pressures in the bearing surface [5]. These influences can lead to high wear rates and result in premature failures and thus should be especially considered in the design process of the planetary journal bearing.
Journal bearings rely on their hydrodynamic oil film for a wear free operation. However, due to high loads, low speeds or varying temperatures the thickness of the oil film can decrease, resulting in an operation within the mixed friction regime. This is characterized by asperity contact, which leads to wear of the bearing. The amount of wear occurring during operation is not only dependant on the asperity contact pressure and speed, but also on several other influences, such as material and lubricant properties [6]. For the prediction of the occurring wear several mathematical models have been developed over time. The use of different wear models and their factors is not part of this publication. This work focuses on the identification of operating con-ditions in wind turbines, which lead to the occurrence of asperity contact in the planetary bearings. However, to indicate the criticality of the occurring asperity contact, the ARCHARD's wear model [7] is used.
The planetary journal bearings considered in this work are located in the first planetary gearbox stage. These bearings operate at maximum speeds of 60 rpm. With maximum specific pressures of about 12 MPa at nominal operation [8].
Due to varying wind speeds and dynamic operating events such as turbine start, idling, noise reduced operation and turbine stop a large variation of loads needs to be considered in the design process of a planetary journal bearing. Some of the aforementioned operation modes are considered potentially critical to the journal bearing in terms of asperity contact and thus wear. The number of loads which need to be considered lead to a high calculation effort for designing a planetary journal bearing. In this design process, computational calculations are usually needed to capture all structural, dynamic and elasto-hydrodynamic effects. Therefore, multi-body simulations (MBS) of the wind turbine are often used to derive the loads acting upon the planetary journal bearing. The journal bearing is then analysed in detail by an elasto-hydrodynamic (EHD). The combination of detailed EHD and MBS simulations on the system level leads to even more calculation effort. For that reason, hydrodynamic (HD) bearing models are often used for quicker evaluation [9]. Lower computational effort is thereby achieved at the cost of accuracy by neglecting elastic deformations of the planetary pin and gear wheel. However, these elastic deformations are of importance for the load distribution inside the gearbox [10] and the pressure build-up in the journal bearings [11]. For this reason, recent research focuses on simulations of journal bearings taking these elastic deformations into consideration [12,13].
So far, this research focuses on the numerical simulation of the planetary journal bearings in wind turbine by considering different aspects, e.g. temperature, profiling, mixed friction, wear and deformation. To validate these simulations only a limited number of operation modes are needed and thus considered. Manufacturers consider in design process of the planetary journal bearing generally more of these operation modes [14]. The number of operation modes which is needed to capture all critical operating areas of the planetary bearing can however differ per bearing design. Therefore, this design process should not be limited to a standard set of load cases, but the entire operating range should be described. Considering this entire operating range leads in return to a high calculation effort. Therefore, a method which captures the critical operating areas by considering the entire operation range efficiently is needed. For which, this paper introduces a simulation-based method for planetary journal bearings in wind turbines. In here, a large variety of operation modes is considered, so the critical operating conditions can be identified. This criticality is visualized in a so-called risk map, which represents the wear risk over the applied load. Next to this, statistical methods [15] are implemented to describe the operation modes homogenously and reduce the number of considered operation modes.
The considered operation modes in this work are derived from an MBS model of a Vestas V52 wind turbine, which is coupled to a detailed EHD planetary bearing model. At the resulting identified critical operating areas for the planetary bearing measures can be taken to reduce operation or bearing design changes can be implemented. Furthermore, the effect of a bearing design change on the critical operating areas is demonstrated.

Method
The simulation methodology is further explained by giving an overview of the coupling of the several simulation models, the method to reduce the simulation effort and the identification of critical operating conditions. The used model, the derivation of bearing loads and the model parameters are described in Sect. 3. The described methodology is demonstrated by using a simplified use case, which only considers the steady production operation.

Simulation overview
The tool chain used in this work consists of an MBS [16,17] system model of a wind turbine that can be co-simulated with a model of the turbine controller. This allows for the evaluation of realistic turbine dynamics. The loads considered in this work are taken from simulations of standardized design load cases for wind turbines according to [18]. The MBS delivers the loads acting on the planetary journal bearing. These loads are applied to an EHD simulation model of the planetary bearing. An overview of the entire tool chain is provided in Fig. 1.
The MBS model of the turbine is fed with forces on the rotor blades simulating wind inflow. These forces are applied using a force-element for rotor aerodynamics based on the blade-element method (BEM). BEM is a widely used method for calculating aerodynamic loads such as thrust and drive torque on wind turbines, and is much more computationally efficient than corresponding CFD calculations [19]. Wind flow speed, direction and turbulence can be parametrized in the model and hereby a realistic environment of the turbine can be simulated. For the sake of simplicity this work focuses purely on laminar and constant wind flow. The loads caused by the windfield are then transferred by the rotor to the gearbox, in which all planetary gears are included. The planet carrier is here condensed as a flexible body by means of the finite-elementmethod (FEM). From this model the forces at each planetary gear can be extracted. Which are fed in the EHD model of the planetary journal bearing, which is build using the AVL Excite MBS software. This model accounts for temperature and pressure dependent fluid properties, cavitation, heat transfer, thermal expansion, structural deformation and hydrodynamic as well as asperity friction work [20,21]. To calculate the pressure in the fluid film of the planetary journal bearing the Reynolds equation is solved, which is extended by flow factors according to Patir and Cheng to consider roughness effects on the pressure and shear flow. The calculation workflow and solved equations are further explained in [20,21]. The flexible planet carrier is included in the MBS-model, therefore the EHD model consist of two condensed flexible bodies of planetary gear and pin. In the EHD model the bearing loads are applied at the planetary gear, which results from the system simulation. The bearing loads consist of a radial load and tilting moment. The EHD model will be further explained in Sect. 3. The validation of this model is currently on-going and only used in this work for demonstration of the simulation methodology. Finally, the EHD results of interest are extracted to perform the wear calculation and create the risk maps.

Reduction of simulation effort
The possible load combinations are wind turbine specific and based on the occurring wind speeds for each location. The wind speeds are site specific and constrain the possible load combinations. The area between these constraints is referred to as the simulation space. An overview of this simulation space based on the torque and rotational speed of the turbine driveshaft is given in Fig. 2.
In this work, only the steady production operation conditions are considered. For these operating conditions the total number of possible bearing load combinations, if integer speed and torque values are considered, is 5720. This leads to high computational effort and therefore sampling techniques are implemented. The simulation space needs to be described as evenly as possible to capture all possible critical operation areas, such approaches are referred to as space-filling designs. Latin Hypercube sampling (LHS) combined with the maximin criterium is a sampling technique which is a most commonly used class for space-filling designs [15] and is therefore used in this work. The number of samples needed to identify the critical operating condi- Simulation Space Samples Fig. 2 The constrained simulation space and high-density evenlyspaced sample points tions is not initially determinable, since the occurrence of asperity contact is a highly non-linear problem influenced by numerous effects. At first, an evenly spaced sampling is implemented with a high density of samples [22]. In this case a sample size of 134 bearing load combinations is used, this leads to step widths of 2 rpm and 20 kNm. These samples are indicated in Fig. 2. The evenly spaced sampling has the advantage that evaluation points can easily be removed to get an evenly spaced sampling with a smaller sample size. For a random sampling the evaluation points locations are very dependent on the used sampling size. With the evenly spaced sampling it is intended to get a minimum sample size which fills the simulation space efficiently and leads to a similar risk map as with a higher sample density. When the needed sample size is identified, the space filling sampling is implemented, which uses LHS with the maximin criterium, to maximize the distance between the sampling points and uniformly distribute these [15]. Using such an algorithm can reduce the needed sample density further and leads to the minimum number of load combinations needed to identify the critical operation condition areas.

Critical operation conditions
Critical operation conditions are defined as load combinations in which asperity contact occurs between the journal bearing and the planet. However, the severity of wear volume removal due to asperity contact is dependent on several other factors. Therefore, in this work the criticality is determined based on ARCHARD's wear model [7,21]. This model is selected, since it's commonly used for journal bearings and standards are in place to derive the required factors [23]. In ARCHARD's wear model material is removed from the surface at areas where asperity contact occurs. The wear K load is calculated by Eq. 1, from which the wear depth hv follows according to Eq. 2.
In where pc is the asperity contact pressure, Tacc is the total accumulation time, u is the relative sliding speed, T is the cycle time, H is the hardness of the material and k is the wear factor. The total wear volume is determined by the wear depth and the asperity contact area. The wear factor is constant for all load combinations [7], just as the hardness. Since static operation conditions are considered, the cycle time is constant as well. The accumulation time Tacc is based on the probability of occurrence for the regarded operating point, which is derived from a statistical wind speed distribution. The sliding speed is derived from the rotational speed. The asperity contact pressure and contact area is a result of the EHD analysis. The criticality of the operating conditions is therefore dependent on the sliding speed, contact pressure, contact area and the accumulation time. The EHD contact is represented in AVL Excite as a surface mesh. At each node the results are extracted, the wear volume for each node is calculated and summed to derive the total wear volume for each operating point. This is done according to Eq. 3.
Where A is the elemental area, n the number of nodes and Vv the total wear volume. To summarize this information in a risk map, the data is normalized so that a value of 0 indicates safe operation conditions and at values larger than 0 asperity contact will occur. Due to the range of speeds, probabilities and the resulting asperity contact areas and pressures the differences in the occurring wear volume is large. Therefore, a logarithmic normalization is used, so that large numbers have less effect on the risk map and thus operating points with a low wear volume will still be visible in the contour. The wear risk is then calculated for each operating condition according to Eq. 4. Here Rv is the risk value ranging from 0 to 1. The wear volume is here a function of the operating points, the maximum and minimum wear volume are used as well by extracting these values from all the operating points.

Model description
As previously described, in this work an MBS-based wind turbine system model and an EHD-model of a planetary journal bearing are used to identify wear critical operating conditions (see also Fig. 1). The gearbox of the wind turbine is originally designed with planetary roller bearings and a possible journal bearing variant of this gearbox is developed for this work. The selected geometry of the journal bearing is intended to generate asperity contact at rated speed. Therefore, the B/D ratio is kept relatively small compared to state-of-the-art bearings. This leads to an intentionally increased wear risk due to the higher specific bearing load and more risk to high edge pressure due to lower bearing widths. The system parameters for both models (wind turbine and planetary journal bearing) are provided by Table 1. Flexible tooth contacts are neglected in this work for the sake of simplicity. Thus, ovalization effects of the planetary gear is also not considered. The inclusion of these deformation effects into the toolchain is on-going, however for now excluded for the demonstration of the simulation methodology.  Fig. 10 & Fig. 11. The journal bearing is profiled with a quadratic recession, the profile is shown in Fig. 3. The loads are applied to the planetary journal bearing according to Fig. 4. The circumferential tooth force and tilting moment, resulting from the helical gearing of the planet, are applied as a lumped force at the centre of the planetary gear. Which are distributed by means of a distributed coupling to the entire planetary gear. The planetary speed is applied, as well, at the centre of the planetary gear. The pin is presented as the stationary body, by restricting the degrees of freedom at both ends.
From the MBS-model, the planetary speeds and bearing loads are extracted. By doing so, the journal bearing load cases that correspond to specific wind speeds, can be derived from the presented simulation tool chain, as schematically shown in Fig. 1 and presented in Fig. 4. In this work, only the static load cases are considered. These are derived by considering the probability of occurrence of stochastic wind speed distribution for a specific wind turbine location. The IEC standard for wind turbines specifies the design load cases (DLC) that have to be considered during the turbine design process [18]. The design load case of steady production operation (DLC1.1) is often referred to in the evaluation of drive train loads.
In a first step, MBS simulations were carried out for stationary production operation at wind speeds at hub height in the range between the minimum and maximum wind speed of the turbine. A location specific wind speed distribution for middle Germany is used to derive a distribution of load cases (torque and rotational speed of the turbine's main shaft) and a corresponding probability of occurrence for these load cases [17]. In Fig. 5a the used Weibull distribution of the on-site wind speed at hub-height and the resulting load-case distribution is shown. In the following only the speed-torque combinations with a probability greater Asperity contact pressure area [%] 6 Rotor torque [kNm] 10 10 200 Rotor rotational speed [rpm] 14   than zero are considered in the EHD study, which are the coloured bars in Fig. 5b. The derived loads cases are sampled, as described in Sect. 2.2.

EHD results
In Fig. 6. the results of the EHD planetary journal bearing simulations are shown. To identify the critical operating conditions, the most interest is in the asperity contact pressure and area. To understand the bearing behaviour the maximum oil film pressure and minimum film thickness are included in Fig. 6 as well. The maximum or minimum value is extracted for each sample point, as shown in Fig. 2, and plotted in the bar graph in Fig. 6.
The results show that two diverging operating areas lead to potential wear critical conditions in the planetary journal bearing. Asperity contact occurs both at low (5 rpm) and high speeds (22 rpm). At these two exemplary operating points the resulting distribution of asperity contact are shown in Fig. 7, the corresponding operating points parameters are given in Table 2. For higher speeds, the applied loads increase, as is indicated in Fig. 6. This leads to higher tilting moments and thus increasing edge pressures, which leads to wear only occurring at the bearing's edges as is shown on the right in Fig. 7. At low speeds, journal   Fig. 6d. The location of the asperity contact changes, therefore wear occurs more towards the bearing centre as well as is shown on the left in Fig. 7. The asperity contact area is here larger, which is shown in Fig. 6b. As expected, the selected bearing geometry leads to a large number of operating conditions at which asperity contact occurs due to the comparably low B/D ratio of the bearing (see Fig. 6). This results in distinct contour gradients in the risk map for the initial sampling (high-density evenly-spaced), which is depicted in Fig. 9a. As described in Sect. 2.2, the number of sampling points is reduced by factors of two relatively to the initial sampling until the contour changes. Reducing the initial sample points by a factor of 4 leads to a contour change, which is indicated in Fig. 9b. The sampled points are shown in Fig. 8b, the  number of samples is 32 for this set. This number is also used for the Latin hypercube sampling, which is shown in Fig. 8c. The same number of samples is used to compare the performance of Latin hypercube sampling versus evenlyspace sampling. The Latin hypercube sampling matches the initial sampling better, as is shown in Fig. 9c. The risk map given in Fig. 9 indicates the criticality of the operating conditions. The wear risk is in here scaled between 0 and 1, a value of 1 indicates the operating condition with the highest wear. The initial sampling, which is used as reference, consisting of 132 samples shows as expected the clearest risk gradients. The reduced-evenly spaced sampling and Latin hypercube sampling both use 32 samples, of which Latin hypercube sampling shows similar critical operating areas compared to the initial sampling. Out of the 32 simulations at 18 of them asperity contact occurs for the Latin hypercube sampling and at 15 of them asperity contact occurred for the reduced evenly-spaced sampling. Thus, the Latin hypercube sampling captures more critical conditions.
To indicate how different designs can affect the risk map, the simulations have been performed for a design in which the edge recession of the bearing profile is reduced to 50 µm. This result is shown in Fig. 11b. This design change leads to more operating conditions at which asperity contact occurs.
The differences in results for one operation point between the two profiles is shown in Fig. 10, the selected point is marked with an asterix in Fig. 11.

Discussion and conclusions
This work presents a method to derive critical operating conditions for planetary journal bearings. Considering the loads and constrains on the system level and the detailed component level. With the detailed bearing model, the asperity contact is extracted, the wear volume is calculated and the risk map is created. The risk map indicates the operating conditions at which asperity contact occurs. For the considered bearing model asperity contact occurs at many operation points, due to the comparably small B/D ratio. The bending moments, which occur due to the helix angle in planetary gears, can lead to asperity contact at the bearing edges. However, asperity contact occurs also more towards the bearing centre, due to the low rotational speeds.
Methods to reduce the number of simulations are implemented and lead to visually comparable risk maps. In this use case the number of simulations could be reduced by a factor of 4, leading to 32 simulations for the risk map generation. It can be concluded that, using the Latin hypercube sampling method is advantageous as it describes the possible load combinations homogenously and therefore captures more critical operating conditions. Small changes in the design can lead to different critical operating conditions, which has been demonstrated by changing the quadratic profile of the bearing. A geometric change can lead to more operating conditions in which asperity contact occurs. Next to this, the location at which the asperity contact occurs can change as well. This should demonstrate the necessity of considering a wide range of operating points in the design phase, since small bearing design changes can lead to different critical conditions. An extension of this work to account for the running in process in the bearing is on-going. It is expected that after running in the amount of asperity contact reduces and can for some operating points even disappear. The presented tool chain can be extended to account for the running in behaviour as well. This also applies for the dynamic operating conditions, in this work only the applicability of the method has been shown based on the static operating conditions. However, the changing torque and speed are likely to have an effect on the asperity contact as well. Identifying the possible torque and speed gradients plus the sequence in which these occur should be extracted from the occurring wind speed information. An extension of this work to identify dynamic operating conditions and their influence on the bearing criticality is also on-going.