Journal topic
Wind Energ. Sci., 5, 503–517, 2020
https://doi.org/10.5194/wes-5-503-2020
Wind Energ. Sci., 5, 503–517, 2020
https://doi.org/10.5194/wes-5-503-2020

Research article 20 Apr 2020

Research article | 20 Apr 2020

# The effects of blade structural model fidelity on wind turbine load analysis and computation time

The effects of blade structural model fidelity on wind turbine load analysis and computation time
Ozan Gözcü and David R. Verelst Ozan Gözcü and David R. Verelst
• DTU Wind Energy, Technical University of Denmark (DTU), Frederiksborgvej 399, 4000 Roskilde, Denmark

Correspondence: Ozan Gözcü (ozgo@dtu.dk)

Abstract

Aero-servo-elastic analyses are required to determine the wind turbine loading for a wide range of load cases as specified in certification standards. The floating reference frame (FRF) formulation can be used to model the structural response of long and flexible wind turbine blades. Increasing the number of bodies in the FRF formulation of the blade increases both the fidelity of the structural model and the size of the problem. However, the turbine load analysis is a coupled aero-servo-elastic analysis, and computation cost not only depends on the size of the structural model, but also depends on the aerodynamic solver and the number of iterations between the solvers. This study presents an investigation of the performance of the different fidelity levels as measured by the computational cost and the turbine response (e.g., blade loads, tip clearance, tower-top accelerations). The analysis is based on aeroelastic simulations for normal operation in turbulent inflow load cases as defined in a design standard. Two 10 MW reference turbines are used. The results show that the turbine response quickly approaches the results of the highest-fidelity model as the number of bodies increases. The increase in computational costs to account for more bodies can almost entirely be compensated for by changing the type of the matrix solver from dense to sparse.

1 Introduction

Modern wind turbine blades are large, slender and flexible composite structures with a complex pre-bent and twisted geometry. Over their operational life blades undergo large deflections and rotations due to external loads (e.g., aerodynamic, inertial and control actuator loads). An aero-servo-elastic code or framework is used to accurately calculate the complex dynamical response of wind turbines with large and flexible blades. This has led to the implementation of geometrically nonlinear structural solvers in wind-turbine-specific aero-servo-elastic codes. For example, the structural solver BeamDyn was implemented in FAST . It uses the geometrically exact beam theory (Hodges1990) based on the Legendre spectral finite element method. Another example is a recent release of Bladed (DNV GL2018) that uses a multibody formulation to capture large structural deflections of the modeled structures. BHawC is another nonlinear aeroelastic wind turbine simulation code which uses a corotational formulation to resolve large deflections accurately.

In this study the turbine responses of DTU10MW and IEA10MW are considered with different structural fidelity levels of the blades for 432 load cases according to design load case (DLC) 1.2 . Deterministic load cases (without turbulent wind) are also considered to evaluate the turbine steady-state response at various wind speeds. The loads at different points on the turbine, controller activity and turbine performance are compared. Section 2.1 introduces the solver (HAWC2) and geometrically nonlinear structure modeling in the multibody (FRF) formulation. Section 2.2 presents the reference wind turbines, load cases and their models as used for this study. Section 3 includes the calculation methods used when post-processing the results, the plots of the computation time, steady-case results, DLC 1.2 blade results, DLC 1.2 tower and performance results, stability results, and a discussion of the results. The conclusions of this study are given in Sect. 4.

2 Method and analysis

Evaluating the aero-servo-elastic response of large and flexible wind turbines using time domain simulations under turbulent inflow conditions requires rigorous analysis. Both the aero-servo-elastic solver and the considered model and load cases are therefore carefully outlined in the following two sections. The applied analysis method presented here is based on a numerical experiment of blades with varying structural model fidelity levels.

## 2.1 Method

The turbine analyses for the presented work were performed with HAWC2 version 12.6, which is a strongly coupled aero-servo-elastic wind turbine simulation tool. The aerodynamic solver of HAWC2 uses the blade element momentum formulation including effects of dynamic stall, dynamic inflow, wind shear on induction, tip loss, tower shadow and large blade deflections. A proportional–integral–derivative (PID) controller algorithm is used to determine the set point of the pitch bearing angle and generator torque. The servo actuators are modeled as a second-order dynamical system with an appropriate given frequency and damping. The structural dynamics of HAWC2 are based on a multibody formulation using an augmented FRF method (Shabana2010). Each structural element has two nodes with 6 degrees of freedom (dof) and is modeled as a linear classical isotropic or anisotropic Timoshenko beam . A body, defined in the FRF formulation, can be composed out of an arbitrary number of elements. Bodies are attached to each other with constraints in any of the 6 dof (three rotations and three translations). The bodies are deflected linearly, but their body reference coordinate system follows the translation and rotation from the last node of the previous body in a continuous structure model.

A general wind turbine structure can be built out of Ne elements and Nb bodies with constraints, but NbNe. The constraints allow the user to capture the correct nonlinear geometrical response of a collection of bodies in a continuous structure as long as the deflections within one body are small . In the limit case where a continuous structure model has the same number of bodies as elements (Nb=Ne), the solution is equivalent to the corotational approach . For example, Fig. 1 shows how the body discretization of a 2D beam structure model captures the nonlinear effect on the beam length as bending deflection occurs. The beam model has nine linear beam elements. The round markers represent the finite element nodes, and the square markers represent the body discretization of the structure. As seen in the figure, the one-body model has linear deflections with fictitious elongation due to lack of large rotations, while the three-body model shows the large rotation effects due to constraints between the bodies.

