NREL is a national laboratory of the U.S. Department of Energy, Office of Energy 

Efficiency & Renewable Energy, operated by the Alliance for Sustainable Energy, LLC.

Contract No. DE-AC36-08GO28308

Atmospheric and Wake 

Turbulence Impacts on Wind 

Turbine Fatigue Loading

Preprint

S. Lee, M. Churchfield, P. Moriarty, J. Jonkman, 

and J. Michalakes 

To be presented at the 50th AIAA Aerospace Sciences Meeting

Nashville, Tennessee

January 9-12, 2012

Conference Paper

NREL/CP-5000-53567

December 2011
NOTICE

The submitted manuscript has been offered by an employee of the Alliance for Sustainable Energy, LLC 

(Alliance), a contractor of the US Government under Contract No. DE-AC36-08GO28308. Accordingly, the US 

Government and Alliance retain a nonexclusive royalty-free license to publish or reproduce the published form of 

this contribution, or allow others to do so, for US Government purposes.

This report was prepared as an account of work sponsored by an agency of the United States government. 

Neither the United States government nor any agency thereof, nor any of their employees, makes any warranty, 

express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of 

any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately 

owned rights. Reference herein to any specific commercial product, process, or service by trade name, 

trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, 

or favoring by the United States government or any agency thereof. The views and opinions of authors 

expressed herein do not necessarily state or reflect those of the United States government or any agency thereof.

Available electronically at http://www.osti.gov/bridge

Available for a processing fee to U.S. Department of Energy

and its contractors, in paper, from:

U.S. Department of Energy

Office of Scientific and Technical Information

P.O. Box 62

Oak Ridge, TN 37831-0062

phone: 865.576.8401

fax: 865.576.5728

email: mailto:reports@adonis.osti.gov

Available for sale to the public, in paper, from:

U.S. Department of Commerce

National Technical Information Service

5285 Port Royal Road

Springfield, VA 22161

phone: 800.553.6847

fax: 703.605.6900

email: orders@ntis.fedworld.gov

online ordering: http://www.ntis.gov/help/ordermethods.aspx

Cover Photos: (left to right) PIX 16416, PIX 17423, PIX 16560, PIX 17613, PIX 17436, PIX 17721

Printed on paper containing at least 50% wastepaper, including 10% post consumer waste.
1

Atmospheric and Wake Turbulence Impacts on Wind 

Turbine Fatigue Loadings

S. Lee1, M. Churchfield2, P. Moriarty3, J. Jonkman4 and J. Michalakes5

National Renewable Energy Laboratory, Golden, CO, 80401

Large-eddy simulations of atmospheric boundary layers under various stability and surface 

roughness conditions are performed to investigate the turbulence impact on wind turbines. In 

particular, the aeroelastic responses of the turbines are studied to characterize the fatigue loading 

of the turbulence present in the boundary layer and in the wake of the turbines. Two utility-scale 

5-MW turbines that are separated by seven rotor diameters are placed in a 3 km by 3 km by 1 km 

domain. They are subjected to atmospheric turbulent boundary layer flow and data is collected 

on the structural response of the turbine components. The surface roughness was found to 

increase the fatigue loads while the atmospheric instability had a small influence. Furthermore, 

the downstream turbines yielded higher fatigue loads indicating that the turbulent wakes 

generated from the upstream turbines have significant impact.

Nomenclature

ABL = atmospheric boundary layer

D = rotor diameter

DEL = damage equivalent load

∆ = filter length scale

δij = Kronecker delta

εijk = alternating unit tensor

Fi = body force for momentum equation

f = frequency

fc = Coriolis parameter

Mip = in-plane blade root moment

Mlss-x = low-speed-shaft torque

Moop = out-of-plane blade root moment

Mtwr-y = pitch moment at the tower base

Mtwr-x = roll moment at the tower base

Myaw = yaw moment at the nacelle yaw bearing

NL = neutrally stable ABL with low surface roughness

NH = neutrally stable ABL with high surface roughness

υ = viscosity

φ = latitude

p = filtered pressure

Pr = Prandtl number

qj = j

