We consider the Large Eddy Simulation (LES) of incompressible flows. It is performed using an in-house developed fourth-order finite differences code (Duponcheel et al., 2008; Duponcheel et al., 2014), formulated in primitive variables (i.e. in velocity and pressure) and using a staggered mesh arrangement. Using a Cartesian coordinate system, the Navier-Stokes equations (supplemented by a subgrid scale (SGS) model) are formulated as ∂ u i ∂ x i = 0 , (1) ∂ u i ∂ t + u j ∂ u i ∂ x j = − 1 ρ ∂ p ∂ x i + ν ∂ 2 u i ∂ x j ∂ x j + ∂ τ i j M ∂ x j + f i ρ , (2) where i , j ∈ 1,2,3 correspond to the streamwise (x, u), vertical (y, v) and lateral (z, w) directions and velocity components, respectively; p is the pressure and f _i is the external body force (per unit volume) acting on the flow to model the wind turbine (discussed in Section 2.3); ρ and ν are the fluid density and molecular viscosity, respectively; τ i j M is the SGS stress tensor model. We here use a classical Smagorinsky model (Smagorinsky, 1963).

The equations were initially discretized in space using the fourth-order finite differences scheme of Vasilyev (2000), which is such that the discretized convective term conserves the kinetic energy on Cartesian meshes. Recent developments saw the addition of wall modeling strategies. This required the convective term to be rewritten; this new version thus preserves the discrete kinetic energy up to fourth order (in the absence of wall modeling procedure, see Thiry (2017) for details) and it conserves the momentum exactly on uniform grids. In this work, only the new version of the code is used as we consider simulations of rotors in an atmospheric boundary layer.

Equations 1 and 2 are solved using a fractional-step method with the “delta” form for the pressure (Lee et al., 2001). The Poisson equation for the pressure is solved using a multigrid solver with a Gauss-Seidel smoother. The time integration is carried using a second-order Adams-Bashforth scheme.

There are two main approaches for generating turbulent inflows in LES. The first one implies using synthetic turbulence field generation techniques, while the second one consists in generating turbulence using a precursor LES.

In the first approach, the synthetic turbulence field is pre-generated from a chosen model and is then imposed over a mean profile at the inlet of the simulation of interest (here denoted the wind farm (WF) simulation). The main advantage of this method is the low computational cost overhead compared to the WF simulation. The most advanced models match physical spectra and moments that are assumed to evolve rapidly enough into realistic turbulence, such as that developed by Mann (1998) and widely used in the wind energy domain. However, Munters et al. (2016) demonstrated that the adaptation length required for the larger structures to develop inside the WF simulation was very important. Those large turbulent structures are known to contribute to the wake meandering process (Espana et al., 2011) and are thus particularly important to capture. This large establishment length sensibly increases the computational domain and thus the computational costs. Consequently, the second approach is used in the present work, i.e. generating the turbulent inflows using a precursor LES.

The precursor technique relies on an auxiliary wind flow simulation, performed on an independent domain without wind turbines, for generating turbulent inflow conditions. This auxiliary simulation is here denoted the Atmospheric Boundary Layer (ABL) simulation. When it reaches a statistically-steady state, either yz-planes of instantaneous velocities are saved to generate an independent data base inflow data, or this auxiliary simulation runs concurrently with the WF simulation (concurrent simulations), which allows inflow conditions to be directly transferred to the inlet of the main domain. The use of precursor techniques has increased in wind energy in recent years (Wu and Porté-Agel, 2015; Munters et al., 2016), and is strongly recommended for wind farm simulations, in order to obtain faithful flow structures impacting the wind turbines (Park et al., 2014).

In this work, concurrent simulations are considered, in order to avoid saving a large amount of inflow velocity planes. For the ABL simulation, periodic conditions are enforced in the horizontal directions, while a no-through flow is imposed at the top. At the bottom, a wall stress model (Thiry, 2017) for a rough wall is applied, to compute the surface shear stress as a function of the LES velocity field at the third vertical grid point and a roughness length, y ₀. The flow is driven by a pressure gradient that is continuously adapted in order to conserve an imposed mass flow, itself preliminary determined to achieve the desired hub velocity and turbulence intensity (TI). For the WF simulation, outflow conditions are then imposed at the outlet boundary, while periodic conditions are kept in the spanwise direction. Again, a wall model and a no-through flow condition are imposed at the bottom and top, respectively. The code allows different resolutions and streamwise domain sizes for the ABL and the WF simulations.