Figure 1Structural modeling of a cantilever beam in floating reference system with multiple bodies, in deflected and undeflected states.

HAWC2 constructs a system of differential equations representing the equations of motion of the system with constraints (see Eq. 1) which is based on a given set of Ne elements and Nb bodies (Shabana2013) for the ith time step “ti”. $\mathbf{M}\in {\mathbb{R}}^{N×N},\phantom{\rule{0.25em}{0ex}}\mathbf{C}\in {\mathbb{R}}^{N×N}$, and $\mathbf{K}\in {\mathbb{R}}^{N×N}$ are the inertia, damping and stiffness matrices, and N is the number of generalized coordinates. The generalized coordinates and their first and second time derivatives (velocities and accelerations) are shown as u, $\stackrel{\mathrm{˙}}{\mathbit{u}}$ and $\stackrel{\mathrm{¨}}{\mathbit{u}}$. Lagrange multipliers and constraint equations are represented by $\mathbit{\lambda }\in {\mathbb{R}}^{{N}_{\mathrm{c}}}$ and $\mathbit{g}\in {\mathbb{R}}^{{N}_{\mathrm{c}}}$, where Nc is the number of constraints in the model. The Jacobian of constraint equations with respect to the generalized coordinates is presented by ${\mathbf{G}}_{\mathrm{u}}\in {\mathbb{R}}^{{N}_{\mathrm{c}}×N}$. Generalized external forces and quadratic velocity vectors, including gyroscopic and Coriolis force components, are shown as f and fv. The solver computes u, $\stackrel{\mathrm{˙}}{\mathbit{u}}$, $\stackrel{\mathrm{¨}}{\mathbit{u}}$ and λ at each time step for known external loads while satisfying the constraint equations. In HAWC2, the computed structural response (u, $\stackrel{\mathrm{˙}}{\mathbit{u}}$, $\stackrel{\mathrm{¨}}{\mathbit{u}}$) is sent to the aerodynamic solver. Based on these state variables, the aerodynamic solver computes the corresponding aerodynamic loads which go into the external force vector (f). This load update procedure takes place at each iteration. Hence, the generalized external forces and inertia matrix are a function of time, deflections, velocities and accelerations.

$\begin{array}{}\text{(1)}& \begin{array}{rl}& \mathbf{M}\left(\mathbit{u}\right)\stackrel{\mathrm{¨}}{\mathbit{u}}\left({t}_{i}\right)+\mathbf{C}\stackrel{\mathrm{˙}}{\mathbit{u}}\left({t}_{i}\right)+\mathbf{K}\mathbit{u}\left({t}_{i}\right)+{\mathbf{G}}_{u}^{T}\left(u,{t}_{i}\right)\mathbit{\lambda }\left({t}_{i}\right)\\ & =\mathbit{f}\left(\mathbit{u},\stackrel{\mathrm{˙}}{\mathbit{u}},{t}_{i}\right)+{\mathbit{f}}_{\mathrm{v}}\left(\mathbit{u},\stackrel{\mathrm{˙}}{\mathbit{u}},{t}_{i}\right)\\ \mathbit{g}\left(u,{t}_{i}\right)=\mathbf{0},& {\mathbit{G}}_{u}\left(u,{t}_{i}\right)=\frac{\partial \mathbit{g}\left(u,{t}_{i}\right)}{\partial \mathbit{u}\left({t}_{i}\right)}\end{array}\end{array}$

As the reference–rigid body (ur) and elastic parts (ue) of the generalized coordinates are separated, Eq. (1) can be written as shown in Eq. (2) for body “j”. The stiffness and damping matrices of the body have only elastic components which are constant for linear elements. Similarly, Mee is also constant and the constant matrices are computed once in a FRF solution process. The rest of the M matrix needs to be computed at each iteration together with g, Gu, f and fv since they are state dependent.

$\begin{array}{}\text{(2)}& \begin{array}{rl}& \left[\begin{array}{ll}{\mathbf{M}}_{\mathrm{rr}}^{j}& {\mathbf{M}}_{\mathrm{re}}^{j}\\ {\mathbf{M}}_{\mathrm{er}}^{j}& {\mathbf{M}}_{\mathrm{ee}}^{j}\end{array}\right]\left[\begin{array}{l}{\stackrel{\mathrm{¨}}{\mathbit{u}}}_{\mathrm{r}}^{j}\\ {\stackrel{\mathrm{¨}}{\mathbit{u}}}_{\mathrm{e}}^{j}\end{array}\right]+\left[\begin{array}{ll}\mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathbf{C}}_{\mathrm{ee}}^{j}\end{array}\right]\left[\begin{array}{l}{\stackrel{\mathrm{˙}}{\mathbit{u}}}_{\mathrm{r}}^{j}\\ {\stackrel{\mathrm{˙}}{\mathbit{u}}}_{\mathrm{e}}^{j}\end{array}\right]\\ & +\left[\begin{array}{ll}\mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathbf{K}}_{\mathrm{ee}}^{j}\end{array}\right]\left[\begin{array}{l}{\mathbit{u}}_{\mathrm{r}}^{j}\\ {\mathbit{u}}_{\mathrm{e}}^{j}\end{array}\right]+\left[\begin{array}{l}{\mathbf{G}}_{{u}_{\mathrm{r}}}^{{j}^{T}}\\ {\mathbf{G}}_{{u}_{\mathrm{e}}}^{{j}^{T}}\end{array}\right]\left[\begin{array}{l}{\mathbit{\lambda }}_{\mathrm{r}}^{j}\\ {\mathbit{\lambda }}_{\mathrm{e}}^{j}\end{array}\right]\\ & =\left[\begin{array}{l}{\mathbit{f}}_{\mathrm{er}}^{j}\\ {\mathbit{f}}_{\mathrm{ee}}^{j}\end{array}\right]+\left[\begin{array}{l}{\mathbit{f}}_{\mathrm{vr}}^{j}\\ {\mathbit{f}}_{\mathrm{ve}}^{j}\end{array}\right]\end{array}\end{array}$