th component of temperature flux

ρo = constant air density

SGS = sub-grid scale

τij = fluid viscous and SGS stress

θ = potential temperature

 1 Postdoctoral Researcher, National Renewable Energy Lab, 1617 Cole Blvd., MS3811, Golden, CO, AIAA member

2 Research Engineer, National Renewable Energy Lab, 1617 Cole Blvd., MS3811, Golden, CO, AIAA member

3 Senior Engineer, National Renewable Energy Lab, 1617 Cole Blvd., MS3811, Golden, CO, AIAA member 4 Senior Engineer, National Renewable Energy Lab, 1617 Cole Blvd., MS3811, Golden, CO, AIAA member

5 Senior Scientist, National Renewable Energy Lab, 1617 Cole Blvd., MS1622, Golden, CO, AIAA member
2

θ0 = reference potential temperature

u



i = i

th component of filtered velocity

UL = unstable ABL with low surface roughness

UH = unstable ABL with high surface roughness

ω = planetary rotation rate

zi = boundary layer height

I.Introduction

ESPITE the rapidly growing wind farm installations worldwide, atmospheric and wake turbulence interactions

are not well understood with respect to the fatigue loads on wind turbines. Various numeric modeling 

approaches have been performed. One by Moriarty et al. [1] employed the stochastic turbulence simulator,

TurbSim [2] which was coupled with the aeroelastic code, FAST [3], generated multiple samples of loading data 

under various wind conditions. Recent work by Sim et al. [4] reported load statistics from FAST using a large-eddy 

simulation (LES), with coarse resolution, and found good agreement with the TurbSim inflow technique. Lavely et 

al. [5] reported that coherent structures are an important element in the fatigue loading and strong spatial 

correlations were found with extreme loading events and the passage of turbulent structures. Both the works by Sim 

et al. and Lavely et al. were based on one-way coupling between the LES and an aeroelastic tool and the turbines 

were not modeled within the LES. In this study, atmospheric and wake turbulence impacts on wind turbines are 

investigated using a two-way coupled aeroelastic tool with LES. The turbine blades are modeled as a body force

within the flow solver so that the wakes can be generated and to collect loading data. Neutral and unstable 

atmospheric boundary layers (ABL) are studied under different surface roughness conditions, land-based and 

offshore, to assess the fatigue impacts on the turbine. Furthermore, the loads from the downstream turbines are

investigated to characterize the wake effects of the upstream turbines.

II. Methodology

A. LES Framework

Large-eddy simulation is performed to simulate the atmospheric boundary layer under various surface roughness 

and stability conditions. An incompressible formulation of the Navier-Stokes equation, with the Coriolis force and 

the buoyancy term using the Boussinesq approximation, with the continuity condition and the potential temperature 

transport equation is shown in the following:



i

i

u

x

∂

∂

= 0

(1)

     0

3 3

0 0 0

1 1 ( ) ij i

j i c ij j i i

j i j

u p u u f u g F

tx xx

τ θ θ ε δ

ρ θ ρ

∂ − ∂ ∂ ∂  

+ =− − − + +   ∂∂ ∂ ∂  

(2)

  ( ) j

j

j j

q

u

tx x

θ θ ∂ ∂ ∂

+ =−

∂∂ ∂

(3)

The tilde on the velocity vector denotes the resolved component. The filtered static pressure with the gravity 

potential term is p and ρo is the constant density. Fluid viscous stress and the subgrid-scale (SGS) stress are 

included in the τij term. The εijk is the alternating unit tensor, fc is the Coriolis parameter defined as fc = ω[0, cos(φ), 

sin(φ)], where ω is the planetary rotation rate (7.27x10-5 rad/s) and φ (45o

) is the latitude. The SGS flux is computed 

using the Smagorinsky [6] model with the Smagorinsky constant of 0.13. The filter length scale is defined as

1/3 ∆= ∆ ∆ ∆ ( ) xyz , where ∆x, ∆y and ∆z are the local cell lengths. The buoyancy effect is calculated using the 

Boussinesq approximation, where g is the gravity, δij is the Kronecker delta, θ is the resolved potential temperature, 

