the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# Near-wake analysis of actuator line method immersed in turbulent flow using large-eddy simulations

### Jörn Nathan

### Christian Masson

### Louis Dufresne

The interaction between wind turbines through their wakes is an important aspect of the conception and operation of a wind farm. Wakes are characterized by an elevated turbulence level and a noticeable velocity deficit, which causes a decrease in energy output and fatigue on downstream turbines. In order to gain a better understanding of this phenomenon this work uses large-eddy simulations together with an actuator line model and different ambient turbulence imposed as boundary conditions. This is achieved by using the Simulator fOr Wind Farm Applications (SOWFA) framework from the National Renewable Energy Laboratory (NREL) (USA), which is first validated against another popular Computational Fluid Dynamics (CFD) framework for wind energy, EllipSys3D, and then verified against the experimental results from the Model Experiment in Controlled Conditions (MEXICO) and New Model Experiment in Controlled Conditions (NEW MEXICO) wind tunnel experiments. By using the predicted torque as a global indicator, the optimal width of the distribution kernel for the actuator line is determined for different grid resolutions. Then, the rotor is immersed in homogeneous isotropic turbulence and a shear layer turbulence with different turbulence intensities, allowing us to determine how far downstream the effect of the distinct blades is discernible. This can be used as an indicator of the extents of the near wake for different flow conditions.

- Article
(12207 KB) - Full-text XML
- BibTeX
- EndNote

An important aspect for the conception of wind farms is the turbine spacing, which depends on the interaction of wind turbines through their wakes. This phenomenon can decrease the wind park energy output by up to 20 % due to the velocity deficit propagated by the wakes (Manwell et al., 2010). Additionally, it can increase the turbine fatigue due to the increased turbulence intensity. In order to study wake interactions, the flow around the rotor has to be modelled correctly. Hence, the model should account for the apparition of turbulent structures of different magnitudes, for instance, the vortices created by the blade tips and their interaction with the ambient turbulence.

As opposed to the far-wake region (Olivares Espinosa, 2017), the near-wake representation in a computational fluid dynamics simulation depends heavily on the applied rotor model. Approaches range from an actuator force representation inserted as momentum sink in the Navier–Stokes equations to full rotor modelling where the attached boundary layers on the blades are simulated (Sanderse et al., 2011). This work will apply the actuator line method (ALM) in order to model the transient behaviour of the rotor by representing distinctly the rotating blades as presented by Troldborg (2009). Each blade is represented by a force line allowing us to reproduce the helicoidal vortical structure in the near wake and allowing us to assess its interaction with the flow.

In order to evaluate the soundness of the present method, a comparative study of the Simulator fOr Wind Farm Applications (SOWFA) framework, from the National Renewable Energy Laboratory (NREL), and EllipSys3D, from DTU, was conducted as initially presented in Nathan et al. (2017). Based on this study, the method used throughout this work will be evaluated before proceeding to establish the base case for the non-turbulent inflow. For establishing the base case, the optimal width of the distribution kernel of the forces of the actuator line is determined. While previous work often focused on numerical stability as in Troldborg (2009) or Ivanell et al. (2010), when choosing the distribution width, Martínez-Tossas et al. (2015) states that, with decreasing distribution width, the line forces are getting too concentrated, resulting in a wrong prediction of the rotor torque. Hence, this work tries to evaluate the optimal width for each mesh resolution by using the predicted torque as a global indicator.

For the introduction of a turbulent inflow, different methods exist for imposing a statistically generated velocity field, such as inserting it via a momentum sink as done in Troldborg et al. (2011) or as boundary conditions as done in Olivares Espinosa (2017). This work adheres to the latter approach, as it was seen as more straightforward than the conversion of the velocity field to a force which then is translated back to the velocity field by the numerical solver.

Then, shear layer turbulence (Muller et al., 2014) is introduced, exposing the rotor model to a more realistic wind flow situation bearing more resemblance to applied wind energy. This novel approach takes into consideration the temporal evolution of the sheared velocity field, hence allowing it to be imposed as boundary condition as well. Also, this method was not yet applied to the actuator line method.

Finally, the numerical results are used to examine the spatial extents of the near-wake region. While in previous work such as Krogstad and Eriksen (2013) or Sarmast et al. (2016) often the profiles of velocity deficits or turbulent kinetic energy are taken into consideration for evaluating the near wake, this work uses the energy spectra to determine how far downstream the discernible effects of the distinct blades are noticeable. While the analysis of the turbulent inflow in previous work (Olivares Espinosa, 2017) has often been conducted using energy spectra, energy spectra including the rotor effects are seldom included in the analysis of the near wake. As with increasing turbulence intensity the statistical convergence tends to be longer, the energy spectra approach in this work permits the analysis of the spatial extensions of the near-wake region even without full convergence of second-order statistics.

The numerical simulations are based on the incompressible Navier–Stokes equations:

with **F** representing the actuator force inserted
as a momentum sink, **U** the velocity, ∇*p*∕*ρ* the modified
pressure and *ν* the kinematic viscosity. In Sect. 2.1 the
rotor model and the derivation of the force term **F** are discussed. Then, in Sect. 2.2 an overview is given on the methods
generating the two different ambient turbulences and how they are imposed as
boundary conditions. Finally, a summary of the numerical framework and its
setup is given in Sect. 2.3.