The main driving factors in computation time of the multibody solver are the simulation time, the size of the problem (matrices) and the number of iterations. The vector ${\mathbit{u}}_{\mathrm{r}}^{j}$ includes six variables to define the position and rotation of the body “j” reference point. The size of ${\mathbit{u}}_{\mathrm{e}}^{j}$ depends on the number of elements in body “j”. As more bodies are defined in a model, the number of generalized coordinates and state-dependent parts of the matrices increase. For example, the one-body case in Fig. 1 has 60 generalized coordinates (six reference coordinates, 54 elastic coordinates), whereas the three-body model has 72 generalized coordinates (18 reference coordinates, 54 elastic coordinates).

In HAWC2 the time integration is performed using the Newmark algorithm (Newmark1959) with β and γ constants. The updates of the current state are performed by u and λ, computed according to Eq. (3). In Eq. (3) rq and rg are the force and constraint residuals at the current iteration step. Keff is the effective tangent stiffness at the current state, which is shown in Eq. (4). The sparsity of the constraint Jacobian matrix (Gu) increases with the number of constraints defined in the model. Different numerical approaches can be used when solving dense or sparse matrix problems. HAWC2 can optionally utilize a sparse matrix solution method in which λ from Eq. (3) is computed using the PARDISO sparse matrix routines . Note that $\left({\mathbf{G}}_{\mathrm{u}}{\mathbf{K}}_{\mathrm{eff}}^{-\mathrm{1}}{\mathbf{G}}_{\mathrm{u}}^{T}\right)$ is symmetric and positive definitive for the considered HAWC2 models.

$\begin{array}{}\text{(3)}& \begin{array}{rl}△\mathbit{\lambda }& =\left({\mathbf{G}}_{\mathrm{u}}{\mathbf{K}}_{\mathrm{eff}}^{-\mathrm{1}}{\mathbf{G}}_{\mathrm{u}}^{T}{\right)}^{-\mathrm{1}}\left({\mathbf{G}}_{\mathrm{u}}{\mathbf{K}}_{\mathrm{eff}}^{-\mathrm{1}}△{\mathbit{r}}_{\mathrm{q}}-△{\mathbit{r}}_{\mathrm{g}}\right)\\ & △\mathbit{u}={\mathbf{K}}_{\mathrm{eff}}^{-\mathrm{1}}\left(△{\mathbit{r}}_{\mathrm{q}}-{\mathbf{G}}_{\mathrm{u}}^{T}△\mathbit{\lambda }\right)\end{array}\text{(4)}& {\mathbf{K}}_{\mathrm{eff}}=\frac{\mathrm{1}}{\mathit{\beta }{h}^{\mathrm{2}}}\mathbf{M}+\frac{\mathit{\gamma }}{\mathit{\beta }h}\mathbf{C}+\mathbf{K}\end{array}$

## 2.2 Analysis

The approach in the study is based on numerical experiments of two turbines: the DTU10MW and the IEA10MW . The corresponding HAWC2 input files used for this study can be found in . However, the versions of the DTU and IEA10MW models used here are slightly different compared to the original models: the DTU10MW has a small offset on the blade twist distribution, and the IEA10MW has a different drivetrain mass and inertia.

The properties of the blade models are shown in Table 1. It should be noted that the IEA10MW rotor has more prebend and lower blade frequencies (see Table 1), which implies a more flexible blade structure and larger geometrical couplings when compared to the DTU10MW. These differences are relevant when considering the nonlinear geometrical response of a wind turbine rotor.

Table 1General properties of the reference wind turbines: DTU10MW and IEA10MW.

It is practical to call the bodies used for a continuous structure or a component the main body, and the bodies defined in a main body are called sub-bodies. A main body can be attached to other bodies or boundaries by constraints in any direction, whereas the constraints between the sub-bodies are always in 6 dof to satisfy the continuity of the structure. In the analyses the number of sub-bodies of the blade varied from 1 (linear response) to 30 (one body for each element, equivalent to a corotational approach). The rest of the turbine model was kept the same for a coherent comparison. The HAWC2 models of the considered turbines for this publication are composed of nine main bodies: tower, tower top, nacelle, three hubs and three blades. Table 2 shows the number of sub-bodies in the turbine models and the number of beam elements in each body. The tower, tower top, nacelle and hubs are modeled via one sub-body; in other words they are modeled as linear structures. Blades are the only parts which are modeled by multiple sub-bodies to capture large deflections. Both turbine models have 50 aerodynamic sections (or calculation points) on each blade, and the open-source Basic DTU Wind Energy controller was used. The turbulence boxes were generated by the Mann turbulence generator (Mann1994). A constant time step of 0.01 s was used for all considered cases. The computational time was recorded for all cases, and both the sparse and dense matrix solvers were considered.

Table 2HAWC2 turbine models' main bodies and number of elements and sub-bodies used in each main body.