In the present work, thermal and stratification effects are not included; the resulting simulations are thus equivalent to truly neutral atmospheric boundary layer (ABL) cases.

The wind turbine model targets the context of accurate and still affordable LES of wind farms, we thus choose a rotating Actuator Disk (AD), discretized in the radial and tangential directions, and for which the aerodynamic forces are estimated from local velocities. The AD method is here improved with a tip-loss factor: this correction is based on the commonly used Glauert tip-loss factor, for which the infinite upstream velocity is estimated based on the 1D momentum theory, applied on each AD elements (Moens et al., 2018). The AD forces are finally translated into bulk forces for the LES equations (f _i in Eq. 2) by using a second-order mollifier, the M ₄ kernel (Monaghan, 1985). This AD approach, supplemented with the tip-loss factor, has been verified with a higher fidelity approach in a finer resolution, a Vortex Particle-Mesh method coupled to immersed lifting lines (Chatelain et al., 2017; Caprace et al., 2020): we refer to Moens et al. (2018) for the comparison results and additional methodological details.

In the present work, we will investigate four different moments: the blade flapwise and edgewise root bending moments and the rotor yaw and tilt moments. They are evaluated based on the aerodynamic contributions only. The flapwise and edgewise moments on blade b are denoted M _f,b and M _e,b , respectively, and are evaluated as M f , b = ∫ R h u b R t i p r − R h u b F b . e n d r , (3) M e , b = ∫ R h u b R t i p r F b . e θ d r , (4) where F _b is the aerodynamic blade force per length; e _n and e _θ are the unit vectors normal and tangential to the rotor plane, respectively; R _hub and R _tip are the radii at the hub and the tip locations. The yawing and the tilting moments, denoted M _y and M _t , respectively, are expressed as M y = ∑ b = 1 N b e y . ∫ R h u b R t i p x × F b d r , (5) M t = ∑ b = 1 N b e z . ∫ R h u b R t i p x × F b d r , (6) where N _b is the number of blades; x is the position from the hub center and the considered force point; e _y and e _z correspond the yaw and tilt axes, respectively. The present study assumes that only the normal force component contributes to the yaw and tilt moments. Indeed, as mentioned above, the moments are calculated based on the aerodynamic contributions only. For the wind conditions considered in this work, the tangential forces will be roughly 10 times lower than the normal forces. Their contribution to M _y and M _t is thus expected to be less important than that related to the normal forces. This is also reinforced by the fact that the tangential forces are multiplied by the shaft length or a fraction of the shaft length in the M _y or M _t calculations (if we consider investigating the fatigue damage on the main shaft). This shaft length is usually small. In comparison, the normal forces are multiplied by the local radius, which is higher than the shaft length for the large part of the blade. If the gravitational loads are added in the computation of M _y and M _t , the contribution of the tangential forces might be no more negligible. However, for a well-balanced rotor, the contribution of the three blades cancels the effect of the gravity on the yawing moment. The gravitational loads play a role on the mean value of the tilting moment, and would add an offset in the temporal signal presented in this paper. As we only account for the cycle amplitudes in the equivalent moments (see Section 4.4), the addition of the gravity should not impact the fatigue analysis of the tilting moment.

For discrete line type approaches, the instantaneous blade loading is directly available and can be integrated along the blade span. For the AD method, we propose to replicate the blade trajectories through the disk and recompute blade-attached aerodynamic forces. At each particular radial position, r, and time, t, the forces F _b acting on the blade b are thus interpolated from the disk forces of the two adjacent disk elements at the same radius r. The instantaneous blade loading is directly available and can be integrated along the blade span to estimate the root bending moments. This procedure can be performed when the simulation is running, in order to avoid the save of a large amount of AD data and additional post-processing computations. Details about the procedure can be found in Moens et al. (Forthcoming 2022b), where it was shown that the fatigue loads estimated using this methodology were in excellent agreement with those obtained using a Vortex Particle-Mesh method coupled to immersed Lifting Lines (Chatelain et al., 2017; Caprace et al., 2020). This methodology has also been used in Moens et al. (Forthcoming 2022a), in order to investigate the individual pitch control within an AD framework, and again, the AD method provides blade moments that are really close to those obtained using a discrete blade type approach. This clearly positions the proposed approach as an efficient method for recovering pertinent and accurate information about fatigue loads.