## 2.1 Rotor model

The force term **F** in Eq. (1) is obtained using

with the lift and drag forces shown in Fig. 1, and defined as

with *c*_{l} and *c*_{d} being the lift and drag coefficient,
*U*_{mag} the sampled velocity magnitude in the blade reference frame,
*c* the chord width, *l*_{s} the length of the actuator segment,
𝒢 the Gaussian kernel and *f*_{tip} the tip correction.

The same unmodified airfoil coefficients were used as shown in
Shen et al. (2012). These data were obtained from 2-D experiments without
rotation in a wind tunnel, and therefore they do not include the stall delay
due to boundary layer stabilizing effects such as Coriolis and centrifugal
forcing, which enhance the lift of the airfoil (Vermeer et al., 2003) near the
root of the blade. As shown in Nathan et al. (2017) using unmodified airfoil
data results in problems predicting the blade forces in the root region for
high wind speed flows. They seem to handle the moderate wind
speeds which are in the focus of this work fairly well. Hence, despite their shortcomings,
the unmodified airfoil data are used throughout this work. Then, the forces in
Eqs. (4) and (5) are projected onto the blade reference
frame using the unit vectors *e*_{L} and *e*_{D} in
Eq. (3).

While the Glauert tip correction *f*_{tip} was originally intended
(Glauert, 1935) to represent the otherwise absent tip vortices in the
actuator disk model, it still proves advantageous for the ALM at lower
resolutions as shown in Nathan (2018). Due to the relatively low
resolution, the shed tip vortices from an ALM are much larger than the ones
observed experimentally. Hence, the induction caused by the simulated vortices
is weaker than in reality, and the Glauert tip correction permits us to
compensate for it (Nathan et al., 2017).

Finally, in order to avoid spurious oscillations around the point of the
inserted force, the punctual force is distributed using a kernel function
$\mathcal{G}\cdot (\cdot )$. As done in previous works such as
Troldborg (2009) or Olivares Espinosa (2017) this work adheres to a
normal distribution with the distribution width *ϵ*=*σ*Δ*x*, with *σ* being the distribution width of the normal function and Δ*x*
the cell width.

## 2.2 Turbulence inflow generation

### 2.2.1 Homogeneous isotropic turbulence

A synthetic velocity field representing homogeneous isotropic turbulence based on the von Kármán energy spectrum (Pope, 2000)

is obtained by using the algorithm proposed by Mann (1998). The
technical details can be found in the article of Mann (1998) or more
recently in Olivares Espinosa (2017). The main parameters for this approach are
the integral length scale *L* and the coefficient *α**ϵ*^{2∕3}, which can be used as a scaling factor to obtain the desired amplitude of the
turbulent structures. The range of wavenumber *κ* depends on the grid
resolution and dimension extents. Hence, these parameters determine the
ability of the numerical mesh to resolve a certain range of turbulent scales.

While several implementations of this method exist, e.g. Olivares Espinosa (2017) or Muller et al. (2014), the implementation of Gilling (2009) was chosen to synthesize the homogeneous isotropic turbulence (HIT) for several reasons. It is in public domain, it corrects for a divergence-free velocity field and allows us to impose HIT at the boundaries at relatively low computational cost.

In Fig. 2 the midplane of a generated turbulent field is shown for
different turbulence intensities. The flow structures are identical apart
from the different scaling of the velocity fluctuations. This results from
using the same seed for the random number generator in the Mann algorithm and
by scaling the obtained velocity field with *α**ϵ*^{2∕3} to obtain
the desired synthetic turbulence intensity TI_{syn}.

Contrary to Troldborg (2009) where the synthetic turbulence is
imposed as a momentum source, in this work it is imposed as a boundary
condition as done in Muller et al. (2014) and Olivares Espinosa (2017). The
velocities are imposed by convecting the velocity field of the synthetic
turbulence through the computational domain by the mean velocity *U*_{∞}
at each simulation time step. They are then projected by trilinear
interpolation onto the computational points. In order to speed up the
statistical convergence, the simulation is also initialized with the
synthetic turbulence field.

### 2.2.2 Shear layer turbulence

Based on the Mann algorithm (Mann, 1998), Muller et al. (2014) developed a
method to impose the synthetic turbulence on a sheared flow as boundary
condition, including the evolution of the vortical structures. Apart from the
turbulence intensity, no further turbulence characteristics were published by
Schepers et al. (2012). Hence, a turbulence length scale of *R*∕4 was
chosen. These scales were larger than the width of the coarsest cells in the
numerical mesh in order to dampen the effect of numerical dissipation in
axial direction. A typical flow field generated by this algorithm can be seen
in Figs. 3 and 4.

The mean velocity profile is obtained via the power law

hence, the velocity at the bottom of the domain has not necessarily got to be
zero. Therefore, the computational mesh can be much smaller than in a
wall-resolved flow, as its mesh has to include the ground and has to have a
high refinement in this region. The reference height *z*_{ref} was set to hub height, and the reference velocity was set to *U*_{ref}=15 m s^{−1}. The parameter *α*_{p}
can usually be deduced from experimental measurements if available. As they
were not available for this experiment, standard conditions are assumed with
${\mathit{\alpha}}_{p}=\mathrm{1}/\mathrm{7}$ (Pope, 2000).