The number of bodies in the model alters the problem size since it changes the number of generalized coordinates and constraints in the equations. The number of generalized coordinates and constraint equations can be determined by Eqs. (5)–(6). In the equations, Nmb is the number of main bodies, and ${N}_{\mathrm{el}}^{i}$ and ${N}_{\mathrm{sb}}^{i}$ are the number of elements and sub-bodies in the ith main body. The number of bodies in each blade model varies from 1 to 30. The 30-sub-body blade model (similar to the corotational model) is the most accurate with the highest N and Nc, whereas the one sub-body blade case is the linear blade case. Table 2 shows the element numbers at each main body in the turbine models. In all cases, the blades dominate the problem sizes. For example, in the one-body case the blades have 558 generalized coordinates and 18 constraint equations. For the 30-sub-body case, the three blades have 1080 generalized coordinates and 540 constraint equations. Although the problem size in the FRF formulation changes with the number of bodies defined in the model, the number of independent coordinates (NNc) is always 648 for this turbine model.

$\begin{array}{}\text{(5)}& N=\sum _{i=\mathrm{1}}^{{N}_{\mathrm{mb}}}\left({N}_{\mathrm{el}}^{i}+{N}_{\mathrm{sb}}^{i}\right)×\mathrm{6}\text{(6)}& {N}_{\mathrm{c}}=\sum _{i=\mathrm{1}}^{{N}_{\mathrm{mb}}}{N}_{\mathrm{sb}}^{i}×\mathrm{6}\end{array}$

Table 3Number of generalized coordinates N and constraint equations Nc for the full turbine model. The number of sub-bodies refers to the sub-bodies for the different blade models; it does not refer to the total number of sub-bodies of the entire turbine.

Table 4Qualitative breakdown of the fatigue load contributions of various design load cases for the IEA10MW (based on the loads reported in ). Hour distribution is based on a 20-year lifetime.

Table 5Design load cases (DLC) 1.2 power production on normal turbulence load case simulation setup.

The stability analysis includes a turbine model which is free to speed up without generator torque and has a fixed blade pitch angle at 0 . A steady-state rotor speed under zero aerodynamic torque is found close the cut-in wind speed. From there, the wind speed is increased following a shallow linear ramp. Consequently, the rotor slowly accelerates. The instability is then determined when significant blade vibrations are observed.

3 Results

The simulation results of the blade models with different numbers of sub-bodies are compared to the blade with the 30-body case (highest fidelity). The loads and total number of iterations are normalized with respect to the highest-fidelity results, while the computation time is normalized with respect to the lowest-fidelity model (one sub-body, linear case) in combination with the dense matrix solver. The computation time and total iteration number are defined here as the total central processing unit (CPU) time and the sum of iterations for all load cases, respectively.

The activity of the pitch bearing is evaluated by integrating the pitch angle signal over time for all load cases; see Eq. (7). The pitch angular speed of the jth blade at the ith time step is shown by ${\stackrel{\mathrm{˙}}{\mathit{\varphi }}}_{i}^{j}$. There are Nt number of time steps in all load cases. In addition to the total pitch angle change ϕtotal, the power needed by the pitch actuator (${P}_{i}^{j}$) of the jth blade at the ith time step is calculated by considering the torsion moment at the blade root (${M}_{i}^{j}$) and angular speed of the pitch bearing ${\stackrel{\mathrm{˙}}{\mathit{\varphi }}}_{i}^{j}$; see Eq. (8). Note that the bearing friction is neglected in the equation. The max power needed by the pitch actuator might determine the size of the component (i.e., actuator, bearing).

$\begin{array}{}\text{(7)}& {\mathit{\varphi }}_{\mathrm{total}}=\sum _{j=\mathrm{1}}^{\mathrm{3}}\sum _{i=\mathrm{1}}^{{N}_{t}}\frac{{\stackrel{\mathrm{˙}}{\mathit{\varphi }}}_{i-\mathrm{1}}^{j}+{\stackrel{\mathrm{˙}}{\mathit{\varphi }}}_{i}^{j}}{\mathrm{2}}\mathrm{\Delta }{t}_{i}\text{(8)}& {P}_{i}^{j}={M}_{i}^{j}×{\stackrel{\mathrm{˙}}{\mathit{\varphi }}}_{i}^{j}\phantom{\rule{0.25em}{0ex}}\text{at}\phantom{\rule{0.25em}{0ex}}i\text{th time step}\end{array}$

Figure 2 shows DLC 1.2 load cases' computation time and number of total iteration ratios of both turbines for dense and sparse matrix solvers. The computation time ratio is calculated with respect to the linear (one sub-body) case using the dense solver, and the ratio of the number of iterations is calculated with respect to the 30-sub-body blade case, which has the lowest number of iterations for both turbine models. The total number of iterations does not change for sparse and dense matrix solver types; therefore there is only one curve for the number of iterations. The dense matrix solver CPU time results are given only for 1-, 2-, 6-, 15- and 30-sub-body cases. The computation time is dependent on the number of iterations observed in a simulation and the number of sub-bodies of the blade. Therefore, it is possible to observe a decrease in computation time as the number of dof's and constraint equations increases. The number of iterations decreases until the 15-sub-body case, which also affects the CPU time accordingly. After the 15-sub-body case, the number of iterations remains approximately constant, and correspondingly the CPU time increases as the number of bodies increases.

Figure 2The total number of iterations is normalized by the result of the 30-sub-body case, and the total CPU time is normalized by the one sub-body dense matrix solver case for the DTU10MW and IEA10MW turbines. The CPU time ratios are given for both dense and sparse solvers.