The simulation tool is not coupled to an aeroelastic code, and the blades thus act as rigid bodies.

We choose the most robust technique of Coudou et al. (2018) to perform tracking of the wake centroid, (y _c , z _c ), in cross-flow planes located downstream of a wind turbine. It is based on finding the minimum of a convolution product y c , z c = arg min p y , z * f G y , z (7) between the available power density in the flow p y , z = 1 2 u | u | 2 (8) and a Gaussian masking function f G y , z = A e x p − y − μ y 2 2 σ y 2 + z − μ z 2 2 σ z 2 (9) with A = 1, μ _y,z is equal to the wind turbine center, and σ _y and σ _z are set to 0.25D and 0.5D, respectively. Physically, this convolution method relies on the computation of the wind power inside a disk (with a Gaussian weighting) located in a cross-flow plane and with a diameter equal to the rotor diameter D. This disk center position can be shifted in the y − and z − directions; the wake center corresponds to the disk position for which the available power is minimum (Vollmer et al., 2016). Finally, by tracking the wake centroid in several downwind cross-flow planes, the wake centerline can be obtained.

For the wind farm simulations, the time-averaged freestream velocity profile is substracted from the velocity field before computing the convolution.

In order to verify the performances of the code for simulating large wind farms, we first compare the time-averaged power output of the present results to available observations in the reference paper of Barthelmie et al. (2009). The results are moreover compared to an other set of LES performed by Wu and Porté-Agel (2015). The Horns Rev measurements, the different LES, as well as the considered size of wind sectors, are detailed in the next paragraphs.

We now investigate the history of the moments for several wind turbines and two representative and “extreme” wind directions: α = 270°, corresponding to the fully-aligned configuration, and α = 262°, where the rotors can be out of the upstream wakes during large periods. We consider four rotors, located in R ₁, R ₂, R ₆ and R ₁₀ in the middle line L ₄ (see Figure 1 for the locations of those wind turbines). These rows are representative of the typical behaviors inside the wind farm, with rotors impacted by an unperturbed wind flow (R ₁) or by a very well-structured wake (R ₂) and with machines located deeper in the wind farm (R ₆ and R ₁₀). For α _w = 262°, the row R ₁₀ may also correspond to situations where the machine is impacted by the wakes of the adjacent line.

We report on Figures 6, 7 the instantaneous signals of the root bending moments and the yaw and tilt moments. The histories of the flapwise and the edgewise root bending moments are shown for one blade, and are denoted M _f,1 and M _e,1, respectively. For the sake of clarity, we only show the instantaneous signals for 15 min. The time is made dimensionless by multiplying it by U ̄ h u b / D and is denoted t ^*. The filtered signals, at the rotation frequency (f _Ω) for M _f,1 and M _e,1, and at three times f _Ω for M _y and M _t , are added on Figures 6, 7 in order to highlight the low frequency fluctuations of the moments.

For rows R ₂, R ₆ and R ₁₀, we also report on Figures 6, 7, the horizontal position of the wake center (denoted Z _wm ) just upstream of the rotor, in order to discuss the moment behaviors with respect to the upstream wake movement. Z _wm is calculated as the relative horizontal position of the upstream wake with respect to the center of the rotor, and normalized by the Vestas diameter. For the tilt moment, it is more interesting to discuss its behavior as a function of the vertical upstream wake movement, denoted Y _wm . Y _wm is calculated as the relative vertical position of the upstream wake with respect to the hub height, and, again, it its normalized by the rotor diameter.

Additionally, we report on Figure 8 snapshots of the instantaneous streamwise velocities at several times, in order to support our discussion of some particular behaviors of the moments for the different machines. We also add the horizontal wake centerlines on those velocity slices.

