The impact of geometric non-linearities on the fatigue analysis of trailing edge bond lines in wind turbine rotor blades

The accurate prediction of stress histories for the fatigue analysis is of utmost importance for the design process of wind turbine rotor blades. As detailed, transient, and geometrically non-linear three-dimensional finite element analyses are computationally weigh too expensive, it is commonly regarded sufficient to calculate the stresses with a geometrically linear analysis and superimpose different stress states in order to obtain the complete stress histories. In order to quantify the error from geometrically linear simulations for the calculation of stress histories and to verify the practical applicability of the superposition principal in fatigue analyses, this paper studies the influence of geometric non-linearity in the example of a trailing edge bond line, as this subcomponent suffers from high strains in span-wise direction. The blade under consideration is that of the IWES IWT-7.5-164 reference wind turbine. From turbine simulations the highest edgewise loading scenario from the fatigue load cases is used as the reference. A 3D finite element model of the blade is created and the bond line fatigue assessment is performed according to the GL certification guidelines in its 2010 edition, and in comparison to the latest DNV GL standard from end of 2015. The results show a significant difference between the geometrically linear and non-linear stress analyses when the bending moments are approximated via a corresponding external loading, especially in case of the 2010 GL certification guidelines. This finding emphasizes the demand to reconsider the application of the superposition principal in fatigue analyses of modern flexible rotor blades, where geometrical nonlinearities become significant. In addition, a new load application methodology is introduced that reduces the geometrically non-linear behaviour of the blade in the finite element analysis.


Introduction
In order to capture enough wind energy modern multi-MW wind turbines are equipped with large rotor blades with lengths of 80 m and more. Such blades are slender and flexible structures, and they experience large deformations in flap-and edgewise direction. These deformations provoke geometrically non-linear effects when it comes to ultimate strength analysis in the design stage [1]. The computational costs for a single geometrically non-linear static ultimate strength analysis is relatively low. However, geometrically non-linear effects within a fatigue load history will result in huge computational efforts in the framework of transient stress analyses for the calculation of accurate stress results.
This paper is concerned with geometrically non-linear effects of 3D finite-element-based stress analyses in the context of fatigue lifetime predictions of trailing edge bond lines in large rotor blades. The simulations are based on the reference wind turbine model IWT-7.5-164 [2], which is a 7.5 MW direct drive wind turbine with blades of 80 m length. The analysis is further focused on the shear stress proof according to the 2010 edition of the Germanische Lloyd (GL) certification guidelines [3] and a multiaxial stress proof as demanded by the latest DNV GL standard [4]. Both standards disregard any potential geometric non-linearity of the blade arising during fatigue analyses.
Geometric non-linearities may appear when the strains in a material exceed more than a few percent (large strain), or when large deflections result in additional offsets of loadings, and thus in additional internal moments and stresses [5]. Modern wind turbine rotor blades are subjected to large structural deformations [6], not only in ultimate load situations but also in fatigue time histories. Since the accuracy in stress response predictions is crucial for a correct estimate of the fatigue lifetime, it is of utmost importance to take into consideration geometric non-linearities in the design phase of blades, as will be pointed out in the framework of this paper. Otherwise, a structural blade design suffers from underestimations of the structural stresses, and consequentially from a nonconservative overestimation of the fatigue lifetime. The consequences will be severe damages up to failure of blades and probably of entire turbines in real life.