θ0 is the reference temperature taken to be 300K. The last term, Fi, is the force field generated by the blade model,

D
3

which is discussed in the following section. The transport equation for the resolved potential temperature equation is 

shown in (3), where qj represents the temperature flux (4) as follows:



Pr

SGS

j

SGS j

q

x

υ ∂θ = − ∂

(4)

where υSGS and PrSGS are the SGS viscosity and turbulent Prandtl number following the Moeng’s [7] formulation.

Further details can be referred to in [8].

B. Actuator Line Method and FAST

The wind turbine used in the present study is the National Renewable Energy Laboratory's (NREL) 5-MW 

reference turbine [9], which is based on the REpower 5M and Multibrid M5000 turbines. This upwind, horizontalaxis turbine has a three-bladed rotor, 126 m in diameter with a hub height of 90 m. For simplicity, the pitch and the 

yaw angles were fixed at 0o relative to their neutral frame of reference, while a variable speed torque controller was

activated [9]. The torque controller follows the optimal torque output, which is proportional to the square of the 

rotor's rotational speed, with a gain coefficient of 0.0255764 N-m/RPM2

. The rated power of 5 MW occurs at a wind 

speed of 11.4 m/s and a rotor speed of 12.1 RPM. At a hub-height-average flow speed of 8 m/s, the upstream 

turbine rotor rotates at approximately 9 RPM producing power in the range of 2 MW.

The turbine model consists of an actuator line representation [10] of the turbine blades coupled with FAST [3]. 

FAST employs a combined modal and multi-body dynamics formulation. The flexible blades and the tower are 

characterized using a modal representation that assumes small deflections. Within FAST, each turbine blade is 

represented by a set of discrete elements along the blade that move in space as the turbine blades rotate and flex. 

FAST uses the local flow field at each blade element to compute aerodynamic forces via the dynamic airfoil 

coefficient lookup tables in the AeroDyn module [11]. To capture dynamic stall, the Leishman and Beddoes [12] 

model is employed. The resulting structural response and the updated blade-element positions are computed and

include the deflection and the rotations of each blade. Using the updated blade-element positions, the aerodynamic 

forces are projected back onto the flow field as a body force field using Sørensen and Shen's [10] actuator line 

method with the Gaussian width equal to the twice the grid cell length. The size of the width is based on a previous 

study [8]. When coupled with the LES, the wake calculation is disabled in AeroDyn [11] since the wake and axial 

induction are captured by the LES.

C. Numerical Approach

The above incompressible Navier-Stokes with the Coriolis force and the potential temperature flux equations are

solved using OpenFOAM [13]. The equations are discretized using an unstructured collocated finite-volume 

formulation. The second-order central-differencing scheme is employed and uses Rhie and Chow [14] interpolation 

to avoid checkerboard pressure-velocity decoupling. For the time advancement, PISO (Pressure Implicit Splitting 

Operation) is used, with three sub-step corrections to maintain a second-order temporal accuracy. FAST’s 

aeroelastic calculation employs the fourth-order Adams-Bashforth predictor and Adams-Moulton corrector time 

integration scheme. Further details are provided in [3].

D. Computational Domain

The size of the computational domain is 3000 m by 3000 m in area and 1000 m in height as shown in Fig 1. The 

two 5-MW NREL turbines are placed in the middle of the domain, where the bisecting point matches with the center

of the bottom wall. The separation distance between the two turbines is seven rotor-diameters (7D), and their 

positions are rotated 30o counterclockwise around the bisecting point relative to the x axis, to be inline with the 

general wind direction. The grid is refined in two steps in the rectangular region and the turbines are placed as 

shown Fig. 1. Outside of the refinement zone, the grid cell size is uniformly 10 m. In the first stage of the mesh 

refinement, the intermediate zone with a length of 19D and a 5D width, the grid cell size is divided in half yielding a 

5 m resolution. The size of the inner zone is indented by 1D on each of the five sides (without the bottom wall side)

and the resolution is increased to 2.5 m grid cells, which yield a total of 26 million cells. As the coherent structure 