We structure our discussion in terms of 1) the depth of the turbines, and 2), the center of partial wake impingement (respectively, 270° or 262°).

Rows R ₁ and R ₂: For α _w = 270°, there is an increase in the amplitude of the high frequency oscillations of the root bending moments between R ₁ and R ₂ (see Figures 6A,B). Those moment fluctuations, denoted here as the 1P oscillations, are characterized by a frequency equal to that of the blade rotation and are due to the blade passage through areas of higher and lower velocities over one rotation. For row R ₁, those oscillations occur when the rotor is partially immersed in a gust or in a lull and when the blade crosses the sheared wind. For row R ₂, those moment fluctuations are mainly due to the presence of an upstream wake, which, at times, partially immerses the rotor. An example of the wind turbine in a half-wake situation is shown in Figure 8B for t ^* = 126, which can be correlated with the high amplitudes of the high frequency moment fluctuations. For the considered wind conditions ( U ̄ h u b of 8 m/s and low TI), the wake downstream of the first wind turbine in a line is very “well-structured”, with a velocity deficit that is clearly more marked than the ambient wind fluctuations (as highlighted in Figure 3). The amplitude of the moment fluctuations for a rotor impacted by such a wake is thus higher than that observed when the machine is partially impacted by a gust. Moreover, Figure 6 shows that the meandering of the wake upstream of row R ₂ is weak, with very small oscillations of the wake center: the rotor of R ₂ is thus often in half wake situations, leading to many episodes of large-amplitude high-frequency moment fluctuations. This is also translated by a decrease of the time-averaged values of M _f,1 and M _e,1.

We also notice on Figures 6A,B, for row R ₂, lower frequency fluctuations of M _f,1 and M _e,1 (highlighted in the filtered signals). This is due to the transition between full and half wake situations, characterized by a decrease and an increase in the mean values of M _f,1 and M _e,1.

For the yaw and tilt moments and α _w = 270° (see Figures 6A,B), there is also an increase in the amplitude of the fluctuations related to the blade rotation (the frequency of those oscillations is here equal to three times the rotation frequency (3P), as opposed to above (1P)). There is also a significant increase in the amplitude of the lower frequency fluctuations (see red curves of Figures 6A,B). For M _y , we see that those oscillations follow the horizontal wake movement, therefore at a frequency that is close to that of the wake meandering. As R ₂ is fully aligned with R ₁, M _y remains centered around zero, and oscillates between positive and negative values, depending on the wind turbine side impacted by the wake. The tilt moment is also correlated to the horizontal wake movement, but the correlation is more obvious when we consider the vertical wake movement (see Figure 6D). So, even if the vertical wake movement is small, its impact remains important, mainly in a configuration where the rotor is most of the time impacted by the wake.

The behavior of M _f,1, M _e,1, M _y and M _t is different for α _w = 262°. As the wind turbine of R ₂ is misaligned with respect to that of R ₁, it is often fully out of the upstream wake, making its behavior more similar to a machine of row R ₁. For example, in Figure 7, between t ^* = 135 and 145, we see that the moments have oscillations with means and amplitudes that are similar to those of row R ₁: this can be linked to the wake centerline position, which is highly off-centered from the rotor hub at this particular time (see Figure 8A for t ^* = 136). However, the wind direction α _w = 262° still leads to situations where the wind turbine of row R ₂ is immersed in the upstream wake. For example, around t ^* = 200, the center of the upstream wake is closer to that of the rotor, resulting in an increase in the amplitude of the high frequency oscillations for all the moments. It is also highlighted in the streamwise velocities at t ^* = 201, reported in Figure 8A. As for α _w = 270°, there is an increase in the amplitude of the lower frequency oscillations, correlated to the horizontal motion of the upstream wake. For M _f,1 and M _e,1, those fluctuations occur at the same frequency as that of the wake meandering, as the wind turbine is off-centered from the upstream wake, and thus oscillates between periods where it is partially immersed in the wake and periods where it is fully out of the wake. The difference of the mean value around which the signal oscillates during these different periods can be large, as the rotor switches from a freestream situation, with high velocities, especially if it is impacted by a gust, to a waked configuration (for example between t ^* = 145 and 155, in Figures 7A,B).