Aero-servo-elastic turbine simulations and finite element model of the blade
The loads for the fatigue and strength analysis are calculated through an aero-servo-elastic simulation that couples a multi-body system of several non-linear beam models [7] for the entire structure of the turbine with an aerodynamic model of the rotor. In this work, we utilized the turbine simulation software HAWC2 [8] to compute the loads. The relevant output for the fatigue analysis of the blade are internal forces and/or deflections in different cross sections along the blade, taking into account large deflections of the blade and the tower.
For the stress analysis a finite element (FE) model of the IWT-7.5-164 reference blade [2] is utilized. The FE model is generated with a Model Creator and Analysis tool called MoCA that has been developed at the Institute for Wind Energy Systems at Leibniz Universität Hannover. MoCA is encoded in MATLAB [9] and creates fully parameterized three-dimensional finite element models of rotor blades for ANSYS [5]. The stress post processing in this paper is limited to a span-wise region between 43 m and 47 m. A mesh convergence study has been performed to guarantee reliable results. The resulting global FE mesh and a local view of the adhesive discretization are shown in Fig. 1.
The boundary conditions have to be applied to the FE model in a reasonable way. The application of internal moments to discrete points along the blade is impossible without cutting the blade into sections, and would result in a stepwise increase of the moments along the blade. Therefore, it is common practice to apply external lateral single substitute forces Q sub,j at different positions along the blade that are denoted by r j , see Fig. 2   Q sub,j result in a stepwise linear distribution of bending moments and are adjusted to obtain a good approximation of the real bending moment distribution, see Fig. 2 (c). This methodology is similar to the loading strategy in a full-scale blade test, which is depicted in Fig. 2 (a). Note that with this load application the applied transverse force distribution is not well approximated, since single forces yield stepwise jumps, which are not realistic in a real model. The normal forces in span-wise direction and the torsional moment are not included in the load application. However, since the bending moments are the major design drivers, this load application is considered sufficiently accurate to reproduce the stresses and strains in the blade.
In order to distribute the single forces Q sub,j across the cross-sections at the respective positions r j along the blade, we introduce master nodes located in the shear centers in order not to create erroneous torsional moments. Those master nodes are connected to slave nodes in the cross-sections via rigid multi point constraints (MPC) in ANSYS. These MPCs represent the load clampings as seen in a full-scale test, see Fig. 2 [10] (a), bending moment approximation with substitute loads on a beam (b), and comparison of approximated and real bending moments plotted against the span-wise position along the blade (c).

Fatigue analysis with substitute forces
After computing the load time series for specific wind speeds and operation modes via aero-servoelastic turbine simulations with HAWC2, the stress histories have to be derived for all nodes in the finite element model. For this purpose the linear superposition principal is normally utilized, which states that the total response of a linear system is the sum of the responses of all contributing individual forces [11].
The fatigue analysis methodology based on the superposition principal works as follows: First, we define the number of single substitute loads n to be applied for the approximation of the bending moments. Then, we apply those loads, at this stage as unit loads, individually to the 3D finite element model. For each of the unit loads we carry out a finite element analysis and save the nodal stresses in an output file. Then we calculate the magnitudes of the substitute loads for each time step of the load time series that are required to have a sufficiently accurate approximation of the bending moments. The next step is to scale the stresses from the n unit load FE simulations by the real magnitudes of the substitute loads, and to superimpose all n scaled stress values for the entire time series. In this way we transfer the loads analysis to the 3D FE model. In other words, we transfer the internal forces and moments from beam theory in the framework of turbine simulations to stresses from 3D FE simulations. These stresses then enter a fatigue life prediction in a postprocessing calculus according to Palmgren & Miner theory, including rainflow counting, Markov matrix generation, Goodman diagram and S-N-curve evaluation, and linear damage accumulation. Note that for this methodology n detailed finite element simulations have to be carried out. Usually, 5 ≤ n ≤ 8, i.e. the number of substitute loads that is required to have a good approximation of the bending moments is quite low. Consequentially, the computational costs for the FE simulations is also limited to an acceptable degree. However, we would like to emphasize that geometrical and material linearity is required for the validity of the superposition principal. It is seen in the next section that this is a rather rough assumption, which is generally not valid for large and flexible rotor blades.