## 2.3 Numerical framework

This work is realized within the open-source framework
OpenFOAM (version 2.2.2)^{1} together with the
SOWFA project,^{2} which contains a similar implementation of the ALM as presented by
Troldborg (2009). A more detailed explanation of the implementation
can be found in Martínez-Tossas et al. (2016). OpenFOAM is a set of libraries
and executables entirely written in C++. While the first released scientific
article about the framework was by Weller et al. (1998), its inner workings are
described in more detail by Jasak (1996).

The computational domain is cubic with an edge length of ${L}_{x}={L}_{y}={L}_{z}=\mathrm{10}\phantom{\rule{0.125em}{0ex}}R$, with *R* being the rotor radius and the rotor positioned at the domain
centre. The cells in the rotor vicinity are refined in the range of $-\mathrm{0.4}\le x/R\le \mathrm{0.4}$ with the size $\mathrm{\Delta}x=D/\mathrm{128}$. Within SOWFA, several
refinement zones are applied, each time halving the cell edge length as also
done in Vanella et al. (2008). The final mesh size consists of
1.9 × 10^{6} cells. The technique used in SOWFA proves highly
advantageous in terms of computational cost. As the mesh dimensions at first
sight seem relatively small compared to other work such as
Troldborg (2009) or Martínez-Tossas et al. (2017), an extensive
sensitivity study was conducted by varying domain extents in axial direction
up- and downstream of the rotor as well as in the lateral direction for
different grid refinements. The findings show that the dimensions used have a
negligible impact as shown in Nathan (2018).

For the boundary conditions, the velocity is imposed as a uniform inflow velocity of $\mathbf{U}=({U}_{\mathrm{\infty}},\phantom{\rule{0.125em}{0ex}}\mathrm{0},\phantom{\rule{0.125em}{0ex}}\mathrm{0})$ for the non-turbulent flow and in the turbulent cases as the synthetic velocity as explained in Sect. 2.2. The lateral boundaries are set as symmetric for the non-turbulent and homogeneous isotropic turbulence case. For the shear layer turbulence, the velocity is also imposed at the lateral boundaries.

The large eddy simulations use the dynamic Lagrangian sub-grid-scale model (Meneveau et al., 1996). For the discretization of the convective term, a linear combination of 75 % central differencing and 25 % second-order upwind scheme is applied as presented by Warming and Beam (1976). In OpenFOAM terminology, this scheme is called linear-upwind stabilized transport (LUST). The choice of the scheme is made as a trade-off between the accuracy of a linear discretization and the stability of an up-winding scheme. This scheme proved to preserve the turbulent structures well (Nathan, 2018). The remaining spatial terms are discretized by central differencing, and for the time discretization the Crank–Nicolson method is used.

The pressure is resolved using a geometric agglomerated algebraic multi-grid solver and the remaining variables are solved for with a bi-conjugate gradient method using a diagonally based incomplete LU preconditioner. The total simulation run-time comprises 60 rotor revolutions (≈8.5 s), and the time step has to be small enough to avoid the actuator point representing the blade tip skipping a computational cell during rotation. It is set to $\mathrm{0.327}\times {\mathrm{10}}^{-\mathrm{3}}$ s. The total run-time is chosen as first- and second-order statistics are deemed to be converged.

For the parameterization of the ALM, different distribution widths are chosen in order to obtain the optimum for the examined case and 40 actuator points are used to represent one blade in accordance with what was found in Nathan (2018).

## 3.1 Validation and verification

The implementation was validated against EllipSys3D and verified against the Model Experiment in Controlled Conditions (MEXICO)
and New Model Experiment in Controlled Conditions (NEW MEXICO) experiment in Nathan et al. (2017). The MEXICO rotor is
a three-bladed rotor with a radius *R*=2.25 m. It rotates at a constant rpm of
424.5 min^{−1} for the three different flow cases with ${U}_{\mathrm{\infty}}=\mathit{\{}\mathrm{10},\mathrm{15},\mathrm{24}\mathit{\}}$ m s^{−1}; hence, their respective tip speed ratios are
10.0, 6.7 and 4.2. As shown in Nathan (2018), the chord-based
Reynolds number $Re={U}_{\mathrm{\infty}}c/\mathit{\nu}$ varies from 0.4×10^{6} for
the low-velocity case towards the hub up to 0.6×10^{6} for the high-velocity case in the tip region. The power coefficients for the observed
cases are ${C}_{\mathrm{P}}=\mathit{\{}\mathrm{0.28},\mathrm{0.36},\mathrm{0.19}\mathit{\}}$. A comparison of the axial
profiles of the velocity components can be seen (Fig. 5) showing that
both codes reproduce very similar results in the near wake further away from
the rotor. Apart for the high-velocity case where 3-D effects become
important, an excellent agreement can be observed for the other two cases.

