Interactive comment on “ Simulation of transient gusts on the NREL 5 MW wind turbine using CFD ” by Annika Länger-Möller

Simulation of transient gusts on the NREL5 MW wind turbine using CFD By Annika Länger-Möller The article describes a CFD approach to gust modeling for wind turbines, presenting some interesting results. As the CFD model is incompressible, any change in velocity in the domain is instantaneously influencing the full domain. As a consequence, the gust is not traveling through the domain, but the velocity is basically increasing everywhere in the domain instantaneously. The consequence of this with respect to the wake development should be discussed. Additionally, the assumption of infinite mass and inertia along with a non-elastic model C1


Introduction
In the past years, the growing computer power enabled the resolved simulation of wind turbines including the sites with Computational Fluid Dynamics (CFD).Studies have been performed for example by by Schulz et al. (Schulz et al., 2016) or Murali et al. (Murali and Rajagopalan, 2017) who used an Unsteady Reynolds-Averaged Navier Stokes (U-RANS) solver to perform according studies.Moreover, hybrid Large Eddy Simulation(LES)-RANS approaches are implemented to analyse the behaviour of wind turbines in a complex terrain as for example presented by Castellani et al. (Castellani et al., 2017) who also considered the unsteady atmospheric inflow conditions.The challenges of correctly predicting uncertainty of the fluc-20 tuating wind loads is a research field on its own.For example, Bierbooms et al. and Suomi et al. (Bierbooms and Drag, 1999;Suomi et al., 2013) investigated wind fields to better understand the shape of wind gusts.Matthäus et al. (Matthäus et al., 2017) argued that a detailed understand-ing of wind fields is not necessary.Matthäus et al. rather respected unknowns of all parts of the wind turbine life cycle as for example changes in the blade shape due to production tolerances, ageing, or the wind field and summarized them in uncertainty parameters to estimate the effective power outcome and rotor loads.Mücke et al. (Mücke et al., 2011) proved that turbulent wind fields are not distributed Gaussian as assumed in the International Electronic Committee Standard 61400-1 (IEC, 2005).A similar conclusion was drawn by Graf et al. (Graf et al., 2017) who investigated whether the 50 year loads as defined in the IEC 61400-1 standard (IEC, 2005) adequately fulfil their purpose by applying different approaches of probability prediction to the generic NREL 5 MW turbine using the FAST rotor code.The aerodynamic interferences between the unsteady wind conditions and wind turbines are of major importance for the prediction of fatigue loads and the annual power production.Therefore it is part of the certification computation for each wind turbine.Nevertheless, the detailed investigation of isolated effects as the 50 year Extreme Operating Gust (EOG) on the flow of a wind turbine using high fidelity methods as CFD is rare even though the blade loads resulting from the extreme load cases are dimensioning load cases.In the case of vertical axis wind turbines Scheurich et al. (Scheurich and Brown, 2013) analysed the power loss of a wind turbine subjected to a sinusoidal fluctuation in wind speed.However, amplitudes were small compared to the EOG.Horizon-tal axis wind turbines which are hit by an EOG as defined in the IEC 64100-1 standard (IEC, 2005) were presented by Uzol et al. in 2013 (Sezer-Uzol andUzol, 2013).The wind turbine under consideration was the NREL phase VI rotor with a wind speed of 7 m/s using the panel code AeroSIM+.
The impact of the gust was then evaluated regarding rotor thrust, torque and wake development.Preceeding this study, Bierbooms (Bierbooms, 2005) examined the flap moment of wind turbine blades which were subjected to a gust with extreme raise in 2005, also using CFD.