Linear versus non-linear finite element simulations
To evaluate the validity of the aforementioned methodology, a comparison of linear and nonlinear FE simulations is carried out. In this context we use the highest shear stress τ yz (the y-axis is pointing in chord direction of the cross-section, and the z-axis in span-wise direction) for the proof according to the 2010 edition of the GL guidelines [3] and the Drucker-Prager equivalent stress [12] for the multi-axial stress proof demanded by the latest DNV GL guideline [4]. The Drucker-Prager equivalent stress is calculated according to Herein, m = R t /R c is the ratio between the tensile and compressive strength of the adhesive, R t and R c , and σ 1 , σ 2 , and σ 3 are the principal stresses, respectively, following the convention σ 1 ≥ σ 2 ≥ σ 3 . The highest deviations between linear and non-linear stresses emerge at large deflections, thus high loads. For the example of the trailing edge bond lines the maximum edgewise load is assumed to be the critical load case, as the strains in the bond lines reach their maxima due to the large distance to the principal bending axis. Table 1 shows the load set for maximum edge wise bending during the simulation of the design load case (DLC) 1.2 according to IEC standard 61400-1 [13], representing the normal power production load case for fatigue life prediction.
To determine the influence of geometric non-linearities, the stresses in the bond line were computed performing both, geometric linear (lin) and non-linear (nlin) calculations using the load application methodology explained in section 3. In the following, the stresses of both calculation methods (linear and nonlinear) are compared. The fatigue proof of the trailing edge according to the GL certification guidelines uses the shear stresses in the adhesive to determine the material effort. Fig. 3 shows the mean deviation of the shear stresses from linear to nonlinear analysis and the range of the three-fold standard deviation for all finite elements in the corresponding bond line cross section. The relative deviation is related to the linear calculations, as these are considered as the state of the art. Disregarding the outer regions of the analyzed bond line, where free edge effects could affect the stress intensity, the mean value of the deviation from non-linear to linear FE-solution of the shear stress is approximately 48 %. Respecting the range of the threefolded standard deviation the value can even reach up to 80 %. For the estimation of the error of the fatigue life prediction according to the new DNV-GL guideline the Drucker-Prager failure hypothesis is used to evaluate the multi-axial stress state. This criterion is developed to model materials with different tensile and compressive strengths [14], which in fact is also the nature of commonly used adhesives (EP and PU) showing a significantly higher compressive than tensile strength [15]. Fig. 4 illustrates the deviation of the computed Drucker-Prager equivalent mean stresses in the bond line cross section with their corresponding three-fold standard deviation range. The deviations for the multi-axial stress proof differ from the shear stress proof. Mean values for the relative difference of geometric linear to non-linear calculation are only about 3 to 4 %. According to the graph in Fig. 4 the maximum expected values reach up to 9 %, disregarding the outer regions of the bond line.
The large deflections of the blade and the large distance between the evaluation region and the load application positions lead to a significant geometrically non-linear effect in the stress space. The shear stresses tend to suffer from geometric non-linearities to a high degree, see Fig. 3. The effect is smaller for the Drucker-Prager euqivalent stress (Fig. 4), which is dominated by normal span-wise stresses, see Fig. 10. Nevertheless, the non-linearity is remarkably high, and could even increase if the torsional moments and normal span-wise forces are applied in addition to the bending moments. We can thus conclude that the linearity assumption is a strong weakness in the fatigue life prediction methodology presented in section 3.