presences are intermittent in an ABL, the turbines are placed in various locations to enhance the statistical
4

Fig. 1 Schematic of the computational domain

independence of the FAST data. The turbine locations are shifted laterally (perpendicular to the wind direction) 5D

up and down from the center point and the refinement zone locations are shifted accordingly as indicated by the 

dashed-lined boxes in Fig. 1. Sets of 1800 s data are collected for each turbine at the three locations (yielding 5400 s

of data), which is well over the minimum limit of one continuous 60-minute or six 10-minute samples recommended

by the IEC-61400-1 standard [15]. Note that for each simulations, a total of 2100 s of data is recorded of which the 

initial 300 s transient data are disregarded. The surface stress and the temperature flux conditions are imposed at the 

surface boundary following the Moeng’s approach [7]. The friction velocity is approximated using the MoninObhukov similarity theory. Zero stress and temperature flux is imposed for the top boundary while a periodic 

condition is set for the north and the south faces and inflow and outflow conditions for the west and east side, 

respectively. One of the disadvantages of periodically reusing the outflow data for the inflow boundary condition is 

that the turbulent structures reappear at the same location [16]. To circumvent this problem, the wind direction is set 

at a skewed angle (30o

) to break up the spatial correlation of the large scale structures by allowing them to “drift”. 

With this condition, the large scale structures exiting the north side reappear on the south side. To provide a 

turbulent flow field at the inflow boundary, a precursor ABL simulation data is used [8].

E. Atmospheric Boundary Layer

Four cases of ABLs are performed and can be divided into two groups: neutrally stable (N) and unstable (U) 

boundary layer condition. The average surface temperature flux, s < > θ ω′ ′ , is set to be 0.04 K-m/s for the unstable 

case while a fixed value of zero is used for the neutral case. Within the different stability cases, two sub-categories 

of low (L) and high (H) surface roughness are investigated. The roughness height, z0, is 0.001 m and 0.2 m for L and 

H are representative of offshore and land-based conditions, respectively [17]. Therefore, the abbreviations for 

“neutral stability and low roughness,” “neutral-high,” “unstable-low,” and “unstable-high” are NL, NH, UL, and 

UH. However, in all cases, the potential temperature profile is set to 300 K from the surface to 700 m, which is 

followed by a prescribed capping inversion of up to 800 m allowing the potential temperature to linearly increase to 

308 K. Above this height, the temperature rises at a constant rate of 0.003 K/m. The boundary layer height, zi, is 

defined as the altitude at which the time- and horizontal-averaged temperature flux is minimum. The horizontal 

average wind speed at the hub height (90 m) is adjusted to remain at 8 m/s at a 30o (northeasterly) direction by

adjusting the driving mean pressure gradient.

Instantaneous flow field of the four ABL cases described above are shown in Fig. 2. The blue streaky structures

are the iso-surface of the x-direction component with a velocity of -1.25 m/s. Compared to NL (Fig. 2a), 

significantly higher numbers of elongated low speed streaks are present for NH due to increased roughness. Updrafts 

due to the temperature flux at the bottom surface are demarcated by the abundantly distributed red structures for UL
5

a) b)

c) d)

Fig. 2 Iso-surface of instantaneous flow of the horizontal component at -1.25m/s (blue) and the vertical 

component 1m/s (red) for a) NL, b) NH, c) UL and d) UH

a) b)

Fig. 3 a) mean streamwise velocity profile and b) energy spectra for various ABL conditions
6

Fig. 4 Two NREL 5-MW turbines subjected to NL atmospheric conditions showing the instantaneous 

streamwise velocity with streamlines

and UH cases, which are the iso-surface of the vertical velocity at 1 m/s. Note that the locations of the updrafts are 

correlated with that of the low speed streaks in the unstable cases.

The x-direction components of the mean velocity profiles for the four ABL conditions are shown in Fig. 3a. In 

general, the vertical mixing, from the updrafts in the unstable ABL cases (UL and UH), results in flatter profiles 

compared to those from the neutral cases (NL and NH). The inflections of the curves above z/zi ~ 0.7 are caused by 