For M _y , the mean of the signal varies between zero (when the rotor is out of the wake) and negative values, as the wake only impacts one side of the rotor. For M _t , the correlation between the low frequency oscillations and the vertical wake movement is less obvious, as the wind turbine is either out of the wake, or partially immersed in the upstream wake.

Row R ₆: For α _w = 262°, the behavior of row R ₆ is comparable to that of row R ₂, as the upstream wakes that impact these two rows have a similar behavior (see Figure 3A for the global view of the wind farm). However, the wake upstream of row R ₆ is slightly less structured, and wider: so, even when the wake centerline is estimated far from the rotor center, the rotor can be still impacted by some flow structures of the wake, as highlighted for t ^* = 126 and row R ₆ in Figure 8A. For α _w = 270°, compared to α _w = 262°, the behavior of the upstream wake differs more between R ₂ and R ₆. The amplitude of the wake meandering increases, leading to a wind turbine that is sometimes out of the upstream wake (for example around t ^* = 139, see Figure 8B). However, as the wake is more destructured and also wider, the rotor blade crosses wake structures most of the time, even when the center of wake is computed as being far from the rotor one. This decreases the amplitude of the moment fluctuations related to the blade passage and, mainly for the yaw moment, those of lower frequencies, correlated to the wake meandering. As the wake also recovers faster deeper in the wind farm, the velocities in the wake upstream of row R ₆ are higher, further attenuating the high frequency moment fluctuations occurring when the rotor is in half wake situation.

Row R ₁₀: For α _w = 270°, the behavior of the wind turbine of row R ₁₀ is globally similar to that of R ₆. For α _w = 262°, as already mentioned, in addition to the upstream wakes of its own line, R ₁₀ is also impacted by the wakes of the adjacent line L5 (see Figure 3 the global view of the wind farm). The rotor is then immersed in a more “uniform flow”, with changes in velocity speed that are less marked, as it is well noticeable in Figure 8A for t ^* = 126 and R ₁₀. This decreases the amplitude of the moment oscillations at the blade passage frequency, but also at lower frequencies, correlated to the large scale wake motion.

It is relevant to compare the LES results to those obtained with the actual measurements available in Hansen et al. (2014). Their work studied the flapwise root bending moment for the wind turbine located in R ₂ and L ₄ as a function of the different atmospheric stabilities and incoming wind speeds and directions. The wind speed was estimated based on the measurement at the nacelle of that wind turbine, while the wind direction was derived from pairs of undisturbed wind turbines assuming no yaw misalignment. The atmospheric stability was estimated based on the Monin-Obukhov theory and wind speed, air and see temperatures, all measured at met mast 7 (see Figure 1).

For wind directions that led it to freestream conditions and for 8 m/s, the wind turbine of R ₂ and L ₄ exhibits a measured time-averaged flapwise bending moment of slightly less than 1,400 kNm. These particular operating conditions can be related to those of R ₁ in our simulations, which corresponds to a freestream turbine. We obtain a time-averaged value that is equal to 1,250 kNm for this particular turbine, which underpredicts the mean M _f by about 10% compared to the measurements. The discrepancies between the LES results and the measurements may be explained in several manners. First, although a large part of the study of Hansen et al. (2014) was dedicated to the dependence of the fatigue loads to the atmospheric stability, it seems that the time-averaged value was computed from all the valid measurements, and not only from those representative of the neutral atmospheric stability. The atmospheric stability impacts the wind shear and the turbulence intensity, which may lead to differences in the resulting time-averaged M _f . Additional discrepancies may come from to the establishment of the inflow flow characteristics. As mentioned, the wind speed is measured at the nacelle of the wind turbine. For wind speeds below 11 m/s, the authors mention that the measurement at the nacelle matches well with the ambient flow when the wind turbine is in freestream condition. However, it is not specified if a correction accounting for the rotor induction, potential pressure gradients around the nacelle or wake effects (because the anemometer is located behind the blades) exists. This contributes to increase the uncertainties, which are not quantified in Hansen et al. (2014).