10
Even though the literature on gust simulations on wind turbines is rare, some research has been conducted in the field of aerospace science.Reimer et al. and Kelleners et al. (Kelleners and Heinrich, 2015;Reimer et al., 2015) presented two approaches: the velocity-disturbance approach and the resolved-gust approach to apply vertical gust on airplanes.The velocity-disturbance approach adds the gust velocity to the surface of the investigated geometry while the resolvedgust approach propagates the gust velocity through the flow field with the speed of sound.The validity of both implemen-20 tations was demonstrated by the time history of the position of the centre of gravity, pitch angles and load factors necessary for keeping the flight path of an aircraft constant.In the so-called field approach Parameswaran et al. (Parameswaran and Baeder, 1997) added the gust velocity to the grid velocity of the computational grid to all cells with wherein x is the coordinate in flow direction, u the transportation velocity and t the physical time.
The simulation of unsteady inflow conditions of wind turbines in CFD implies several challenges: The simulation of a wind turbine including the tower is, itself, an instationary problem which needs the computation of several rotations to obtain a periodic solution.Superposed by instationary (stochastic) inflow conditions, periodicity can never be gained because the same flow state never occurs twice.Moreover, a computation in which the rotor motion is adapted to the actual rotor forces using a strong coupling approach as proposed by Sobotta (Sobotta, 2015) should be included in the computation.By using strong coupling between CFD and a pitch control algorithm for the rotor-motion Sobotta has been able to implement a simulation-procedure of turbine start-up.Heinz et al. (Heinz et al., 2016) performed the computation of an emergency shut-down of a turbine also by using CFD and by neglecting the tower throughout the aerodynamic computations.Anyway, he considered the rotor mass and inertia by coupling of the CFD solver with the aeroelastic code HAWC2.The present study aims for filling the gap in the gust simulation of wind turbines by applying the resolved-gust approach to DLR's U-RANS solver THETA (Länger-Möller, 2017;Löwe et al., 2015).The generic NREL 5 MW wind turbine at rated operating conditions is chosen as test case for a 1−cos()-shaped gust which lasts about 7 s and the EOG following the international standard IEC 61400-1 definition.
Rotor mass and inertia are set to infinity.Thus a computation respecting a pitch controller is not necessary, as the rotor will not accelerate during the gust strike.As results, rotor thrust and torque, pressure distributions, friction force coefficients and wake vortex transport are evaluated.To further simplify the evaluation of aerodynamic effects due to the gust, the atmospheric boundary layer is not considered during the study.

Flow solver THETA
The DLR flow solver THETA is a finite volume method which solves the incompressible Navier-Stokes (NS) equation on unstructured grids.The grids can contain a mix of tetrahedrons, prisms, pyramids and hexahedrons.The transport equations are formulated on dual cells, which are constructed around each point of the primary grid.The transport equations are solved sequentially and implicitly.The Poisson equation which links velocity and pressure is either solved by the SIMPLE algorithm for stationary problems or the projection method for unsteady simulations.Pressure stabilization is used to avoid spurious oscillations caused by the collocated variable arrangement.
The technique of overlapping grids (Chimera) is used to couple fixed and moving grid blocks (Kessler and Löwe, 2012).The interpolation between the different blocks at interior boundaries is integrated in the system of linear equations on all grid levels of the multi-grid solver leading to an implicit formulation across the blocks which was identified to be crucial for achieving fast convergence of the Poisson equation.Implicit time-discretization schemes of first order (implicit Euler) or second order (Crank Nicolson, BDF) are implemented.A variety of schemes from first order upwind up to second order linear or quadratic upwind or a central scheme and the low dissipation, low dispersion second order scheme (Löwe et al., 2015) are implemented.Throughout this study, the second order central scheme is used.
The THETA code provides a user interface for setting complex initial and boundary conditions using the related C functions.All physical models are separated from the basis code.Therefore, new physical models can be implemented without modification of the base code.The user-interface is used extensively in the present study to prescribe the gust profile at the inflow boundary.
For turbulence modelling the commonly used Spalart-Allmaras, k − , k − ω or Menter-SST models are available.
In former studies by the author (Länger-Möller, 2017), the Menter-SST lead to most promising results for wind energy applications.Hence, this model is applied throughout the present study.The time step equivalent to a rotor advance of Ψ = 0.5 • is chosen together with the the Eulerian implicit scheme for the temporal discretization.