For the high-velocity case (*U*_{∞}=24 m s^{−1}), the vortex sheets
shed from the blades become visible by the oscillations in the axial velocity
component *U*_{x}.

## 3.2 Non-turbulent flow

When refining the grid using the actuator line method, the distribution
parameter *ϵ* has to be adjusted to obtain a global torque 〈*T*〉 close to the reference value *T*_{ref}. In the following, only
the case for *U*_{∞}=15 m s^{−1} will be examined as the other
cases in Nathan et al. (2017) served as extreme cases for determining how the
model behaves at its limits.

Instead of relying on a constant *ϵ*∕Δ*x* for different grid
resolutions, this work adapts *ϵ*∕Δ*x* depending on the grid
resolution or number of cells across the rotor diameter $N=\mathrm{2}R/\mathrm{\Delta}x$.
The results are shown in Eq. (6). A confidence interval of ±1 % was established around the reference torque value *T*_{ref}.
Through iterations, an optimal distribution parameter is found to fall in
this range.

The lower bound for the distribution parameter here is *ϵ*=1.7Δ*x* for the sake of the numerical stability of the applied method. Other
frameworks applying a different numerical discretization can go even lower, e.g. in Ivanell et al. (2010). By doing so, it can be seen that the best
solution in terms of global torque for a resolution of $N=D/\mathrm{\Delta}x=\mathrm{32}$
is off by around 4 % in Fig. 6.

As a general trend, it can be seen that *ϵ*∕Δ*x* has to be
increased with increasing resolution. This stems from the fact that by
refining the mesh with a constant *ϵ*∕Δ*x*, the punctual induction
caused by the blade would be too high and eventually the torque would be
below the reference value, e.g. for *ϵ*=2Δ*x* for *N*≥64. By contrast, when having a very low resolution, a constant *ϵ*
distributes the force too widely, causing a lower induction around the rotor
resulting in an overestimation of the torque, e.g. for *ϵ*=2Δ*x* for *N*≤48.

The optimal distribution parameter *ϵ*∕Δ*x* found in
Fig. 6 is now shown in dependence on the grid resolution $N=D/\mathrm{\Delta}x$ in Fig. 7. It seems as if this value would eventually reach an asymptotic limit for higher resolutions. Looking at a more
theoretical approach in Martínez-Tossas et al. (2017) suggests that the
optimal distribution width *ϵ* lies between 0.14 and 0.25 of the chord
*c* whereas in this case for $D/\mathrm{\Delta}x=\mathrm{128}$ the *ϵ*∕*c* lies between
0.5–8.9 depending on the span-wise location. The observation made by
Martínez-Tossas et al. (2017) is backed by Shives and Crawford (2013), where
*ϵ*∕*c* falls in the same range. But it should be kept in mind that
Shives and Crawford (2013) use a much higher grid resolution, allowing
$\mathit{\u03f5}/\mathrm{\Delta}x\ge \mathrm{4}$ and in the case of Martínez-Tossas et al. (2017)
even $\mathit{\u03f5}/\mathrm{\Delta}x\ge \mathrm{5}$.

The curvature in Fig. 7 also confirms the findings of Jha et al. (2013)
that keeping the relation $\mathit{\u03f5}/\mathrm{\Delta}x=\mathrm{const}$ while
increasing the resolution is not a very good solution. While
Ivanell et al. (2010) suggests choosing the smallest possible distribution
width *ϵ* in order to minimize interactions with the vortical
structures ($\mathit{\u03f5}/\mathrm{\Delta}x=\mathrm{1}$), the observations made here fall more in
line with work such as Shives and Crawford (2013), suggesting *ϵ*∕Δ*x* to
be adapted to the physical model in order to distribute the force over a
meaningful length scale.

For an excerpt of the resolutions presented in Fig. 6, the radial
profiles of the velocity components can be found in Fig. 8. It can be
seen that the method seems to converge towards a solution when refining the
mesh. As shown in Fig. 6 the lowest resolution at *N*=32
over-predicts the torque by distributing the force to widely, which is also
reflected in the low axial induction downstream at $x/R=\mathrm{0.13}$. Despite
following the trend of the experimental values well, the method seems to
converge towards radial profiles which are especially off in the tip and hub
region where the strongest vortices are shed. These are limitations intrinsic
to the ALM, which is less apparent when using high-fidelity approaches such as
full rotor simulations (Carrión et al., 2015). In order to ameliorate the
results at the tip, a non-isotropic kernel could be investigated
(Rullaud et al., 2018).

In Fig. 9 the shed vortical structures can be seen in dependence on the grid resolution. While the root vortex is rather diffuse, a clear tip vortex can be noticed. It is interesting to note the vortices shed around mid-span due to the suboptimal choice of the airfoils of the blade causing a sudden change in circulation.