The maximum dense solver computation time is observed for the 30-sub-body case. It is approximately 62 % and 70 % (see Fig. 2) slower compared to the linear case for DTU10MW and IEA10MW. Due to a sharp reduction in the number of iterations between the one- and three-sub-body cases, the computational time decreases as well, even though the complexity of the model increases. Hence, the dense solver computational cost due to the increase in model complexity increases more slowly compared to the time gained by having fewer iterations. The number of iterations decreases only moderately between the three- and 15-sub-body cases, which is followed by a modest increase in computational time. It is only after the 15-sub-body cases, for which the total number of iterations is roughly constant, that a continuous increase in the computational time is observed as function of the number of sub-bodies. It is further interesting to note that there is no significant difference in terms of computational cost between the one- and 15-sub-body cases due to the fact that approximately 36 % and 41 % fewer iterations were observed for DTU10MW and IEA10MW, respectively.

Since the sparsity of the matrices in Eq. (3) increases with the number of bodies, the sparse matrix solver becomes computationally more efficient for models with many constraints or bodies . Although not shown here, no difference was observed between the results of the dense and sparse matrix solvers. For the linear case, the CPU time is almost the same for both solver types. The sparse solver is significantly faster for the nonlinear (multibody) cases. The computational speedup for the 15-sub-body case is about 11 % and it is actually faster than the linear case with the dense matrix solver. Obtaining the sparse solution of 30-sub-body cases for the IEA turbine is about 36 % faster than using dense matrix techniques. The highest-fidelity model with a sparse matrix solver is just 9 % slower than the linear case for the IEA turbine, and this number goes down to 4 % for DTU turbine.

Figure 3Linear and nonlinear blade model power, blade pitch (left axis in the figures) and effective blade radius change (right axis in the figures) results with respect to wind speeds for steady-wind-load cases. Panel (a) shows the DTU10MW results and IEA results are given in (b).

Figure 4 shows blade torsion deformation results at 75 % blade span and blade tip for the linear and nonlinear blade models. Since the IEA blade is more flexible and longer than the DTU blade, the IEA torsional deflections are up to 1 larger than DTU deflections. The blade torsion deflections are large enough (up to 2.4 at IEA blade tip) to alter the turbine loads and performance. The deflections become significant in particular after the rated wind speed where the pitch activity is high. Although the torsional deformation curves of the linear and nonlinear models with respect to wind speeds look similar for the DTU10MW, the IEA10MW blade deformation curves for the linear and nonlinear models look quite different after the rated wind speed.

Figure 4Linear and nonlinear blade model torsion deformations at 75 % blade span and blade tip with respect to wind speeds for steady-wind-load cases. Panel (a) shows the DTU10MW results and IEA results are given in (b).

## 3.2 DLC 1.2 blade results

Figure 5Normalized blade tip minimum tower clearance, maximum effective blade radius and blade edgewise deflection results. Values are normalized with respect to the results of the 30-sub-body case.

Figure 6 shows the lifetime damage equivalent load (DEL) ratios between the linear (one sub-body) and nonlinear (30 sub-bodies) blade models over the normalized blade span for the DTU10MW and the IEA10MW turbines. The IEA turbine has a larger difference between linear and nonlinear cases in edgewise and flapwise DEL moments than the DTU turbine, but not so for the torsion DEL. A significant difference between the linear and nonlinear cases (30 sub-bodies) of more than 20 % can be observed for certain outboard radial stations. The flap- and edgewise DELs are consistently overestimated for the linear case, while the torsion DEL is underestimated with respect to the 30-sub-body nonlinear case.

Figure 6Flapwise, edgewise and torsion moment DEL ratios between linear (one sub-body) and nonlinear (30 sub-bodies) blade model over the normalized blade span for the DTU10MW and IEA10MW turbines.

Figure 7Normalized flapwise, edgewise and torsion moment DEL ratio variations with respect to the number of blade model sub-bodies at blade stations where the maximum deviations between linear and nonlinear cases occur.

Figure 7 shows flapwise, edgewise and torsion moment DEL ratio variations by model fidelity (number of sub-bodies in blade model) at blade stations where the maximum deviations between linear and nonlinear cases occur for each load component. The results are normalized with respect to the highest-fidelity blade model. The maximum deviation of the IEA turbine in flapwise deviation is 24 %, and it increases to 26 % for the edgewise direction. The DTU turbine has 9 % and 5 % deviations in the flapwise and edgewise directions. The results of both turbines in the flapwise and edgewise directions have a similar trend, meaning that after 15 sub-bodies the deviations become very small. The torsion DEL has the largest deviations for both turbines, and only after the nine-sub-body case can a consistent reduction in difference between the linear and nonlinear cases be observed. The deviations become quite small for cases with 27 sub-bodies or more.

Figure 8 shows the absolute maximum moment load result ratio between linear and nonlinear blade models over the normalized blade spanwise locations. The IEA results generally have larger deviations than the DTU results. The largest difference occurs in torsion moments for both turbines. The difference in the flapwise direction reaches up to 30 % for the IEA turbine and 10 % for the DTU turbine. The edgewise deviations of both turbines reach up to 12 %. The torsion moment deviation hits 50 % in some blade regions for the IEA turbine. The torsion moments are underestimated by linear models, whereas the flapwise and edgewise moments are generally overestimated by linear models.

Figure 8Linear and nonlinear blade absolute maximum moment load result ratio variation in DTU10MW and IEA10MW turbines with respect to blade span location.

Figure 9Cross-section flapwise and edgewise load envelopes at 43.6 m blade radius of the DTU turbine and 51.1 m blade radius of the IEA turbine for 1-, 2-, 6-, 15- and 30-sub-body cases.

## 3.3 DLC 1.2 tower and performance results