the gravity waves in the upper boundary layer. The two-dimensional horizontal velocity spectra taken at the turbine 

hub height, shown in Fig. 3b, capture the inertial range similar to the works by Moeng and Wyngaard [18], with an 

explicit Gaussian filter. Moeng and Wyngaard [18] reported that a longer inertial range was captured when a sharp

wave number cut-off filter was used. However, due to the finite volume nature of the present solver, further 

improvements in capturing the inertial range was not feasible without significantly increasing the grid resolution.

III.Results

An overview of the flow field in a low-roughness neutral stability condition is shown in Fig. 4. A wake downwind 

of each turbine is formed, and the shear layer at the edge of the wake quickly becomes unstable. This instability acts 

to break down the organized spiral of the tip and root vortices within a rotor diameter downwind of the rotor plane. 

High momentum fluid is entrained into the low momentum wake by the rotational flow, which promotes wake 

recovery. Turbulent structures at the edge of the wake appear to dissipate as they convect downstream and the shear 

layer weakens. The streamlines (green) are shown to illustrate both the stronger velocity fluctuations in the wakes 

and the rotation of the wakes counter to the clockwise-turning (as viewed from upwind) turbine rotors. In this 

section, the impact on structural loading from atmosphere and wake turbulence coupled with the various ABL

conditions is investigated. In particular, the spectral analysis is performed to characterize the dominant frequencies 

and the fatigue load is estimated from the Rainflow cycle count [19].

A. Time History of Blade Root Moment and Tower Bending Moment

The schematic of the measured moments in various parts of the wind turbine are shown in Fig. 5. Time histories

of the blade-root out-of-plane moment (Moop) and the tower-base fore-aft moment (Mtwr-y) ranging from 300 s to 

2100 s of data are shown in Fig. 6. For the four ABL conditions, strong variations in the load signals are observed. 

Note that the cyclic signals with a period of approximately 6.5 s are due to the shear in the flow and the gravity,

which deflects the blade further in the upper hemisphere of the rotational cone. The impact of the atmospheric 

turbulence is shown in the upstream Moop case (Fig. 6a). Both the surface roughness and the atmospheric stability
7

Fig. 5 Moments measured at various 

locations

a) b)

c) d)

Fig. 6 Time history of moment signals for a) Moop for WT1, b) Moop for WT2, c) Mtwr-y for WT1 and d) Mtwr-y

for WT2
8

generate large fluctuating signals, which are due to the convections of large scale coherent structures. The unstable 

atmospheric cases (UL and UH) tend to yield large undulations at lower frequencies, which is attributed to the 

updrafts of the warm plumes. The downstream turbines yield different load signals due to the altered flow field 

perturbed with vortices from the blades and the meandering wakes. However, in general, the load signals are lower 

in magnitude due to the velocity reductions. Similar signal patterns are observed with Mtwr-y, without the high 

frequency component due to signal damping and the stationary structure of the tower. However, the out-of-plane 

moments at the blade root exert forces at the tower base at greater magnitude. The moments imparted on other parts 

of the turbine exhibit similar fluctuating signals. Thus, a spectral analysis is performed to characterize the structural 

response to the flow in the next section.

B. Spectral Analysis

The peaks beginning from the lower frequencies of in-plane (Mip) and out-of-plane (Moop) moment at the blade 

root, shown in Fig. 7, correspond to the harmonics of the rotor rotational frequency (1P, 2P, 3P, and so on). The 

definitions of the moments can be referred to Fig. 5. The amplitudes of these harmonics decrease with increasing 

frequency. The turbulent flow field exerts unsteady forces onto the blade with a large range of frequencies so that 

the harmonics of the rotational speed are captured in these spectra. Note that the peak frequencies are lower for the 

downwind turbines due to the reduced velocity field decreasing the rotor speed. The fact that the peak corresponding 

to 1P is largest indicates that the greatest load cycle amplitudes are due to asymmetric blade loading. This 

asymmetry is due to the vertical wind shear or the passage of turbulent structures through one side of the rotor disk. 

The peak corresponding to 1.1 Hz is the first edgewise pitch mode natural frequency (not present in Moop) and the 