Direct application of internal forces and moments
In order to reduce the influence of geometric non-linearities from the load application, a new method is proposed to introduce the loads which is based on a sectionwise analysis of the blade. As we receive non-linear internal forces and moments from the aero-servo-elastic simulations these can be applied directly to the model for the purpose of modelling the boundary conditions as accurate as possible. Here the internal forces and moments are applied behind the analysis  region in span-wise direction to the tip. The analysis region is reduced to a short span-wise portion of the blade dR Analysis to guarantee a good approximation of the loads in this region. In order to avoid influences from the load application modelled with rigid MPCs a distance dR LA Dist is kept to the analysis region. A big advantage of this method is that all forces and moments can be applied as they were calculated in the loads analysis. The applied internal forces with the distance dR LA Dist to the analysis region would introduce additional bending moments in the structure, therefore the applied internal moments are reduced by the amount of the bending moment produced by the internal forces as shown in Fig. 5. The resulting bending moment of the superposition of the applied internal moment and the moment provoked by the internal force is triggered to the most conservative point at the edge of the analysis region that is closer to the blade root. Using this method the finite element analysis is only valid for the analysis region. Only there the internal forces and moments are accurately approximated. For the superposition principle in the fatigue analysis the necessary stress response can be derived from the results of the elements inside the analysis region. To perform the fatigue analysis for the whole blade the analysis region and load application has to be shifted repeatedly in order to obtain information for all blade regions. The larger the analysis region is the less computational effort is necessary to examine the whole blade. However, the larger the analysis region is the poorer the approximations of the internal forces and moments becomes. Hence, the user has to find a compromise between computational costs and accuracy.
In order to compare both load application methods we have chosen the same span-wise region from 43 m to 47 m for the analysis. All internal loads except the torsional moment and the force in span-wise direction were applied to obtain the same scenario as for the analysis described in section 4. Fig. 6 shows the deviation of the shear stress from geometric non-linear to linear analysis. The mean deviation increases in span-wise direction from 33 % to 43 %, but respecting the standard deviation the maximum value can reach from 45 % to 65 %. The deviation is very high and indicates a significant impact of geometrical non-linearity on the stress response. However, in comparison to the former load application method described in section 4 the new method reduces the mean deviation value by 5-15 % and the maximum by 15-35 % for the shear stress analysis of the highest edgewise loads in fatigue load histories. Regarding the Drucker-Prager equvialent stress for the multi-axial stress proof shown in Fig. 7 according to the new DNV GL standard [4] the mean deviation between the geometrically nonlinear and linear analysis is only about 2 % and the maximum expected deviation is up to 3.5 %, neglecting the outer regions of the bond line. These results show that the new application method can halve the deviation between geometrically non-linear and linear results in comparison to the former load application. The new load application method allows the application of torsional moment and spanwise normal forces. The inclusion of these additional loads intensifies the non-linear effect for the shear stress dramatically as shown in Fig. 8. The mean value of the deviation increases to 150 %. While the shear stress is very sensitive to geometric non-linearities introduced by torsional moments the Drucker-Prager equivalent stress for the multi-axial stress proof increases only slightly in its maximum mean deviation up to 3 %, see Fig. 9. This might be due to normal stresses in span-wise direction, which dominate the Drucker-Prager equivalent stress, being relatively insensitive to geometric non-linearity. However, more likely a low impact of geometrical non-linearity on the absolute values of all stresses, see Fig. 10, results in a dramatic relative change in shear stresses (SXY,SXZ,SYZ) (as the absolute values are very small), but a moderate relative change in Drucker-Prager equivalent stress (as that is governed by the normal stress in span-wise direction (SZ), of which the absolute value is orders of magnitude larger than the shear stresses).

Conclusion
This paper has demonstrated that large wind turbine rotor blades behave geometrically nonlinear in finite element analyses, not only in extreme load cases but also in fatigue load histories. This is a fundamental information that every blade designer has to be aware of, because the standard methods for fatigue life estimations are no longer valid. Neglecting the nonlinearities can result in nonconservative designs, and thus to severe damages in real rotor blades that may probably lead to the failure of blades or even entire turbines.
Depending on the type of stress proof for the bond lines in the rotor blade, the non-linear effects can have a high influence on the stress response, e.g. up to 80 % with a shear stress proof as shown in the paper. The utilization of a multi-axial stress proof, in this case based on the Drucker-Prager equivalent stress, reduces the impact of geometric non-linearities to a maximum of 9 %, at least in the analyzed trailing edge bonding region. A side effect is that the major failure mode, namely cracks due to normal stresses in span-wise direction, are captured -in contrast to former design guidelines that only asked for a simplified shear strength proof. It has further been shown that load introduction methodologies have a significant impact on the stress response. Although computationally very attractive, a methodology based on scaled and superimposed unit loads, as described in section 3, is principally not applicable due to the requirement of linearity for its mathematical and physical validity. A novel load introduction technique has been proposed where all internal loads can be considered. The novel concept significantly reduces nonlinearity effects in the stress space. The disadvantage is that the computational costs are quite high, because the evaluation regions are strongly limited in span-wise direction, and thus a sectionwise analysis has to be performed.
The loading state chosen for the investigations in this paper was taken from a normal operation simulation (design load case 1.2 according to the IEC standard), where an instance of time was chosen where the bending moments are maximum. For other time instances, the loads may be much lower, and the nonlinearity effect may be less significant. An evaluation of the impact of geometrical non-linearity on the fatigue lifetime prediction, and not just one instance of time, is ongoing. That study will include the scatter in loading over time, and will be published in the near future.
Principally, the results of this paper indicate that further investigations are needed in order to develop an accurate and computationally not too expensive load application methodology for finite element analyses in the context of fatigue life assessments for wind turbine rotor blades.