Turbine tower damage equivalent load, maximum and cross-section load results are given for the DLC 1.2 load cases. The turbine performance results are also mentioned here. Figure 10 shows the normalized maximum tower-top (yaw-bearing) torsion moments and maximum tower-top accelerations. In the case of excessive tower-top accelerations, the controller starts an emergency stop procedure. The difference in yaw-bearing torsion moment can reach up to 10 % for the IEA turbine. The results approach the highest-fidelity results very fast, and after the nine-sub-body case the deviations become very small compared to the 30-sub-body cases. The difference in tower-top accelerations can be more than 4 % between the linear and nonlinear cases.

Figure 10Normalized tower-top maximum torsion moments, side–side and fore–aft accelerations. Values are normalized with respect to the results of the 30-sub-body case.

Figure 11 shows the DELs of the fore–aft (moment force vector perpendicular to wind direction) and side–side (moment force vector aligned with the wind direction) moments at the tower-top position where the yaw actuator and bearings are located. There is a negligible deviation between the linear and nonlinear case for the side–side DEL moments for both turbine models. However, the deviations in fore–aft and torsion DELs exceed 4 % for the IEA turbine and reach 3 % for the DTU turbine. The results approach the highest fidelity model results smoothly and the deviation becomes very small after 15-sub-body cases for all channels. Figure 12 shows the tower bottom side–side and fore–aft moment load envelopes of the turbines for 1-, 2-, 6-, 15- and 30-sub-body cases. The deviations between linear and nonlinear cases are more explicit in the IEA10MW turbine than the DTU10MW turbine. In contrast to the blade moment envelopes, the linear case is not always the more conservative approach compared to the nonlinear cases.

Figure 11Normalized tower-top torsion, side–side and fore–aft DEL moment with respect to number of sub-bodies in blade model.

Figure 12Tower bottom fore–aft and side–side load envelopes of the DTU and the IEA turbines for 1-, 2-, 6-, 15- and 30-sub-body cases.

Figure 13 shows the normalized blade pitch actuator DEL, total pitch angle change of the turbines in all simulations computed by Eq. (7) and maximum power at pitch actuator computed via Eq. (8). The IEA turbine has a deviation of about 3 % in cumulative pitch angle results. This indicates that the controller activity is also affected by the fidelity of blade modeling. The maximum pitch actuator power depends on both blade root torsion moment and pitch angle speed. A very large deviation is observed in the pitch power results, which are 38 % and 34 % for the DTU and IEA turbines, respectively. The deviations in the IEA turbine results are generally higher than the DTU10MW turbine results; however the DTU turbine has larger deviations in terms of percentage than the IEA turbine in pitch power results. Although the highest-fidelity model causes slightly less pitch activity compared to the linear model, the actuator power increases significantly with the fidelity of the blade model. This is due to significantly increased blade torsional moments with increasing blade model fidelity.

Figure 13Normalized blade pitch actuator (blade root torsion moment) DEL, total pitch angle change for all load cases and maximum power at pitch actuator with respect to number of blade sub-bodies.

The difference in annual energy production (AEP) between the different blade models is well below 1.0 %. This difference is relatively small when compared to the loads since the controller tracks the optimal operating conditions below rated wind speed and maintains the rated power above rated wind speed. Consequently, only in below rated conditions can a very small difference in power output be observed whereby the linear case results in small increase in power output compared to the nonlinear 30-sub body case.

## 3.4 Stability results

The stability of DTU and IEA turbines is evaluated by considering the linear (one body) and nonlinear (30 bodies) blade models. Blade tip torsion deformation depicts the blade vibrations and instability (flutter) clearly. Figure 14 shows the rotor speed and blade tip torsion results with respect to the wind speed. Turbines have zero aerodynamic torque at the initial wind speed, and wind speed acceleration is 0.0145 m s−2. Results show that the DTU turbine has much higher flutter speeds than the IEA turbine for both blade models. The DTU blade linear model (blue curve) shows the flutter instability at almost the same rpm's with the nonlinear model and 1 m s−1 higher wind speed compared to the nonlinear model (red curve). However, for the IEA linear model, flutter occurs at a wind speed which is 3.6 m s−1 lower compared to the nonlinear model. Furthermore, the rotational speed difference between the linear and nonlinear models for the flutter instability is more than 8 rpm for the IEA blade. This shows that the linear models do not always overestimate the flutter speeds.

Figure 14DTU10MW (a) and IEA (b) rotor speed and blade tip torsion deformation results with respect to wind speed. The flutter wind speeds are shown in the figures.

4 Discussion and conclusion

The effects of blade structural model fidelity on the turbine response, loads, stability and computation time are investigated in this study. The blades are modeled by different numbers of sub-bodies in the multibody formulation of HAWC2. The blade model geometric nonlinearity is changed from linear to the highest available fidelity level, which is equivalent to a corotational formulation. The effects of blade geometric nonlinearities are compared by exploring the results of two different blade designs with otherwise identical tower and shaft configuration. The normal power production load cases are selected according to the IEC 61400-1 standard (DLC1.2) but considering 12 instead of six turbulent seeds. In addition, the computational speed of the dense and sparse matrix solvers as used by HAWC2 is compared for different blade model fidelities.

Although there are significant differences between the linear and the nonlinear blade models (with 30 sub-bodies), the results generally approach the highest-fidelity results fast as the number of blade sub-bodies increases. In most of the studied cases the deviations in results become insignificant after 15 sub-bodies. This is also the point after which the total number of iterations does not decrease any further significantly with increasing number of sub-bodies.