We also compare the time-averaged M _f for R ₂ and α _w = 270° to that obtained with the actual measurements of Hansen et al. (2014) for the same machine. For a mean wind speed of about 5.5 or 6 m/s, corresponding to the upstream velocity of R ₂ in our LES results, Hansen et al. (2014) obtained a value between 700 and 800 kNm for wind directions that set the wind turbine in fully aligned conditions. The LES results produce a time-averaged value of about 900 kNm, which corresponds to an error between 13 and 22%. This difference can be explained in several manners. First, the measurements were given for fully aligned conditions, which led to the lowest separation between wind turbines (i.e., 7D), but it was not explicitly mentioned in Hansen et al. (2014) whether these conditions accounted for α _w = 270° only. Other wind directions satisfy that criterion (e.g., α _w = 90°), and result in a different relative wind turbine position within the wind farm, and thus a different M _f behavior. Again, the uncertainties related to the estimation of the inflow characteristics may reinforce the discrepancies. The atmospheric stability also plays an important role, as it impacts the upstream wake behavior. We recall that no distinction was made between the stability types for the computation of the time-averaged value given in Hansen et al. (2014).

We now investigate the number of cycles, n _i , per moment amplitude M _i , computed using the Rainflow counting (RFC) method (ASTM A. E, 2003). We sill consider the wind directions α _w = 262° and 270° and the rows R ₁, R ₂, R ₆ and R ₁₀. In order to have consistent n _i in the large moment amplitudes, which are marginally captured by the simulations, we accumulate the cycles over six lines, from L ₂ to L ₇. We deliberately omit L ₁ and L ₈ in the computation, as they potentially have a different behavior, due to their position at the extremities of the wind farm. The resulting cycle distributions are reported on Figure 9, for the four moments. For the root bending moments, the cycles are moreover accumulated over the three blades: those moments are now denoted M _f and M _e . Again, the discussion follows a row-wise order.

Rotors of row R ₁: Globally, as expected, the cycle distribution for row R ₁ is similar for the two considered wind directions. The small discrepancies observed in n _i are due to the change of the rotor location according to the wind direction. We notice a peak in n _i for moments of middle amplitudes. As already mentioned, this results from the moment fluctuations when the blades crosse long-lasting areas of higher and lower velocities over their rotation.The cycles of large amplitude are due to more global changes in the incoming wind speed, i.e. larger gust or lulls.

Rotors of row R ₂: For the four investigated moments and the two wind directions, the distribution of n _i varies between R ₁ and R ₂. The number of load cycles of middle and large amplitudes increases, and the peak related to the blade rotation is shifted towards higher moment amplitudes, and spread over a wider range of moments. As already highlighted in the discussion about the instantaneous signals, it is due to the presence of an upstream wake, characterized by a velocity deficit that is clearly more marked than the ambient wind fluctuations.

Compared to α _w = 262°, the root bending moments obtained for α _w = 270° have a higher number of cycles around the middle amplitudes for M _f and M _e , with a peak that is slightly shifted towards higher amplitudes. This is due to the high frequency moment oscillations that have, for most of the simulation time, higher amplitudes than for α _w = 262°, as highlighted in Section 4.2. For α _w = 270° and M _e , the amplitude of the lower frequency oscillations, which are the same order of magnitude as those of higher frequency, contributes to increase the peak. For M _f in particular, α _w = 262° seems to lead to higher maximal amplitudes, although the number of those cycles remains small. It is due to the transitions between periods where the rotor is immersed in the wake, and periods where the rotor is out of the wake, as highlighted in Section 4.2. For α _w = 270°, M _y and M _t also present a non-zero number of cycle in the very large amplitudes, which may contribute to increase the differences in the resulting fatigue loads between the two different directions.