In order to estimate the resolution necessary to obtain tip vortex radii as
seen in the MEXICO experiment, the vortex radii are shown in Fig. 10.
The vortex radius *r*_{core} is defined as the limit containing 99 %
of the circulation. The Gaussian distribution is used as an approximation for
the vorticity distribution within the vortex. This assumption is normally
applied for low Reynolds number flows, while this case exhibits a Reynolds
number of $R{e}_{\mathrm{\Gamma}}=\mathrm{\Gamma}/\mathit{\nu}=\mathrm{220}\phantom{\rule{0.125em}{0ex}}\times \phantom{\rule{0.125em}{0ex}}{\mathrm{10}}^{\mathrm{3}}$. As these findings
are related to the vortex dynamics of the flow *R**e*_{Γ} is used instead
of *R**e* defined earlier. Nevertheless, this approximation is used in order to
be able to draw an analogy between the experimental and the numerical
results. It holds fairly well when comparing the Gaussian distribution and
the vorticity for *N*=128 as shown in Nathan (2018). Hence, by assuming
a Gaussian distribution for the vortices in the MEXICO experiment, a
corresponding distribution parameter *ϵ* can be deduced as shown in
Fig. 10.

This would necessitate a resolution of *N*≥4096 for the presented case, which would result in a computational grid beyond any justifiable
computational scope. Full rotor calculations as conducted by
Carrión et al. (2015) allowed us to obtain tip vortices of ${r}_{\mathrm{core}}/R\approx \mathrm{0.012}$ for *N*≈900 in the tip region, which corresponds very
well to results in Fig. 10. Another result for the vortex radius can
be found in Nilsson et al. (2015), where for $\mathit{\u03f5}/\mathrm{\Delta}x=\mathrm{1}$ and *N*≈244 in the tip region a vortex core radius of ${r}_{\mathrm{core}}/R\approx \mathrm{0.055}$ was found. Despite the radii in this work and the references
being calculated based on three different methods, the results fall within the
same range.

## 3.3 Homogeneous isotropic turbulence

In Fig. 11 the longitudinal evolution of the turbulence intensities
can be seen. There is a stronger decay for higher turbulence intensities, which was also found in Olivares Espinosa (2017). In that work EllipSys3D was
compared to a solution based on OpenFOAM, and it was found that over the same
longitudinal distance of 10*R* an absolute difference in the turbulence
intensity of 48 % and 44 % occurred for each framework respectively. This
stands in a stark contrast to the 4 % in this case for the high turbulence-intensity case. This huge decay, which is even more significant for
EllipSys3D, necessitates approaching the introduction of the turbulence close
to the turbine for high turbulence-intensity cases (Olivares Espinosa, 2017).

An important aspect when imposing a synthetic turbulence as boundary
conditions of a Computational Fluid Dynamics (CFD) simulation is respecting the Nyquist–Shannon sampling
theorem (Shannon, 1949) as also mentioned by Muller et al. (2014). Hence, a
study considering different ratios between the grid resolution of the
synthetic turbulence and the simulation was undertaken. It is found that the
higher the computational resolution is compared to the one of the synthetic
turbulence, the less the turbulence intensity decays in longitudinal
direction. While the criterion of Nyquist–Shannon states that the resolution
of the computational domain should be d$x/\mathrm{\Delta}x>\mathrm{2}$, this work uses the
ratio of d$x/\mathrm{\Delta}x=\mathrm{2.5}$, with d*x* being the cell width of the synthetic field
and Δ*x* the cell width of the computational mesh.

When taking the case for TI_{syn}=5 %, it is interesting to
note that while the resolved TI (green dashed line) is around 4.2 % at
the rotor position $x/R=\mathrm{0}$, a huge part of the difference in relation to the
imposed turbulence falls in the SGS model with $\langle {\mathrm{TI}}_{\mathrm{res}}\rangle +\langle {\mathrm{TI}}_{\mathrm{sgs}}\rangle =\mathrm{4.8}$ % and finally just a
relatively small amount of the turbulent intensity or turbulent kinetic
energy is “lost” by numerical dissipation.

Despite the fact that the computational grid respects the Nyquist–Shannon criterion for signal sampling in respect to the synthetic grid, immediately at the inlet a part of the turbulence falls in the sub-grid range. Due to the numerical dissipation caused by the differencing schemes and turbulence modelling, the energy cascade hands down its energy to lesser scales than the resolved ones.

It should be kept in mind, that the turbulence intensity the rotor model is
experiencing through velocity sampling is the resolved turbulence intensity
TI_{res}, and the sub-grid turbulence intensity TI_{sgs} is
therefore only felt indirectly, by an augmentation of the effective
viscosity. When looking at the fraction of the resolved turbulent kinetic
energy over the total turbulent kinetic energy, it can be seen that the
resolved scales exceed 96 %, which lies well above the criterion of 80 %
proposed by Pope (2004). For the flow case with *U*_{∞}=15 m s^{−1}, the total simulation run-time *t*_{tot} results in
roughly ${t}_{\mathrm{tot}}\cdot {U}_{\mathrm{\infty}}/\left(\mathrm{10}R\right)\approx \mathrm{5.5}$ flow through
times. The synthetic turbulence field is large enough that it does not need
to be recycled during one simulation.

In Fig. 12 the effects of the ambient turbulence on the turbine wake
are shown. While there are no noticeable impacts for the low-turbulence case
with TI_{syn}=0.1 %, the beginning of strong non-linear
interactions can be observed for TI_{syn}>0.1 *%*. For
TI_{syn}=15 % the inflow turbulent structures seem to outgrow
the structures created by the wind turbine.