The work outlined here confirms earlier studies that the nonlinear geometrical effects are significant for wind turbine blades, even more so for new turbine designs (DTU10MW vs. IEA10MW). The geometrically nonlinear effects are model dependent and are related to the size, prebend shape and flexibility of the considered blade model. The authors conclude that users are recommended to model blades with as many sub-bodies as there are structural elements, while also using a sparse matrix solver for models that have symmetric effective stiffness matrices in HAWC2. In doing so within the context of HAWC2, no increase in CPU time is noted while at the same time having the blade model with the highest structural fidelity.

Data availability
Data availability.

The corresponding HAWC2 input files used for this study can be found in in addition to all the statistical properties of all channels and time series that have been used to generate the figures in this publication. Due to the large number of simulations, and the corresponding data storage requirements, the time series are not made available for redistribution. However, these can be reproduced using the given input files in combination with HAWC2 version 12.6.

Author contributions
Author contributions.

OG conducted the study as part of his PhD research. The idea was developed by OG and DRV. OG ran the analyses and performed the post-processing including results and figures with tools developed by DRV. The manuscript was written jointly by OG and DRV.

Competing interests
Competing interests.

DTU Wind Energy develops, supports and distributes HAWC2 on commercial terms.

Acknowledgements
Acknowledgements.

The authors are grateful to Sergio González Horcas and Anders Melchior Hansen for their contribution to the HAWC2 solver. Furthermore, the authors would like to acknowledge their colleagues, Sergio González Horcas, Mathias Stolpe, Suguang Dou, Riccardo Riva and Georg Pirrung for their comments on the manuscript.

Financial support
Financial support.

This study was funded by DTU Wind Energy.

Review statement
Review statement.

This paper was edited by Lars Pilgaard Mikkelsen and reviewed by two anonymous referees.

References

DVN GL: BLADED Theory Manual version 4.9, Tech. rep., Garrad Hassan and Partners Ltd, Bristol, UK, 2018. a

Bak, C., Zahle, F., Bitsche, R., Kim, T., Yde, A., Henriksen, L. C., Hansen, M. H., and Natarajan, A.: The DTU 10-MW Reference Wind Turbine, Tech. Rep. Report-I-0092, DTU Wind Energy, 2013. a, b

Beardsell, A., Collier, W., and Han, T.: Effect of linear and non-linear blade modelling techniques on simulated fatigue and extreme loads using Bladed, J. Phys. Conf. Ser., 753, 042002, https://doi.org/10.1088/1742-6596/753/4/042002, 2016. a

Bortolotti, P., Tarres, H. C., Dykes, K. L., Merz, K., Sethuraman, L., Verelst, D., and Zahle, F.: IEA Wind TCP Task 37: Systems Engineering in Wind Energy – WP2.1 Reference Wind Turbines, Tech. Rep. NREL/TP-5000-73492, National Renewable Energy Lab. (NREL), Golden, CO (United States), https://doi.org/10.2172/1529216, 2019. a, b, c

Cardona, A. and Géradin, M.: Flexible multibody dynamics: a finite element approach, John Wiley, Chichester, UK, 2001. a

De Vries, O.: Fluid dynamic aspects of wind energy conversion, Tech. Rep. AGARDograph No. 243, AGARD Advisory Group for Aerospace Research and Development, Amsterdam, The Netherlands, 1979. a

Dibold, M., Gerstmayr, J., and Irschik, H.: On the accuracy and computational costs of the absolute nodal coordinate and the floating frame of reference formulation in deformable multibody systems, in: Proceedings of the ASME 2007 International Design Engineering Technical Conferences & Computers and Information in Engineering Conference, vol. 34756, 1–10, Las Vegas, Nevada, USA, 2007. a

Gözcü, O. and Verelst, D. R.: HAWC2 input data and statistics of time series for the DTU10MW and IEA10MW wind turbines for DLC1.2, https://doi.org/10.11583/DTU.9771722, 2020. a, b

Guntur, S., Jonkman, J., Sievers, R., Sprague, M. A., Schreck, S., and Wang, Q.: A validation and code-to-code verification of FAST for a megawatt-scale wind turbine with aeroelastically tailored blades, Wind Energ. Sci., 2, 443–468, https://doi.org/10.5194/wes-2-443-2017, 2017. a

Hansen, M. H. and Henriksen, L. C.: Basic DTU Wind Energy Controller, Tech. Rep. No. 0028, DTU Wind Energy, 2013. a

Hansen, M. H., Thomsen, K., Natarajan, A., and Barlas, A.: Design Load Basis for Onshore Turbines, DTU Wind Energy, No. 0174, 2015. a

Hansen, M. O. L., Sørensen, J. N., Voutsinas, S., Sørensen, N., and Madsen, H. A.: State of the Art in Wind Turbine Aerodynamics and Aeroelasticity, Prog. Aerosp. Sci., 42, 285–330, 2006. a

Hodges, D. H.: A mixed variational formulation based on exact intrinsic equations for dynamics of moving beams, Int. J. Solids Struct., 26, 1253–1273, 1990. a

IEC: 2005 IEC 61400-1 3rd edition Wind turbines – Part 1: Design requirements, International Electrotechnical Commission, IEC, 2005. a

Jonkman, J., Butterfield, S., Musial, W., and Scott, G.: Definition of a 5-MW Reference Wind Turbine for Offshore System Development, National Renewable Energy Laboratory, Technical Report No. NREL/TP-500-38060, 2009. a

Jonkman, J. M. and Buhl Jr., M. L.: FAST User's Guide, Tech. Rep. NREL/EL-500-38230, National Renewable Energy Laboratory (NREL), 2005. a