Rotors of row R ₆: For α _w = 262°, the n _i distribution for R ₆ is comparable to that observed for row R ₂, as the upstream wakes impacting these two rows globally have a similar behavior (see Figure 3A). For α _w = 270°, compared to row R ₂, R ₆ sees an upstream wake with a lower velocity deficit and that is wider, which reduces the amplitude of the 1P and 3P moment fluctuations when it is fully and partially immersed in wake and thus shifts the distribution towards lower moment amplitudes. The tail of the M _f cycle distribution also spreads over larger amplitudes: as the amplitude of the wake meandering upstream increases, the wind turbine located in R ₆ is sometimes out of the wake, increasing the amplitude of the low frequency oscillations when the it transits between waked and freestream situations. In contrast, the cycle distribution of M _e remains compact and does not extend to larger cycle amplitudes. Indeed, for M _e , compared to M _f , the amplitude of the low frequency fluctuations is closer to that of the 1P oscillations. In fact, it is even of the same order of magnitude as that of the 1P oscillations of row R ₂ (see Figure 6B). Consequently, the increase in the amplitude of the lower frequency oscillations between R ₂ and R ₆ does not contribute to spread the cycle distribution over large amplitudes. For M _y , if the wind turbine is more often completely out-of-wake, this will decrease the amplitude of the lower frequency oscillations due to an imbalance in the rotor loads (highlighted Figure 6C), moving the cycle distribution towards smaller amplitudes. For M _t , the tail of the distribution still extends to cycles with amplitudes similar to those of R ₂. However, Figure 6D highlights that, between R ₂ and R ₆, there is a net decrease in the frequency of the filtered moment fluctuations. Those fluctuations were correlated to the vertical wake movement for R ₂. As R ₆ is more often out of the wake, this is translated by a decrease in the number of cycles of about 400 or 500 kNm, and a decrease in the resulting fatigue damage (this will be discussed in Section 4.4).

Rotor of row R ₁₀: For α _w = 270°, R ₁₀ presents a similar cycle distribution as for row R ₆. For α _w = 262°, as R ₁₀ is impacted by the wakes of its own line, but also those of the adjacent line L ₅, the rotor is then immersed in a more “uniform flow”, decreasing number of cycles n _i in the large amplitudes.

In order to investigate the effects of the wind direction and the wind turbine location on the resulting fatigue damages, a quantitative fatigue indicator, the equivalent moment M _eq , is here computed. M _eq is estimated based on the cycle distribution computing using the RFC method, combined with a Palmgren-Miner rule Miner (1945). Assuming that the moment is proportional to the stress, M _eq is computed as M e q = ∑ i = 1 L c n i M i m N e q 1 m (10) where i and L _c represent the load case number and the total number of load cases, respectively; n _i is the number of cycles for the particular moment amplitude M _i , determined using the RFC method; m is a parameter dependent on the material; N _eq is an arbitrary number of cycles. As we only investigate M _eq as ratios, N _eq does not need to be defined in this work. In the present study, M _i corresponds to a particular M _f , M _e , M _y or M _t amplitude. The exponent m is here defined as 10 for M _f and M _e , which is a widely used value for blades made of fiberglass (Mandell and Samborsky, 1997), such as those of the Vestas 2 MW. The yawing and tilting moments impact, among other components, the main shaft, which is usually made in forged steel (AWEA, BlueGreen, GLWN, and NIST, 2011): m is chosen to be equal to 4, which is within the range between 3 and 5 typically used for the metallic wind turbine components (Guideline and Lloyd, 2010). In each line, the fatigue estimate is normalized by M _eq of the first row. We still accumulate the cycles from line L ₂ to line L ₇, and also over the three blades for M _f and M _e .

We organize the discussion around these four equivalent moments shown in Figure 10.

Flapwise moment M _eq,f : Globally, for the two wind directions, there is a marked increase in M _eq,f between the first row and the downstream rows. This is an expected behavior after the discussions of Sections 4.2 and 4.3. Compared to α _w = 270°, α _w = 262° leads to a higher increase in the equivalent M _f for rows R ₂ to R ₆, despite its smaller n _i in the small and middle amplitudes (see Figure 9A). Indeed, for this wind direction, the number of cycles in the very large moment amplitudes is non-zero. Even if their numbers are small, those cycles are particularly significant in the fatigue estimates, notably due to the presence of the m exponent in the computation of M _eq . From row R ₇ for α _w = 262°, the rotors are additionally impacted by the wakes of the adjacent lines, slightly reducing the number of large amplitude cycles, and shifting the peak related to the blade rotation towards lower amplitudes (see discussion of Section 4.3). As a consequence, the resulting fatigue damage decreases, and becomes smaller than that obtained for α _w = 270°. For α _w = 270°, although the cycle distribution is slightly shifted towards lower amplitudes for rows downstream of R ₂, it is also spread over larger amplitudes (as highlighted in Figure 9A). As a consequence, the ratios M e q , f , R i / M e q , f , R 1 still increase after row R ₂.