In Fig. 13 it can be seen how the strength of the vortical structures
of the ambient fluid increases with higher turbulence intensity up to the
point for *T*_{syn}=15 % where its amplitude almost equals the one
emitted by the rotor model. In Fig. 14 the impact of the rotor
presence on the energy spectrum can be seen. The wavenumber
*κ*_{p} relating to the frequency of a blade passage (3 times
rotor frequency) obtained by ${\mathit{\kappa}}_{\mathrm{p}}=\mathrm{2}\mathit{\pi}\mathrm{3}f/{U}_{\mathrm{\infty}}$
shows a very distinct peak and its higher harmonics at the multiples of
*κ*_{p}.

As the velocity time series obtained from the simulations do not exhibit periodicity, the Welch method (Welch, 1967) is used to generate the energy spectra.

It is interesting to note that the distinct peaks in the spectra occur at the
wavenumber relating to the frequency of the blade passage and its harmonics.
The harmonics are caused by the strong excitement of the fluid by the blade
passage and its interaction with the non-linear term in the NS equations. As
the blade forces and hence the strength of the tip vortices are very
comparable, the peaks are very similar among the different cases for $-\mathrm{0.4}\le x/R\le \mathrm{0.2}$. The higher the turbulent kinetic energy content stemming
from the ambient flow, the faster the peaks are dampened and blend into the
ambient flow. For example there is almost no discernible effect by the blade
at $x/R=\mathrm{0.4}$ for TI_{syn}=15 % while for TI_{syn}=0.1 % the velocity oscillations are still very noticeable. Although it is
of lesser amplitude, the upstream region is also under the influence of the
distinct blades to a certain extent.

## 3.4 Shear layer turbulence

When looking at the vertical planes in Fig. 15, the influence of the sheared flow can be seen by a higher velocity deficit in the wake on the lower half of the rotor. In the vertical plane in Fig. 16, it can be seen that there is an increase in the vorticity magnitude towards the ground. While the increase appears to be rather subtle, it is shown in Nathan (2018) that the TI increases significantly towards the ground as expected in a shear layer flow. Looking at the energy spectra in Fig. 17 reveals a similar picture as shown above for the case of homogeneous isotropic turbulence. As the blade forces and hence the strength of the tip vortices are very comparable, the peaks are very similar among the different cases for $-\mathrm{0.4}\le x/R\le \mathrm{0.2}$. Due to the dissipation caused by the ambient turbulence, these peaks dampen at a different pace as seen at $x/R=\mathrm{0.4}$.

While before and at the rotor position for $-\mathrm{0.4}\le x/R\le \mathrm{0}$ the peaks remain very distinct, vortical structures by the ambient fluid and emitted by the blade cause the injected peaks to dampen and distribute energy to adjacent wavenumbers as seen clearly for $x/R\ge \mathrm{0.2}$. Depending on the level of the ambient turbulence the peak gets attenuated up to a point where it almost completely blends in with ambient turbulence as seen for TI${}_{x,z/R=\mathrm{0}}=\mathrm{15}$ % at $x/R=\mathrm{0.4}$. This relates to the observation made earlier that ambient structures are almost as important as the structures emitted by the blade.

This means that in the near wake in a turbulent flow with an ambient turbulence intensity of TI${}_{x,z/R=\mathrm{0}}\ge \mathrm{10}$ %, the velocity fluctuations at $x/R=\mathrm{0.4}$ already seem to have only a weak relation to the injected turbulence by the rotor but a much stronger one to the ambient turbulence. This means that for this kind of flow, an actuator disk method would probably also be sufficient when looking at the flow characteristics beyond $x/R=\mathrm{0.4}$.

By using a validated actuator line implementation (Nathan et al., 2017), it was
shown that the distribution width *ϵ*∕Δ*x* has a non-linear
dependence on the grid resolution and probably converges towards values
suggested in Martínez-Tossas et al. (2016). The rotor torque is used as a global
indicator for determining the distribution width, but the rotor thrust
followed the same trend. Hence, it is interesting to see that while the rotor
induction is predicted well, the velocity deficit agrees well only for $x/R>\mathrm{5}$ but not in the ultimate rotor vicinity.

It is also shown that with increasing grid resolution the spatial profiles seem to converge. This would be one aspect of a grid-independent solution, but it is still very far away from resolving the shed tip vortices correctly. Although it seems to converge towards a value of $\mathit{\u03f5}/\mathrm{\Delta}x\approx \mathrm{4}-\mathrm{5}$ for which the dimensions of the experimental vortices would be attained, this causes excessive computational costs due to the large mesh.

When looking at the turbulent inflow, a synthetic turbulence generated by the Mann algorithm (Mann, 1998), it was shown that the decay of the turbulence intensity in longitudinal direction is much less pronounced than in previous work. As shown for the axial decay of the turbulence intensity, a significant part of the difference between the resolved turbulence intensity and the imposed one from the synthetic field resides within the sub-grid scales. Hence, there is very little loss due to numerical dissipation, which is also reflected in the energy spectra, which are the better the higher the turbulent content is.