Kallesøe, B. S.: Effect of steady deflections on the aeroelastic stability of a turbine blade, Wind Energy, 14, 209–224, https://doi.org/10.1002/we.413, 2011. a

Kim, T., Hansen, A. M., and Branner, K.: Development of an anisotropic beam finite element for composite wind turbine blades in multibody system, Renew. Energ., 59, 172–183, 2013. a

Krenk, S.: Nonlinear Modelling and Analysis of Structures and Solids, Cambridge University Press, Cambridge, 2005. a

Larsen, T., Hansen, A., and Buhl, T.: Aeroelastic effects of large blade deflections for wind turbines, in: TORQUE 2010: The science of making torque from wind, Delft, The Netherlands, 238–246, 2004. a

Larsen, T. J. and Hansen, A. M.: How 2 HAWC2 the user's manual, Tech. Rep. Risø-R-1597(ver. 12.8)(EN), DTU Wind Energy, Roskilde, Denmark, 2019. a

Madsen, H. A., Larsen, T. J., Pirrung, G. R., Li, A., and Zahle, F.: Implementation of the blade element momentum model on a polar grid and its aeroelastic load impact, Wind Energy Science, 5, 1–27, https://doi.org/10.5194/wes-5-1-2020, 2020. a

Mann, J.: The spatial structure of neutral atmospheric surface-layer turbulence, J. Fluid Mech., 273, 141–168, https://doi.org/10.1017/S0022112094001886, 1994. a

Manolas, D. I., Riziotis, V. A., and Voutsinas, S. G.: Assessing the Importance of Geometric Nonlinear Effects in the Prediction of Wind Turbine Blade Loads, J. Computat. Nonlin. Dyn., 10, 041008, https://doi.org/10.1115/1.4027684, 2015. a

Natarajan, A. and Verelst, D. R.: Outlier robustness for wind turbine extrapolated extreme loads, Wind Energy, 15, 679–697, https://doi.org/10.1002/we.497, 2012. a

Newmark, N. M.: A method of computation for structural dynamics, J. Eng. Mech. Div.-ASCE, 85, 67–94, 1959. a

Pavese, C., Wang, Q., Kim, T., Jonkman, J. M., and Sprague, M. A.: HAWC2 and BeamDyn: Comparison Between Beam Structural Models for Aero-Servo-Elastic Frameworks, in: Proceedings of the EWEA Annual Event and Exhibition 2015, Paris, France, 2015. a

Petra, C. G., Schenk, O., and Anitescu, M.: Real-Time Stochastic Optimization of Complex Energy Systems on High-Performance Computers, Comput. Sci. Eng., 16, 32–42, https://doi.org/10.1109/MCSE.2014.53, 2014.a. a

Petra, C. G., Schenk, O., Lubin, M., and Gäertner, K.: An Augmented Incomplete Factorization Approach for Computing the Schur Complement in Stochastic Optimization, SIAM J. Sci. Comput., 36, C139–C162, https://doi.org/10.1137/130908737, 2014b. a

Pirrung, G. R., Madsen, H. A., and Kim, T.: The influence of trailed vorticity on flutter speed estimations, J. Phys. Conf. Ser., 524, 1–11, https://doi.org/10.1088/1742-6596/524/1/012048, 2014. a

Rezaei, M., Zohoor, H., and Haddadpour, H.: Aeroelastic modeling and dynamic analysis of a wind turbine rotor by considering geometric nonlinearities, J. Sound Vib., 432, 653–679, 2018. a

Riziotis, V., Voutsinas, S., Politis, E., Chaviaropoulos, P., Hansen, A. M., Madsen, A., and Rasmussen, F.: Identification of structural nonlinearities due to large deflections on a 5 Mw wind turbine blade, in: Proceedings of the European Wind Energy Conference and Exhibition 2008, 102–112, Brussels, Belgium, 2008. a

Rubak, R. and Petersen, J. T.: Monopile as part of aeroelastic wind turbine simulation code, in: Proceedings of Copenhagen Offshore Wind, 1–10, 2005. a

Shabana, A. A.: Computational dynamics, John Wiley & Sons, Chichester, UK, 2010. a, b

Shabana, A. A.: Dynamics of multibody systems, Cambridge University Press, Cambridge, 2013. a, b

Tibaldi, C., Henriksen, L. C., and Bak, C.: Investigation of the dependency of wind turbine loads on the simulation time, in: Proceedings of the EWEA Annual Event and Exhibition 2014, Barcelona, Spain, 2014. a

Verelst, D. R., Hansen, M. H., and Pirrung, G.: Steady State Comparisons HAWC2 v12. 2 vs HAWCStab2 v2. 12, Tech. Rep. E-0122, DTU Wind Energy, 2016. a

Wang, Q., Sprague, M. A., Jonkman, J., Johnson, N., and Jonkman, B.: BeamDyn: a high-fidelity wind turbine blade solver in the FAST modular framework, Wind Energy, 20, 1439–1462, 2017. a

Wilson, R. E. and Lissaman, P. B.: Applied aerodynamics of wind power machines, Tech. Rep. PB-238595, Oregon State University, USA, 1974.  a

Zierath, J., Rachholz, R., Woernle, C., and Müller, A.: Load calculation on wind turbines: validation of Flex5, Alaska/Wind, MSC. Adams and Simpack by means of field tests, in: Proceedings of the ASME 2014 International Design Engineering Technical Conferences & Computers and Information in Engineering Conference, 2014–34670, 1–13, Buffalo, New York, USA, 2014. a