Grid characteristics
The computational grid consists of 3 parts.The first part contains the three rotor blades, stubs and the rotor hub.On the blade surface, a structured grid with 156 × 189 elements in spanwise and chordwise direction was generated.The boundary layer mesh of the blades consists of 49 hexahedra layers in an O-O-topology.The height of the wall-next cell is δ = 3 • 10 −6 m along the entire blade, ensuring y + ≤ 1.
The second part of the grid has the shape of a disk and contains the entire rotor.The disk measures D = 166.7 m in diameter and has a depth of 26.7 m.It is filled with tetrahedrons with an edge length between 0.002 m and 0.9 m.The entire disk is used as chimera child grid for the overlapping grid technique and contains approximately 11.63 • 10 6 points.The disk is displayed in figure 1(a) wherein the black coloured blade number 1 is at zero-azimuth position (Ψ = 0).
The chimera parent grid has the dimensions of 504 × 504 × 1512 m 3 in width, height, and length.The 54 prism layers, used to resolve the boundary layer of the viscous floor, have a total height of H = 5 m with a wall-next cell height of δ = 3 • 10 −5 m.The floor is defined as viscous wall.
In the chimera parent grid, the edge length of cells continuously grow from very small in the rotor-tower and wake region to rather large close to the farfield boundaries.The entire chimera parent grid contains approximately 13.25 • 10 6 points.
In figure 1 The procedure of applying the gust to the flow field starts by computing the flow field around the wind turbine in question until the flow field and the global rotor loads have become periodic.For the NREL 5MW turbine in the given setup 9 revolutions are required.Then, the inflow velocity on the in-60 flow boundary is modified according to the velocity change described in section 4.2 or 4.3.The computation is continued so that the gust is propagated through the flow field.In the approach by Reimer et al. and Kelleners et al. (Kelleners and Heinrich, 2015;Reimer et al., 2015) using CFD to solve the 65 compressible RANS equations, the gust front is transported with the speed of sound.In their approach as well as in the present paper the computation has been run at least until the gust has entirely passed the geometry in question but can be continued as long as wished by the user.

70
The restrictions on the resolved-gust approach named by Reimer et al. and Kelleners et al. (Kelleners and Heinrich, 2015;Reimer et al., 2015) to ensure a loss-free transport of the gust velocity are a fine grid upstream of the geometry in question

75
a fine time step.
As THETA is an incompressible solver, the speed of sound is infinite.In addition to the strong implicit formulation and the choice of boundary conditions that prevent the flow from escaping sideways, this leads to a spread of the gust velocity 80 through the flow field instantaneously.If the same gust velocity is added to the constant inflow condition in every point in the inflow plane and the boundary conditions are chosen as specified in section 3.2, the transport of the gust velocity will be loss-free through the entire domain.Hence, the 85 restrictions by Reimer et al. and Kelleners et al. (Kelleners and Heinrich, 2015;Reimer et al., 2015) regarding the grid resolutions are obsolete while a fine time step is required to ensure numerical stability.

Cosines gust
The 1 − cos() gust is modelled in analogy to the EASA certification standard (EASA, 2010) as with u(t) and u g as the time dependent velocity and the gust velocity, respectively and H as gust gradient.In the work presented, H is defined as (3) 10 wherein t represents the actual physical time and T S is the time at which the gust started.Inserting equation 3 in equation 2 the following definition of the gust results: Wherein T g is the duration time of the gust.The gust velocity u g is defined as +0.25 m/s representing a gust and −0.25 m/s representing a sudden calm.In both cases, the maximum change in wind speed is 0.5 m/s or 4.4 % at rated wind speed of the NREL 5MW turbine.The resulting inflow velocities are displayed in figure 2. (5) Wherein T g = 10.5 s is the gust characteristic time, t the computational time, u(z) the velocity profile depending on the height and the gust speed v g .The latter is defined as The resulting gust profiles of the EOG in comparison to the moderate 1 − cos()-gust is visualized in figure 2.

Constant inflow conditions
As described in section 4.1 a periodic flow field with periodic rotor loads is mandatory as initial conditions for computing gusts that act on wind turbines.The resulting history of rotor thrust F x and torque M x for constant inflow conditions over revolution 6 to 10 are displayed in figures 3 and 4 with the red line.In both figures the periodic behaviour of a periodic flow field is visible as well as the typical 3/rev-characteristic of a rotor-tower configuration of the wind turbine.The average value of F x and M x over the four revolutions displayed is 738.9N and 4.15 MNm, respectively.Compared to the reference (Jonkman et al., 2009), rotor thrust and torque at rated conditions, the values deviate about approximately −3.77 % and −0.98 %, respectively.Imiela et al. (Imiela et al., 2015) achieved a rotor thrust of 786 kN and a torque of 4.4 MNm in their studies with the compressible U-RANS solver TAU for the NREL 5MW turbine.The agreement between the CFD computation performed with THETA, TAU and the values obtained with linearised models which were used to generate the NREL 5MW documentation (Jonkman et al., 2009) is excellent.If