2.02 Hz frequency is the second blade flapwise collective mode as shown in Fig. 7a and 7b. The harmonics of 3P 

correspond to the peaks in the nacelle bearing yaw moment, Myaw, as shown n Fig. 8a. Yaw moments of frequency 

3P occur when an unbalanced aerodynamic force is imposed laterally upon the three blades. For instance, as a 

turbulent structure passes through one side of the rotor disk, the thrust loads on the three blades are not causing a 

cyclically-varying yaw moment equally. The two right-most peaks, with corresponding frequencies of 1.85 Hz and 

2.16 Hz, are due to the second asymmetric flapwise pitch and yaw modes, respectively. Perhaps the overhang 

distance from the yaw bearing to the hub makes the yaw moment induced by these asymmetric flapwise modes large 

enough that they are discernable. Similar mechanism of the 3P excitation peaks are observed for the low-speed-shaft 

torque; Mlss-x (Fig. 8b). The peak corresponding to the first drivetrain torsion mode (f=2.26Hz) is larger than the 3P 

counterpart, which indicates that the variable speed controller should be improved to dampen this pronounced peak. 

The second edgewise mode is not modeled in FAST. The side-to-side moment and the for-aft moment at the tower 

base, Mtwr-x and Mtwr-y, are largely affected by the 3P harmonics. These excitations can impose large moments due to 

the tall tower. In particular, three natural frequencies associated with the first edgewise yaw mode (f=1.25Hz), which 

appears to have the strongest influence, the second asymmetric flapwise pitch mode (f=1.85), and the first drivetrain 

torsion mode (f=2.26Hz) appear in the energy spectra as shown in Fig. 8c. Only the second asymmetric flapwise 

pitch mode is captured in the Mtwr-y spectra since other natural frequencies are related to the edgewise mode.

a) b)

Fig. 7 Energy spectra for blade root moments in a) in-plane: Mip and b) out-of-plane: Moop directions
9

a) b)

c) d)

Fig. 8 Energy spectra of a) yaw moment at the tower-top yaw bearing: Myaw, b) torque at the low-speed-shaft 

of the rotor, c) side-to-side moment at the tower base: Mtwr-x and d) fore-aft moment at the tower base: Mtwr-y

The difference regarding the dominant peaks relating to 1P and 3P frequencies in Figures 7 and 8 is because the 

individual blades can only experience the dominant 1P event for an asymmetric aerodynamic loading while the 

remaining structure receives the impacts from all three blades. While the spectral peaks related to the harmonics of 

1P and 3P and the natural frequencies are found, the distinction between various ABL conditions is difficult to 

analyze. Therefore, the probability of extreme loading events and the damage-equivalent loads are discussed in the 

following sections.

C. Short-term Load Distributions

The probability of exceedance for Moop, Myaw, Mtwr-x and Mtwr-y (Fig. 9) are computed from the raw rainflow cycle 

counted data. In general, surface roughness has a significant impact on the occurrence of extreme loads, especially 

for the neutrally stable atmosphere with high roughness (NH). For the upstream turbine, the probabilities of the 

strong loads were roughly equal regardless of the stability of the atmosphere in the low roughness cases (NL-WT1 

and UL-WT1). The downstream turbines yielded stronger loads than the upstream turbines (with the exception of 

NH-WT2 for Moop), which indicates that the added turbulence by the upstream wake can have a significant impact 

on the downstream turbines, even at the separation distance of 7D.

D. Fatigue Load Estimation

The damage equivalent loads (DEL) computed from stress range histograms derived from the rainflow cycling 

counting algorithm are shown in Fig. 10, where the Wöhler exponent for the blade is 10, which corresponds to a 
10

a) b)

c) d)

Fig. 9 Probability of exceedance for rainflow cycle counts of a) Moop, b) Myaw, c) Mtwr-x and d) Mtwr-y

typical composite material used for wind turbine blades, while the remaining components have exponents of three, 

assuming steel material. The Goodman method [20] was employed for the cycle midpoint correction with ultimate 

loads based on maximum values from the FAST data for each component in question. The maximum loads are 