As expected the wake does recover at a faster pace for a higher turbulence
intensity. It is very interesting to note that the turbulent structures of
the ambient flow eventually catch up with the amplitude of the structures
emitted by the rotor. This is already noticeable in the instantaneous
velocity fields but becomes even clearer when evaluating the spectra. When
considering the velocity fluctuations in the downstream flow caused by the
blade passages for determining the near wake, it can be observed that in this
case for TI_{syn}≥10 %, the near wake already ends at $x/R=\mathrm{0.4}$. This is particularly interesting as a turbulence intensity of 15 % at
hub height is still considered to be low turbulence intensity according to
ISO 61400 (2015), and many real sites exhibit even higher turbulence
intensities. Hence, for some cases, the limit of the near wake would be $x/R=\mathrm{0.4}$ and even lower.

The SOWFA framework on which this work is based is made available by NREL https://github.com/NREL/SOWFA/ (last access: 1 March 2018) and the turbulence generator for the homogeneous isotropic turbulence can be obtained via http://vbn.aau.dk/en/publications/tugen(3e097a90-b3d8-11de-a179-000ea68e967b).html (last access: 1 March 2018). The results for the NEW MEXICO experiments were provided upon request by Gerard Schepers.

Christian Masson is a member of the editorial board of the journal.

This work is partially supported the Canadian Research Chair on the Nordic
Environment Aerodynamics of Wind Turbines and the Natural Sciences and
Engineering Research Council (NSERC) of Canada. Thanks for the great work
done by Matthew Churchfield and colleagues at National Wind Technology
Center, Boulder, CO, by establishing the open-source framework SOWFA. The
data used have been supplied by the consortium which carried out the EU FP5
project Mexico: “Model rotor EXperiments In COntrolled conditions”. Thanks a lot also to Gerard Schepers for providing results of the NEW MEXICO
experiment.

Edited by: Sandrine
Aubrun

Reviewed by: two anonymous referees

Carrión, M., Steijl, R., Woodgate, M., Barakos, G., Munduate, X., and Gomez-Iradi, S.: Computational fluid dynamics analysis of the wake behind the MEXICO rotor in axial flow conditions, Wind Energy, 18, 1023–1045, 2015. a, b

Gilling, L.: TuGen, Tech. rep., Department of Civil Engineering, Aalborg University, 2009. a

Glauert, H.: Airplane Propellers, in: Aerodynamic Theory, Vol. IV, Division L, edited by: Durand, W. F., 169–360, Springer, New York, 1935. a

ISO 61400: IEC 61400 – ONLINE COLLECTION – Wind turbines, International Standard, 23 April 2015, Edition 1.0, 1000 pp., 2015. a

Ivanell, S., Mikkelsen, R. F., Sørensen, J. N., and Henningson, D.: Stability analysis of the tip vortices of a wind turbine, Wind Energy, 13, 705–715, 2010. a, b, c

Jasak, H.: Error analysis and estimation for finite volume method with applications to fluid flow, Ph.D. thesis, Imperial College, 1996. a

Jha, P. K., Churchfield, M. J., Moriarty, P. J., and Schmitz, S.: Guidelines for Actuator Line Modeling of Wind Turbines on Large-Eddy Simulation-type Grids, J. Sol. Energ.-T. ASME, 136, 1–11, https://doi.org/10.1115/1.4026252, 2013. a

Krogstad, P.-Å. and Eriksen, P. E.: “Blind test” calculations of the performance and wake development for a model wind turbine, Renew. Energ., 50, 325–333, https://doi.org/10.1016/j.renene.2012.06.044, 2013. a

Mann, J.: Wind field simulation, Probabilist. Eng. Mech., 13, 269–282, https://doi.org/10.1016/S0266-8920(97)00036-2, 1998. a, b, c, d

Manwell, J., McGowan, J., and Rogers, A.: Wind Energy Explained: Theory, Design and Application, Wiley, 2010. a

Martínez-Tossas, L. A., Churchfield, M. J., and Leonardi, S.: Large eddy simulations of the flow past wind turbines: actuator line and disk modeling, Wind Energy, 18, 1047–1060, 2015. a

Martínez-Tossas, L. A., Churchfield, M. J., and Meneveau, C.: A Highly Resolved Large-Eddy Simulation of a Wind Turbine using an Actuator Line Model with Optimal Body Force Projection, J. Phys. Conf. Ser., 753, 082014, https://doi.org/10.1088/1742-6596/753/8/082014, 2016. a, b

Martínez-Tossas, L. A., Churchfield, M. J., and Meneveau, C.: Optimal smoothing length scale for actuator line models of wind turbine blades, Wind Energy, 20, 1083–1096, 2017. a, b, c, d

Meneveau, C., Lund, T. S., and Cabot, W. H.: A Lagrangian dynamic subgrid-scale model of turbulence, J. Fluid Mech., 319, 353, https://doi.org/10.1017/S0022112096007379, 1996. a

Muller, Y.-A., Masson, C., and Aubrun, S.: Large Eddy Simulation of the meandering of a wind turbine wake with stochastically generated boundary conditions, J. Phys. Conf. Ser., 524, 012149, https://doi.org/10.1088/1742-6596/524/1/012149, 2014. a, b, c, d, e