Cosines gust
The impact of the 1 − cos()gust is evaluated regarding rotor thrust F x and rotor torque M x during the gust by comparison to uniform inflow conditions.F x and M x are displayed in figure 3 and 4, respectively.Therein, the period between 30 s and 50 s is displayed while the gust operates between T S = 37 s and T S +T g = 44 s.Thus, it last approximately 1.5 rotor revolutions.As expected, the gust velocity spreads over the entire field immediately and also affects rotor thrust and torque instantaneously.In the case of gust and calm no hysteresis effect is found as rotor thrust and torque recover immediately after the gust.This is visible in both figures 3 and 4 as the curve of constant inflow conditions is matched right after 44 s.
During the gust, F x and M x follow the modification of the inflow condition.Hence, for a positive gust velocity u g (equation 4) rotor loads increase in a 1 − cos()-shape while they decrease in the same manner for a negative gust velocity.Additionally, the tower blockage effect is superposed on the blade loads and remains detectable in the blade load development.Furthermore, the tower blockage effect reduces the time the rotor experiences maximum loads in the case of u g = +0.25 m/s, as is visible at approximately t = 41 s.
A sharp drop in both rotor thrust and torque is visible.This Table 1 lists the relative differences in F x and M x during the gust, computed following equation 9.In equation 9 the sub-40 script max indicates the extreme rotor loads and the over-line the average of rotor loads under constant inflow conditions over 4 revolutions Additionally, the relative difference in the averaged blade load is computed by first integrating rotor thrust and torque during the gust excitation and then computing equation 9.
The result of the gust peak load is listed in table1 while the integrated loads are contained by table 2.

50
In both tables it can be seen that a reduction of the wind velocity due to a calm or the increase of the wind velocity with the same amplitude leads to almost identical absolute It is also important to note that the rotor loads return to the values of constant inflow conditions right after the gust ended.This indicates that there are no reflections or numerical oscillations in the flow field which lower the numerical accuracy.In summary, the behaviour of F x and M x is as expected.The increased wind velocity causes higher thrust and momentum and vice versa while the amplitude is identical for the increase and decrease in wind velocity. 10