multiplied by integers until convergence is achieved. In the case of Mip, the damage equivalent load (DEL) displayed 

less sensitivity to the ABL conditions. The error bars are twice the standard deviation yielding 95% certainty. 

However, the DELs were predicted to be smaller for the downstream turbines, which may be attributed to the lower 

wind speeds that reduce the rotor rotational speed. Relatively, these may decrease the turbulence intensity for the inplane component. Conversely, strong variations are observed for the other cases as shown in Fig. 10b ~ 10f. In 

general, the surface roughness increases the DELs, especially for the downstream turbines. However, the NH case 

includes exceptions for Moop and Mlss-x, which has lower DEL for the downstream turbines suggesting that the high 

atmospheric turbulence mitigates the loading impacts by the upstream wake turbulence.

The DELs for the downstream turbines are consistently higher for Myaw, Mtwr-x and Mtwr-y, which suggests that the 

wake effect on the fatigue loading on the tower is quite significant. The atmospheric stability, on the other hand, 

appears to have weak correlation with the increase or decrease of DELs. The increase in the downstream turbine 

DELs for the high roughness cases is smaller than that of the low roughness suggesting that higher turbulence 

intensity in the atmosphere due to roughness may reduce the relative loading impacts in the wake.

IV.Conclusion

Large-eddy simulation of incompressible Navier-Stokes with Coriolis force coupled with temperature flux 

equations is performed under four different atmospheric conditions. In particular, atmospheric cases of neutrally 
11

a) b)

c) d)

e) f)

Fig. 10 Damage equivalent load for a) Mip, b) Moop, c) Mlss-x, d) Myaw, e) Mtwr-x and f) Mtwr-y

stable, low roughness (NL); neutrally stable, high roughness (NL); unstable, low roughness (UL); and unstable, high 

roughness (UH) are studied to investigate the impact of surface roughness and the stability of the atmosphere on 

structural loadings of wind turbines. In addition, studies on wake effects on the downstream turbines are conducted 

in this study. In all four cases, the average hub-height wind speeds are set at 8 m/s which are below the rated wind 

speeds.

Two NREL 5-MW turbines, separated by seven rotor-diameters (7D), are placed in a 3 km by 3 km by 1 km 

domain and subjected to the atmospheric boundary layer. Structural load data are collected using the aeroelastic tool, 
12

FAST, which is coupled with the LES solver. Large fluctuations are observed in the moment signal at the blade root 

and the tower-base locations, which are tied to the convection of turbulent structures. Spectral analysis indicates that 

for individual blades, the dominant energy peaks corresponded to the harmonics of 1P frequencies while the other 

components of the turbine structure (i.e., the yaw moment at the nacelle bearing, low-speed shaft and tower base) 

yield harmonics of 3P frequencies. The peak frequencies for the downstream turbines are smaller compared to that 

of the upstream ones due to the reduced velocity caused by the upstream turbine wake. However, the peak 

frequencies occur at the same natural frequencies of the structure.

Histograms of stress range calculated using the rainflow cycle count algorithm indicate that the surface 

roughness has a strong influence on the extreme loadings on the turbine, where higher roughness leads to increased 

damage equivalent loads (DELs). Pronounced increases are observed in the NH case with respect to moments at the 

low-speed shaft, yaw bearing, and the tower base. The atmospheric instability appears to have a weak impact on the 

DELs. In general, the downstream turbines yield higher DELs, indicating that the turbulent wakes from the upstream 

turbines can have significant impact even at 7D separation distance. However, the in-place blade root moments (Mip) 

result in lower DELs for the downstream turbines due to the wind speed reduction. In addition, both the out-of-plane 

blade-root moment (Moop) and the low-speed-shaft torque (Mlss-x) yield lower DELs for the downstream turbines in 

the NH case, which suggests that the atmospheric turbulence mitigates the loading impacts by the upstream turbine 

wakes. Note that the present study is conducted at below rated wind speed such that different conclusions can be 

expected for wind speeds above-rated conditions. Further investigations will be reported in future studies.

Acknowledgments