Nathan, J.: Application of the actuator surface concept in large-eddy simulations with regards to the near wake of wind turbines, Ph.D. thesis, ÉTS Montréal, 2018. a, b, c, d, e, f, g

Nathan, J., Forsting, A. R. M., Troldborg, N., and Masson, C.: Comparison of OpenFOAM and EllipSys3D actuator line methods with (NEW) MEXICO results, J. Phys. Conf. Ser., 854, 012033, https://doi.org/10.1088/1742-6596/854/1/012033, 2017. a, b, c, d, e, f

Nilsson, K., Shen, W. Z., Sørensen, J. N., Breton, S.-P., and Ivanell, S.: Validation of the actuator line method using near wake measurements of the MEXICO rotor, Wind Energy, 18, 499–514, 2015. a

Olivares Espinosa, H.: Turbulence Modelling in Wind Turbine Wakes, Ph.D. thesis, ÉTS Montréal, 2017. a, b, c, d, e, f, g, h, i

Pope, S.: Turbulent Flows, Cambridge University Press, 2000. a, b

Pope, S. B.: Ten questions concerning the large-eddy simulation of turbulent flows, New J. Phys., 6, 1–24, https://doi.org/10.1088/1367-2630/6/1/035, 2004. a

Rullaud, S., Blondel, F., and Cathelain, M.: Actuator-Line Model in a Lattice Boltzmann Framework for Wind Turbine Simulations, J. Phys. Conf. Ser., 1037, 022023, https://doi.org/10.1088/1742-6596/1037/2/022023, http://stacks.iop.org/1742-6596/1037/i=2/a=022023, 2018. a

Sanderse, B., van der Pijl, S., and Koren, B.: Review of computational fluid dynamics for wind turbine wake aerodynamics, Wind Energy, 14, 799–819, 2011. a

Sarmast, S., Segalini, A., Mikkelsen, R. F., and Ivanell, S.: Comparison of the near-wake between actuator-line simulations and a simplified vortex model of a horizontal-axis wind turbine, Wind Energy, 19, 471–481, 2016. a

Schepers, J., Boorsma, K., Cho, T., Gomez-Iradi, S., Schaffarczyk, P., Jeromin, A., Shen, W., Lutz, T., Meister, K., Stoevesandt, B., Schreck, S., Micallef, D., Pereira, R., Sant, T., Madsen, H., and Sørensen, N.: Final report of IEA Task 29, Mexnext (Phase 1): Analysis of Mexico wind tunnel measurements, Tech. rep., ECN, 2012. a

Shannon, C. E.: Communication in the presence of noise, Proceedings of the IRE, 37, 10–21, 1949. a

Shen, W. Z., Zhu, W. J., and Sørensen, J. N.: Actuator line/Navier–Stokes computations for the MEXICO rotor: comparison with detailed measurements, Wind Energy, 15, 811–825, 2012. a

Shives, M. and Crawford, C.: Mesh and load distribution requirements for actuator line CFD simulations, Wind Energy, 16, 1183–1196, 2013. a, b, c

Troldborg, N.: Actuator Line Modeling of Wind Turbine Wakes, Ph.D. thesis, DTU, Copenhagen, DK, 2009. a, b, c, d, e, f

Troldborg, N., Larsen, G. C., Madsen, H. A., Hansen, K. S., Sørensen, J. N., and Mikkelsen, R. F.: Numerical simulations of wake interaction between two wind turbines at various infl ow conditions, Wind Energy, 14, 859–876, 2011. a

Vanella, M., Piomelli, U., and Balaras, E.: Effect of grid discontinuities on large-eddy simulation statistics and flow fields, J. Turbul., 9, 37–41, https://doi.org/10.1080/14685240802446737, 2008. a

Vermeer, L. J., Sørensen, J. N., and Crespo, A.: Wind turbine wake aerodynamics, Prog. Aerosp. Sci., 39, 467–510, https://doi.org/10.1016/S0376-0421(03)00078-2, 2003. a

Warming, R. and Beam, R. M.: Upwind second-order difference schemes and applications in aerodynamic flows, AIAA J., 14, 1241–1249, 1976. a

Welch, P. D.: The Use of Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging Over Short, Modified Periodograms, IEEE T. Acoust. Speech, 15, 70–73, https://doi.org/10.1109/TAU.1967.1161901, 1967. a

Weller, H. G., Tabor, G., Jasak, H., and Fureby, C.: A tensorial approach to computational continuum mechanics using object-oriented techniques, Comput. Phys., 12, 620–631, https://doi.org/10.1063/1.168744, 1998. a

^{1}

OPENFOAM (Open source Field Operation And Manipulation) is a registered trade mark of OpenCFD Limited, producer and distributor of the OpenFOAM software via https://www.openfoam.com/ (last access: 1 March 2018).

^{2}

NWTC Design Codes (SOWFA (Simulator fOr Wind Farm Applications) by Matt Churchfield and Sang Lee) http://wind.nrel.gov/designcodes/simulators/SOWFA/ (last access: 1 March 2018). NWTC (National Wind Technology Center) is part of NREL (National Renewable Energy Laboratory) based in Golden, CO, USA.