Extreme operating gust
In figure 5 rotor thrust F x and torque M x are presented during the extreme gust excitation and constant inflow conditions.The extreme gust lasts about 10.5 s or 2 rotor revolutions.In comparison to the rotor loading during the gust, described in section 5.2, the tower blockage effect becomes negligible.Hence, the starting position of the rotor is less important.decreased by 36 % during the calm that precedes or follows the velocity peak and increased about 100 % during the gust peak.The changes in rotor thrust are smaller even though the amplitudes of load change are significant as well.Similar changes to those in rotor thrust and torque are, of course, also visible in pressure distributions and friction forces.To analyse the flow state on the blade during the gust two instances have been chosen: after t min = 2.5s + T S = 42 s (minimum gust velocity) and t max = 5.6s + T S = 45 s (maximum gust velocity).Figure 6a) and b),respectively, display the rotor positions in the instances investigated.In both figures blade number 1 is coloured in black.At t min , when the gust is at its minimum speed, blade number one is right in front of the tower and additionally experiences the tower blockage effect.Conversely, at t max , when the gust is at its maximum speed, blade number 1 is in freestream conditions while the flow on blade number 3 is influenced by the tower blockage.For a meaningful comparison with pressure and friction coefficients from constant inflow conditions, it was ensured that the investigated sections result from blades at the same azimuth positions.
The pressure distributions along with the friction forces are investigated at three radial sections: an inboard section at r/R = 10 %, a mid-section r/R = 50 % and an outboard section at r/R = 90 %, displayed in figure 7 to 9, respectively.In all three figures the pressure is displayed in the upper half and is normalized with the vector sum of the tip speed velocity and the constant inflow velocity.
A noticeable difference in c p at minimum gust velocity is found in all sections (figures 7 to 9) when compared to the constant inflow conditions.The decrease in the stagnation point of c p to lower values due to the calm at t min = 42 s is higher than that of the tower blockage effect while otherwise the pressure distributions keep the general shape.Conversely the difference between constant inflow conditions and the maximum gust is significantly higher.The maximum c p increased about 58 % in the blade tip section.Moreover, the shape of the pressure distribution is changed over the entire blade.This is visible especially in the mid-and outboard blade sections (figures 8 and 9).In both sections, the pressure increases rapidly in the rear half of the upper blade surface and even reaches positive values in the last 20 % of the profile.This behaviour is a first indication of a separation region and reversed flow around the trailing edge.
The indications about separation that were made in the pressure distributions are enhanced by analysing the friction coefficients on the blade.The friction coefficients on the blade sections in undisturbed flow are displayed in the lower half of figures 7 to 9. In all sections, strong fluctuations are visible at the trailing edge which result from the truncated geometry at x/c = 1.The friction force in the inboard section (figure 7) shows large differences between all considered time instances.The oscillations around x/c = 50 % at constant inflow conditions indicate a small separation region with otherwise attached flow.At minimum gust velocity the overall friction force level is increased around the leading edge Conversely to the inboard section, changes in separation at the mid-section appear due to the gust only.In figure 8 the friction force level is decreased at minimum gust velocity and the curve is very smooth.At maximum gust velocity, small oscillations appear around x/c = 90 %, indicating separation in that region.By increasing the rotor radius, the behaviour is enforced.Thus, at the outboard section (figure 9) the local maximum in c f in the last 20 % of the profile almost reaches the level of the leading edge during the maximum gust velocity.The same analysis of pressure coefficients is performed for the blade that is situated right in front of the tower or in the influence of the tower blockage 30 • behind the tower.The according figures are displayed in figures 10 to 12.For the pressure distributions (upper half) the same effects as described for the blades in undisturbed flow are found.The only difference is that due to the tower blockage the overall pressure level is decreased about 1, %. Conversely, the characteristics of friction forces (figures 10 to 12 (lower half)) in the inboard and outboard sections differ from those of the blades in undisturbed flow.In the inboard section in figure 10 Figure 13.Tip vortex transportation in horizontal plane through rotor centre the oscillations in c f reduce to a minimum while the overall friction force level remains constant.Only at maximum gust velocity, some fluctuations around the trailing edge appear.Thus, the tower seems to suppress separation and it takes some time until the separation state is fully recovered.The friction force coefficient in the mid-section behaves similar to the inboard section with the difference, that the separation region is increased (figure 11).In the outboard section in figure 12, the friction force is increased significantly at the maximum gust velocity and only a small separation region around the trailing edge is found.Finally, the transport of the tip vortices is investigated.In figure 13 and 14 three instances of the flow field are compared.The vortices are made visible with the λ 2 criterion (Jeong and Hussain, 1995).In both figures, the black lines represent the vortex transport at constant inflow, the blue lines are extracted at maximum gust velocity and the green ons at the end of the gust.By comparing the vortex transport at the beginning of the gust (black curve) and at the end of the gust (green curve), a compression and stretching of the vortex transport is found.This stretching relates closely to the time history of the inflow profile.The vortices that are shed the latest are positioned at x = −5 m; y = −60 m as well as at x = +5 m; y = +60 m in figure 13.Both vortices are transported identically at the beginning and end of the gust when the inflow velocity is identical.The two previously shed vortices around x = 13 m; y = −65 m and x = 27 m; y = −65 m were generated while the gust was at a lower inflow velocity (end of gust -green) or higher inflow velocity (maximum gust speed -blue).Thus they are not transported 30 as far (end of gust) or even further (maximum gust speed) and a slight compression and stretching towards the latest vortex