Edgewise moment M _eq,e: The increase in M _eq,e is similar in both wind directions, excepted for rows R ₇ to R ₁₀, where the fatigue estimates is lower for α _w = 262°. For both wind directions, R ₂ presents the highest increase in M _eq . Indeed, for rows located downstream, the cycle distribution is shifted towards lower amplitudes, without an increase in the number of cycles of large amplitude, as highlighted in Figure 9B for R ₆ and R ₁₀ and in the discussion of Section 4.3.

Yaw and tilt moments, M _eq,y and M _eq,t: For α _w = 270°, the increase in the fatigue damage due to M _y and M _t is more important than that obtained for α _w = 262°. It was highlighted in the discussion of the instantaneous signals in Section 4.2. The fatigue damages also slightly decrease for rows located deeper in the wind farms. Again, this was highlighted in Figure 9 and in the discussion of Section 4.3. For α _w = 262°, the equivalent moments slightly increase between row R ₃ and R ₆, and, again, decrease for the last rows.

It is relevant to compare the increase in fatigue equivalent loads to results available in the literature. To that end, the ratio of flapwise equivalent loads between R ₁ and R ₂ as a function of the incoming wind direction is reported in Figure 11. We here consider all the simulated wind directions; in order to make clearer the comparison with other results, the wind direction is expressed as the offset with respect to α _w = 270°. We also report on Figure 11 the results of Lee et al. (2018). Lee et al. (2018) numerically investigated the damage equivalent loads in a pair of floating wind turbines. The first wind turbine and its resulting wake were simulated using LES coupled to an AL model. The behavior of the second wind turbine was modeled by extracting time series of a planar data at 7D downstream from the turbine to use as an input to FAST. They made vary the lateral position of the downstream wind turbine, in order to simulate fully and partially immersed situations. They consider a turbine of 3MW, and thus with a radius larger than the Vestas 2MW, but, as mentioned, the spacing between the two wind turbines is 7D, which is equal to our configuration for α _w = 270°. Moreover, the wind speed is about 9 m/s, meaning that the rotor operates in the second region of the power curve, with a thrust coefficient close to that of our LES results. The TI is set to less than 4%, which is lower to the TI used in the Horns Rev simulations. This may lead to discrepancies between the two sets of results, as the TI affects the wake meandering. The comparison between the results should be also cautiously done, as Lee et al. (2018) consider floating wind turbines, capturing the interactions between the platform, the wake and the turbine. However, we expect a similar behavior in terms of the evolution of the fatigue damage.

Globally, we see the same tendency for both sets of results, with a lower increase of M _eq for fully aligned conditions than for small offsets. This tendency is also observed in other research works (Schmidt et al., 2011). Indeed, when the wind turbine is fully aligned with the upstream machine, there are periods where the rotor remains fully immersed in the upstream wake, decreasing the flow variations that the blade crosses over one rotation. The amplitude of the load cycles associated to those 1P oscillations is thus smaller, which decreases the resulting damage equivalent loads. For small offset (e.g., + or -5°), as the wind turbine is most of time partially immersed in the upstream wake, the increase in fatigue loads is more significant. When the offset increases, the rotor is more often out of the wake (this was well highlighted in the discussion of Sections 4.2 and 4.3), which decreases the damage equivalent loads. However, for an offset of about 8°, we see in Figure 11 that there are more discrepancies between the two sets of results. Lee et al. (2018) explained their low increase in M _eq for that offset by a wake deficit region that was non-symmetric due to its interaction with the vertical shear and the presence of the Coriolis effect. We do not account for the Coriolis forces in our LES, which may explain that the wake remains more symmetric, and that we do not observe a net decrease for the positive offset of about 8°. However, for the other offsets, the order of magnitude of the increase of M _eq is globally similar in the two sets of results.