The authors gratefully acknowledge the financial support provided by the NREL’s LDRD program. All of the 

computations were performed on Red Mesa at the National Renewable Energy Laboratory. Discussions and 

guidance by Marshall Buhl at NREL is greatly appreciated.

References

[1] Moriarty, P. J., Holley, W. E., and Butterfield, S. P., “Extrapolation of Extreme and Fatigue Loads Using 

Probabilistic Methods”, NREL/TP-500-34421, (2004)

[2] Kelley, N. D., and Jonkman, B. J., “Overview of the TurbSim Stochastic Inflow Turbulence Simulator”, 

NREL/TP-500-41137, (2007)

[3] Jonkman, J. M., and Buhl, M. L. Jr., “FAST User’s Guide”, NREL/EL-500-38230, (2005)

[4] Sim, C., Manuel, L., and Basu, S., “A Comparison of Wind Turbine Load Statistics for Inflow Turbulence 

Fields based on Conventional Spectral Methods and Large Eddy Simulation”, AIAA 2010-829, (2010)

[5] Lavely, A. W., Vijayakumar, G., Kinzel, M. P., Brasseur, J. G., and Paterson, E. G., “Space-Time Loadings on 

Wind Turbine Blades Driven by Atmospheric Boundary Layer Turbulence”, AIAA 2011-635, (2011)

[6] Smagorinsky, J., “General Circulation Experiments with the Primitive Equations”, Monthly Weather Review, 

91 (1963), pp. 99-164

[7] Moeng, C. H., “A Large Eddy Simulation Model for the Study of Planetary Boundary Layer Turbulence”, J. 

Atmospheric Sciences, 41 (1984), pp. 2052-2062

[8] Churchfield, M. J., Lee, S., Michalakes, J., and Moriarty, P. J., “A Numerical Study of the Effects of 

Atmospheric and Wake Turbulence on Wind Turbine Dynamics”, J. Turbulence (submitted)

[9] Jonkman, J., Butterfield, S., Musial, W., and Scott, G., “Definition of a 5-MW Reference Wind Turbine for 

Offshore System Development”, NREL/TP-500-38060 (2009)

[10]Sørensen, J. N., and Shen, W. Z., “Numerical Modeling of Wind Turbine Wakes”, J. Fluids Engineering, 124, 

(2002), pp. 393-399

[11]Moriarty, P. J., and Hansen, A. C., “AeroDyn Theory Manual”, NREL/EL-500-36881, (2005)

[12]Leishman, J. G., and Beddoes, T. S., “A Semi-Empirical Model for Dynamic Stall”, J. American Helicopter 

Society, 34 (1989), pp. 3-17

[13]OpenFOAM – The Open Source CFD Toolbox, User’s Manual, Version 1.7.1; OpenCFD Ltd., 9 Albert Road, 

Caversham, Reading, Berkshire RG4 7AN, UK

[14]Rhie, C. M., and Chow, W. L., “Numerical Study of the Turbulent Flow Past an Airfoil with Trailing Edge 

Separation”, AIAA J. 21 (1983), pp. 1525-1532

[15]International Electrotechnical Commission, Wind Turbines – Part 1: Design Requirements, IEC-61400-1, 

Edition 3.0, (2007)
13

[16]Ghosh, S., Choi, J.-I., Edwards, J. R., “Numerical Simulations of Effects of Micro Vortex Generators Using 

Immersed Boundary Methods”, AIAA J. 48 (1) (2010), pp. 92-103

[17]Stull, R. B., “Meteorology for Scientists and Engineers”, 2nd edition, Brooks Cole (1999)

[18]Moeng, C. H., and Wyngaard, J. C., “Spectral Analysis of Large Eddy Simulations of Convective Boundary 

Layer”, J. Atmospheric Sciences 45 (1988), pp. 3573-3587

[19]Downing, S. D., and Socie, D. F., “Simple Rainflow Counting Algorithms”, Int. J. Fatigue 4 (1) (1982), pp. 31-

40

[20]Sutherland, H. J., “On the Fatigue Analysis of Wind Turbines”, SAND99-0089, Sandia National Laboratories, 

(1999)