Conclusions
The study presented the simulation of the generic 5MW wind turbine operating under an extreme 50 years gust as defined in the IEC-64100-1 standard using the U-RANS CFD solver THETA.The gust has been introduced with a resolved-gustapproach by introducing the changing velocity at the inflow boundary condition.The gust velocity was then transported through the field with infinite speed of sound.The results obtained represented the effects that are expected during the in-stationary inflow condition in combination with the given boundary conditions.Rotor thrust and rotor torque follow the gust shape closely.An analysis of the time history of rotor thrust and torque during the gust show an increased rotor loading of about 100 % as during constant inflow.Pressure distributions and friction force coefficients reveal that the flow on the rotor blades during maximum gust speed is separated and thus highly in-stationary.Moreover, the effect of accelerating wind speeds was found in the rotor wake as the distance between the vortices is stretched and compressed according to the changes in wind velocity.
As no experiment is available for validation, a final statement about the accuracy of the presented procedure has not been possible.The first mandatory step for future investigations is to find means of validation.Moreover, to achieve a gust transport with other than infinite speed of sound has to be the aim of further research.This may either be realized by adjustments of the resolved-gust-approach presented herein or by implementing the field approach of, for example, Parameswaran et al. (Parameswaran and Baeder, 1997) or the velocity disturbance approach of Reimer et al. (Reimer et al., 2015).A third possibility would be to introduce the fluctuating gust velocities obtained from LES runs which themselves fulfil the continuity conditions.Only then, the procedures of gust computations for wind turbines in THETA are prepared to be extended to respect atmospheric boundary layer flows or for aero-elastic analysis.The method would then enable to gain enhanced knowledge on the flow development and load distribution on rotor loads during gust excitations Bierbooms, W. and Drag, J.: Verification of the Mean Shape of Extreme Gusts, Wind Energy, 2, 137-150, 1999. Castellani, F., Astolfi, D., Mana, M., Piccioni, E., Bechetti, M., and Terzi, L.: Investigation of terrain and wake effects on the performance of wind farms in complex terrain using numerical is a three bladed wind turbine with a rotor radius of 63.0 m and a hub height of 90 m.The rotor has a cut-in wind speed of v ci = 3 m/s and a rated wind speed of v rated = 11.4 m/s.Cut-in and rated rotational speeds are ω ci = 41.4 • /s and ω rated = 72.6 • /s, respectively.Blades are pre-coned and the rotor plane is tilted about β = 5.0 • and not yawed.Along the non-linearly twisted blade 7 different open-access profiles are used.For detailed information please refer to Jonkman et al.(Jonkman et al., 2009).Due to meshing issues the nacelle of the NREL 5MW turbine is neglected while the tower is respected.This approach leads to an error in the flow prediction behind the rotor hub but is supposed to have no impact on the blade loads.The gust simulation is based on the rated wind speed v rated = 11.4 m/s.Air density of ρ = 1.225 kg/m 3 and the kinematic viscosity of ν = 1.82 • 10 −5 m 2 /s is used.

Figure 2 .
Figure 2. Inflow velocity in dependence of physical time t operating gustThe time dependent velocity change of the extreme gust is modelled with(IEC, 2005) 6) and is 5.74 m/s in the given case.In equation 4 σ 1 = 0.11 • v hub is the standard turbulence deviation, λ 1 = 42 m the turbulence scale parameter and D the rotor diameter.u hub represents the velocity at hub height.The velocity u e1 is the average over 10 minutes with a recurrence period of 1 year.It is defined as u e1 = 1.12 • u ref z z hub 0.11 (7) with the reference velocity v ref = 50 m/s as defined in the IEC standard (IEC, 2005) for a wind turbine for the A1 wind class.By respecting a height-independent flow profile, evaluating equation 6 for the given flow conditions and entering equation 7 in equation 5 one obtains the final gust definition

Figure 3 .
Figure 3. Rotor thrust Fx during the gust

Figure 4 .
Figure 4. Rotor torque Mx during the gust

Table 2 .Figure 5 .
Figure 5. Rotor thrust Fx and torque Mx during the gust

Figure 10 .
Figure 10.Pressure distribution of blade with tower blockage effect -inboard section

Figure 14 .
Figure 14.Tip vortex transportation in vertical plane through rotor centre

Table 3 .
Gust induced load during the EOG on rotor in relation to the constant average blade load δu [%] δFx [%] δMx [%] ug(min)[ m/s] Figure 6.Rotor position at minimum (left) and maximum (right) gust velocity; Black blade is blade number 1