## Abstract

The actuator line technique was introduced as a numerical tool to be employed in combination with large eddy simulations to enable the study of wakes and wake interaction in wind farms. The technique is today largely used for studying basic features of wakes as well as for making performance predictions of wind farms. In this paper, we give a short introduction to the wake problem and the actuator line methodology and present a study in which the technique is employed to determine the near-wake properties of wind turbines. The presented results include a comparison of experimental results of the wake characteristics of the flow around a three-bladed model wind turbine, the development of a simple analytical formula for determining the near-wake length behind a wind turbine and a detailed investigation of wake structures based on proper orthogonal decomposition analysis of numerically generated snapshots of the wake.

## 1. Introduction and state of the art

Although offshore wind power only constitutes a small part of the total wind turbine capacity, it is expected to develop rapidly in the near future, and the potential offshore market is today the main driver for the development of large wind turbines. Offshore located wind turbines are always clustered in wind farms in order to limit the overall installation and maintenance expenses. Therefore, wake effects due to the mutual interaction between the wakes of the turbines are an inherent, albeit unwanted, effect, which the developer has to deal with. Wake effects cause both a decrease in total power production and an increase in self-generated turbulence, which, as a consequence, results in lifetime reductions [1–3]. The modelling of flow fields in wind farms is in many cases based on simple single-wake calculations combined with some simplified assumptions regarding the superposition of merging wakes [4–6]. As this type of method is typically based on simple analytical expressions, they constitute simple and fast analysis and design tools, which to a large extent have been implemented in commercial packages for the design of wind farms. However, they do not capture the actual turbulence characteristics in the wake, which may deviate significantly from that of the undisturbed flow field, and in most cases the models need to be calibrated against measurements in order to capture the local power production correctly. Furthermore, they disregard the important effect of wake meandering, which has been shown to have a major impact on wind turbine loading [7]. A more sophisticated type of simulation tool is based on the parabolized Reynolds-averaged Navier–Stokes equations [8,9]. This type of model is three-dimensional and capable of predicting effects such as the position of the wake centre line, the added turbulent kinetic energy (TKE) and the wake width. However, this type of model is limited by the use of the Reynolds-averaged eddy-viscosity approach, which is not capable of capturing correctly the spatio-temporal behaviour of the turbulence. There is therefore a need for developing more accurate simulation tools for describing the flow characteristics in wind farms. Such a tool, however, needs to cope with all the inherent scales in wind farms, where the wind turbines are subject to wind conditions of an atmospheric boundary layer. The size of this challenge is illustrated in table 1, which shows the various time and length scales associated with the problem. Here, it is seen that typical length scales range from millimetres, corresponding to the thickness of the boundary layer on the wind turbine rotor, to about 10 km, corresponding to the size of a wind farm. In principle, there are only two velocity scales, which are the tip speed, corresponding to approximately 100 m s^{−1}, and the wind speed, with a typical value of 10 m s^{−1}. Dividing length scales by velocity scales results in time scales that range from microseconds to more than 10 min. From table 1, it can be clearly seen that a full simulation involving a resolution of the detailed characteristics of the aerofoil boundary layer of each wind turbine for a complete wind farm results in a minimum of seven decades of length scales and eight decades of time scales. Hence, a full-blown direct numerical simulation solving the Navier–Stokes equations would involve a mesh consisting of more than 10^{21} mesh points and in the order of 10^{10} time steps, which is clearly beyond the capability of present-day computer technology. Instead, one has to come up with some simplified solution techniques that do not jeopardize the main characteristics and dynamics of the flow.

A way to do this is to employ large eddy simulation (LES), which reduces the computational costs by modelling the small-scale turbulence through a sub-grid scale (SGS) model, and to replace the actual geometry of the rotor blades by lines carrying body forces corresponding to the loading of the rotor blades. In this way, it is not necessary to make a detailed mesh for the rotor geometry, and the mesh points can be devoted to coping with the turbulence characteristics in the wake. This technique, referred to as the actuator line (AL) technique, was originally developed by Sørensen & Shen [10] and later implemented in a general aeroelastic simulation tool by Mikkelsen [11] and Troldborg [12]. The AL technique combines a three-dimensional Navier–Stokes solver with a technique in which body forces are distributed radially along each of the rotor blades. Thus, the kinematics of the wake is determined by a full three-dimensional Navier–Stokes simulation, whereas the influence of the rotating blades on the flow field is included using tabulated aerofoil data to represent the loading on each blade. The aerofoil data and subsequent loading are determined iteratively by computing local angles of attack from the movement of the blades and the local flow field. The concept enables the detail of the wake dynamics and the tip vortices and their influence on the induced velocities in the rotor plane to be studied. The AL and similar models have in the past decade made it possible to compute wake interaction, turbulent flow fields in wind farms and the interaction between wind farms and the atmospheric boundary layer [13–19].

The group of Porté-Agel [20–22] have successfully combined actuator disc and AL models with LES to study the effects of very large wind farms on the flow structures, including turbulent fluxes of momentum and heat in stably stratified atmospheric boundary layers. Other recent simulations of full wind farms include the work by Ivanell *et al.* [23], Nilsson *et al.* [24] and Yang *et al.* [25,26]. In spite of the important insights already obtained from the use of the AL model, there is a need for further simulations of wind turbine wakes for studying basic features of wakes in order to develop less complex engineering models for design purposes. This study shows some recent results from simulations combining instability theory and flow diagnosis tools with simulated wake characteristics.

The scope is here confined to studying wake fields for rotors subject to atmospheric turbulent inflow conditions, as this corresponds to a wind turbine located in the front row of a wind farm. The investigated wake will subsequently hit the following turbines. The basics of the AL technique is introduced and a number of findings from several detailed wake studies are presented in this article. The main features of the AL technique will be presented in §2, including its coupling to an aeroelastic tool, and then we will briefly describe how the technique is extended to incorporate atmospheric turbulence. A recent experimental study of the flow around a model wind turbine was carried out at the Norwegian University of Science and Technology (NTNU) wind tunnel at Trondheim, Norway [27]. This work provides some of the first detailed experimental results of near-wake characteristics of wind turbine wakes. The results are here employed as a validation of the AL model to reproduce near-wake characteristics. The validation is presented in §3. An important issue regarding wind farm aerodynamics is related to the properties of the near wake and the distance it takes before the wake has reached far-wake behaviour, i.e. a diffusion-dominated Gaussian velocity profile. This question has been elucidated by carrying out detailed AL simulations of a wind turbine rotor. The initial stable part of the wake is examined by employing the stability of the tip vortices. After the breakdown of the tip vortices, the wake evolves into small-scale turbulence. This process occurs in the transitional or near-wake region and has completed once the far-wake state is reached. This analysis, plus the derivation of a simple engineering formula for determining the distance, will be presented in §4. This work is a continuation of the investigation carried out by Sarmast *et al.* [28]. In §5, we go more into the details of the wake dynamics using the proper orthogonal decomposition (POD) technique to investigate large coherent structures and as an alternative means to determine the near-wake length, which corroborates the simple engineering formula. Finally, conclusions are drawn in §6.

## 2. Numerical modelling

### (a) Basic equations and numerics

The numerical simulations are performed using the three-dimensional flow solver EllipSys3D, developed in collaboration between the Technical University of Denmark (DTU) [29] and the former National Laboratory for Sustainable Energy (Risø) in Denmark [30]. This code solves the discretized incompressible Navier–Stokes equations in general curvilinear coordinates using a block structured finite-volume approach. EllipSys3D is formulated in primitive variables (pressure–velocity) in a collocated grid arrangement. The pressure correction equation is solved with the PISO algorithm and pressure decoupling is avoided using the Rhie–Chow interpolation technique. The convective terms are discretized using a hybrid scheme combining the third-order QUICK scheme and the fourth-order CDS scheme. This technique was employed as a compromise to limit unphysical numerical wiggles related to a pure fourth-order scheme as well as numerical diffusion due to the upwind biasing nature of the QUICK scheme.

LES applies a low-pass filter on the Navier–Stokes equations, which results in a filtered velocity field. The implemented LES is based on the mixed scale model by Ta Phuoc *et al.* [31]. The scales resolved by the grid are simulated directly by the filtered Navier–Stokes equations, whereas scales below the grid scale are modelled through an SGS model, which delivers the turbulence closure. The employed SGS model is based on an eddy-viscosity approach, where the SGS eddy viscosity is determined by
2.1
where *ρ* is the density of air, is the filtered velocity and is the filter cut-off length defined as (ΔVol)^{1/3}, where ΔVol is the volume of a given cell. *β* is a parameter taking values between zero and 1. In this study, *β*=0.5 is employed. is the kinetic energy and evaluated as
2.2
where yields the high-frequency part of the resolved velocity field found by subtracting the resolved velocity using a spatial test filter . The performance of the employed SGS model has recently been examined by Chivaee [32].

The flow field is thus approximated by solving the filtered three-dimensional incompressible Navier–Stokes equations for
2.3
and
2.4
where *p* is pressure and *ν* is eddy viscosity. A number of body forces (**f**_{WT} and **f**_{turb}) are explicitly introduced in the simulations to model the effect of the wind turbine and atmospheric turbulence. The individual body forces are described in detail in the following.

The flow solver is implemented in a normalized form, i.e. the directions () are scaled by the rotor radius (*R*) and the velocity components () are scaled by the freestream velocity *U*_{0} in the computations. Hence, *U*, *V* and *W* will be used interchangeably with *U*/*U*_{0}, *V*/*U*_{0} and *W*/*U*_{0} for the streamwise, lateral and vertical velocity components, respectively.

### (b) Wind turbine modelling

The influence of the wind turbine is modelled using the AL method. The AL imposes body forces distributed along rotating lines representing the aerodynamic loads of the wind turbine blades (see [10,33] for details). The body forces imposed in the Navier–Stokes equations are calculated through a coupling with the Flex5 code. Flex5 is an aeroelastic modal-based tool with 28 d.f. used for calculating deflections and load responses of the entire wind turbine structure [34]. Flex5 and AL are coupled through the aerodynamic loads in Flex5. These are usually calculated using blade element momentum theory, but have here been replaced with the Navier–Stokes-based flow field in a fully coupled manner in which the ALs deflect and rotate. The fully integrated model has mainly been used to analyse flows about fixed offshore wind turbines, but Flex5 also includes a module for floating tension-leg platform offshore wind turbines [35]. The three-dimensional hydrodynamic loads can be formulated using a frequency- and direction-dependent spectrum with wave loads computed from the wave kinematics, so the dynamics of the full coupled system is included in the wake aerodynamics. Coupling CFD and Flex5 has recently been done by Schløer *et al.* [36]. However, this study only concerns wake effects from fixed offshore located turbines. As such, Flex5 and EllipSys3D are two independent codes and the coupling is conducted by transfer of line coordinates and forces from Flex5 to EllipSys3D with velocities being transferred back to Flex5 for each time step. In fluid structure interaction methods, loose versus strong numerical coupling is a key issue, but for wind turbine aeroelasticity [37] loose coupling may provide good accuracy due to the limited magnitude of the reduced frequencies for wind turbine structures. The coupling applied in the present form is characterized as a loose numerical coupling. The fully coupled set-up gives additional dynamics to the wake simulations and their interaction, as the results are more realistic, yet less symmetric. Aeroelastic results from such fully coupled simulations are not described in this work, but reference is made to Andersen [38] for an analysis of the aeroelastic outputs from such coupled simulations.

The line velocities transferred serve to compute the local velocity triangle and relative blade velocity. The local velocity triangle relative to the rotating aerofoil segment is determined from the velocity triangle, where *V* _{x}(=*U*) and *V* _{θ} are the velocities in the axial and tangential directions, respectively, and the relative velocity at the blade section is given as
2.5
The flow angle is given as , where *Ω* denotes the rotational frequency. This gives a local angle of attack, *α*=*ϕ*−*γ*, with *γ* denoting the local pitch angle. Lift and drag forces are found from
2.6
where *C*_{L}(*α*,*Re*) and *C*_{D}(*α*,*Re*) are the lift and drag coefficients, respectively, as functions of the angle of attack (*α*) and Reynolds number (*Re*) given as tabulated data. The Reynolds number is based on chord length and relative velocity. **e**_{L} and **e**_{D} are unit vectors parallel to the lift and drag forces, respectively. Projection gives the axial and tangential forces, *F*_{x} and *F*_{θ}. The body forces are numerically smeared across a few cells to avoid singularities. The smearing is Gaussian distributed by applying a convolution to the local load **f**_{2D} using a regularization kernel (*η*_{ϵ}) as follows:
2.7
Here, *N*_{b} is the number of blades and *r*=∥**x**−*r***e**_{i}∥, where ∥∥ is the Euclidean norm, i.e. the distance between the grid point and the force points on the *i*th AL denoted by the unit vector **e**_{i}.

The advantage of representing the individual blades by line-distributed loads is that much fewer grid points are needed to model the influence of the blades, when compared with resolving and simulating the actual geometry of the blades. The AL model allows for detailed studies of the dynamics of the different wake structures, such as the tip and root vortices, using a reasonably low number of grid points. Furthermore, the model benefits from being applicable with simple structured grids and therefore issues connected to grid generation do not occur. The drawback of the method is that it relies on the quality of tabulated aerofoil data.

### (c) Applying atmospheric turbulence

Atmospheric turbulence is introduced by imposing small body force fluctuations (**f**_{turb}) into the flow (e.g. Gilling *et al.* [39] for implementation). The body forces are usually applied in a plane upstream of the wind turbine. The body forces are derived from turbulent fluctuations generated by the so-called Mann model (see [40,41] for details). The Mann model is based on a linearization of the incompressible Navier–Stokes equations and generates a three-dimensional field of all three velocity components. The generated turbulence assumes Taylor's frozen turbulence hypothesis, which links time and space. Second-order statistics (variance, cross-spectra, etc.) are matched to those occurring in a neutral atmosphere and the generated turbulence is homogeneous, anisotropic and stationary.

## 3. Comparison with experiments

The EllipSys3D code has previously been validated against experimental data in [11,12]. Here, a recent, not yet published, comparison with experiments is presented. The performance and the wake development behind a scale model wind turbine are investigated numerically. The investigated set-up refers to a series of experimental measurements of a model scale turbine conducted by NTNU in a low-speed wind tunnel [27]. The modelled wind turbine is a three-bladed horizontal axis wind turbine. The wind turbine blades use 14% thick NREL S826 aerofoils along the span and operate at an effective sectional Reynolds number around *Re*=10^{5}. The details on performance of the aerofoil at low Reynolds numbers can be obtained from [42]. The experiment was performed as a blind-test comparison and gave a unique comparison with some very detailed and accurate experimental wake data.

The computational domain is constructed to replicate the size of the NTNU wind tunnel. It is a regular Cartesian grid composed of a total of 24.5 million mesh points with each blade represented by 43 AL points. The grid is distributed equidistantly around the wind turbine and its resolution is chosen to be similar to the distance between the consecutive points along the AL, e.g. Δ*x*=*R*/43. The effect of the tower is considered as a fixed AL in the simulations. Considering the boundary conditions, the velocity at the inlet is assumed to be uniform in the axial direction. A convective outflow boundary condition is used at the outlet. The wind tunnel walls are included using slip Euler boundary conditions. The wind tunnel turbulence (*T*_{i}≈0.23%) is modelled by introducing synthetic Mann turbulent fluctuations at approximately 2*R* upstream of the wind turbine.

The performance of the wind turbine is presented in figure 1*a*,*b*. The power coefficient () increases with the tip speed ratio (λ=*RΩ*/*U*_{0}) before it reaches a maximum value at the design tip speed ratio (λ=6) and then decreases for increasing tip speed ratio. The power predictions are generally in good agreement with the experiments, although small deviations can be observed for tip speed ratios of 4−5 where the computations under-predict the power output production by 7%. However, the estimated thrust coefficient () is consistently low by a shift of Δ*C*_{T}=0.06 in the normal operating tip speed range, which means that the thrust coefficients are underestimated by approximately 9%.

Considering the wake development behind the rotor, it is expected that both experiments and predictions produce fairly symmetric profiles. Figure 1*c*,*d* shows the velocity deficit and streamwise turbulent stress for λ=6 at three different positions behind the rotor, i.e. *x*/*R*=2,6 and 10. The rotor operates at its optimum, because the geometry was laid out to produce a constant pressure drop across the rotor. Hence, it is expected that the wake deficit profile close to the rotor, i.e. at *x*/*R*=2, is almost top-hat shaped. Comparing the profiles, it is observed that small differences exist which are related to the differences in the measured and computed thrust forces applied to the flow (figure 1*c*). Further downstream, a bell-shaped wake deficit is observed at *x*/*R*=10, indicating that the flow is close to a fully turbulent state, where the wake is dominated by small-scale turbulence structures. The streamwise turbulent stress is small close to the wind turbine, except in the thin region where the tip and root vortices are present. It is also observed how the distinct peaks in the streamwise turbulent stress are smoothed further downstream by the spanwise turbulent diffusion.

## 4. Near-wake analysis

It is very important to determine the optimum distance between neighbouring turbines, when clustering wind turbines in wind farms. A too large distance will result in additional installation costs, whereas a too short distance will result in increased fatigue problems for turbines located in the wake of other turbines. The latter is in particular a problem for the second turbine in a row as it potentially risks being hit by the tip vortices from the turbines located at the front of the wind farm. It is therefore important to determine the lifetime of the tip vortices and the parameters that govern their breakdown into small-scale turbulence. Figure 2 shows the wake development after a wind turbine in which the organized tip vortices break down and transform into small-scale turbulence structures.

From stability analysis [43], it is known that helical tip vortices are inherently unstable and that they will break down at some distance behind the wind turbine. However, the analysis does not predict the location of the breakdown point or the distance from the breakdown point to the location where the wake reaches its fully turbulent state. By combining numerical stability analysis with results from basic rotor aerodynamics, however, it is possible to arrive at a simple analytical expression for determining the length of the near wake behind a wind turbine. The analysis in this section will be carried out in two steps. First, the distance from the rotor to the location where the tip vortices break down will be determined as a function of the operating conditions of the rotor. Next, the distance from the breakdown point to the location where the flow reaches the far wake will be determined.

### (a) Stability analysis of tip vortices

It is possible to study the stability properties of the tip vortices by introducing small perturbations. The initial derivations of the presented analysis follow the work of Ivanell *et al.* [23] and Sarmast *et al.* [28]. In these studies, the simulations are performed using the Tjæreborg wind turbine operating at the optimum power condition (*C*p=0.49) for a wind speed of *U*_{0}=10 m *s*^{−1} and tip speed ratio of λ=7. The computations are performed in a polar grid with each blade represented by 97 grid points along the AL. The computational grid is equidistant in the proximity of the rotor, and the majority of the grid points are distributed in the near wake in order to capture the gradients and resolve tip and root vortices. The aim of the stability analysis is to study which stability modes are present and to what extent they grow in order to quantify frequencies leading to a breakdown of the vortex spirals. The perturbations are introduced in the axial direction only. The disturbance is positioned close behind the tip and results in a spatially developing disturbance wave on the spiral vortices. Using a number of velocity fields taken equidistantly over one period, the response from a specific perturbation frequency is evaluated using Fourier analysis. First, a steady-state solution is found, after which a time-resolved computation is performed to reach a periodic solution. *N* flow fields are extracted equidistantly spaced in time for one period and used to calculate the Fourier coefficients. The amplitude of the perturbation is extracted in order to evaluate the growth along the vortex spiral. This is determined by identifying the maximum amplitude of the perturbation frequency at each *x*-position, i.e. at each position in the flow direction. Figure 3 is a semi-logarithmic plot showing the development of the amplitude of the first harmonic due to disturbances of dimensionless frequencies *f*_{c}=2 and *f*_{c}=5 using different perturbation amplitudes. In figure 3*a*,*b*, the abscissa is the axial direction made dimensionless with the radius of the rotor. The rotor is located at *x*/*R*=0. From the figure, it can be clearly seen that the amplitudes develop with an exponential growth rate. More information regarding the computations can be obtained from Ivanell *et al.* [23] and Sarmast *et al.* [28].

The curves are normalized by the amplitude at *x*/*R*=0.5 to compare the growth from computations with different perturbations. The result shows that the growth rate is independent of the amplitude of the perturbation, as the slopes of the data from all perturbation amplitudes are equal up to the point where nonlinear effects become important. Denoting the initial perturbation at *x*=*x*_{0} by *A*_{0} and the amplitude at position *x* as *A*, the amplification reads
4.1
where *U*_{c} is the convective propagation velocity of the vortices and is the dimensional growth rate of a mutual induction instability. The origin of the instability essentially results from the mutual induction of the vortices in a form that it is analogous to the leapfrogging motion of two inviscid vortex rings. This phenomenon is known to occur in a row of equidistant and identical vortices, where amplifications of small perturbations cause the vortices to oscillate in such a way that neighbouring vortices approach each other and start to group in pairs [44]. By analysing the stability of rows of identical vortices, Lamb [45] found the maximum non-dimensional temporal growth rate to be σ=*π*/2. Similar conclusions were reached by Sarmast *et al.* [28] and Leweke *et al.* [46] for the temporal growth rate of rows of helical vortices.

Following the analysis of Lamb [45] (§156, pp. 225–226), making the growth rate dimensionless by the distance between the vortices, *h*, the strength of the vortices, *Γ*, and the convective velocity yield
4.2
Here, the amplification rate is related to the strength and mutual distance of the vortices. However, it is more practical to express it in terms of the operational conditions of the rotor. To accomplish this, some basic results from rotor theory will be exploited. Following Okulov & Sørensen [43], by assuming the wake consists of a system of tip vortices of strength *Γ* and a single root vortex of strength −*N*_{b}*Γ*, where *N*_{b} is the number of blades, the following relationship can be obtained:
4.3
Assuming that the rotor is loaded with a constant circulation, the following expression for the thrust is obtained [47]:
4.4
where *δ* denotes the inner non-loaded part of the rotor, which is here modelled by a solid body rotation of the core of the root vortex. Introducing further the thrust coefficient *C*_{T}, and dimensionless circulation (*q*=*N*_{b}*Γ*/2*πRU*_{0}), gives
4.5
Assuming a root vortex core of 10–20% of the rotor radius and a typical *q*-value of the order of 0.15, it is readily seen that the last term in the bracket is an order of magnitude smaller than the tip speed ratio. Therefore, this term will be neglected in the following. Combining equations (4.1)–(4.5), the following expression for determining the axial position of a given amplitude amplification can be derived:
4.6
where .

### (b) Breakdown position of tip vortices

The breakdown position of tip vortices is assumed to take place where linear amplification reaches its maximum. From figure 3, it is seen to be reached when the amplitude amplification equals the ratio between the original perturbation and the undisturbed wind speed, i.e. when
4.7
Assuming further that the turbulence intensity *T*_{i} is proportional to the initial perturbation, the following relationship holds:
4.8
where *C*_{1} is a proportionality constant. Combining equations (4.6)–(4.8), the following expression for the position where the tip vortices break down can be obtained:
4.9
It should be noted that the minus sign is due to the fact that the turbulence intensity is always less than 1 and the natural logarithm as a consequence is negative. This expression gives a measure of the position where the helical tip vortices break down as a function of the intensity of the ambient turbulence level and of parameters depending uniquely on the turbine's operational characteristics. In order to evaluate the unknown parameters and *C*_{1}, numerical simulations were conducted on the small-scale wind turbine using different intensities of the inflow turbulence. The modelled turbine specifications are those mentioned in §3. The rotor operates at optimum performance for *U*_{0}=10 m s^{−1} corresponding to a tip speed ratio of λ=6. Four simulations with different inflow turbulence intensities were performed and an overview of the results are presented in table 2. The constants *C*_{1} and can be determined from equation (4.9) to be *C*_{1}≈0.30 and . An additional approach to obtain the convective velocity is using the data from figure 3 together with equation (4.6). As previously mentioned, the Tjæreborg wind turbine is used to study the stability properties of the wind turbine wake. The wind turbine operates at its optimum condition (*U*_{0}=10 m *s*^{−1}, λ=7.07). The convective velocity is then computed to be . The variation in the values is related to the different operating condition of the wind turbines.

### (c) Length of near wake

This section is dedicated to the analysis leading to the determination of the transition part of the wake, which is here defined as the distance from the breakdown of the vortices (with approximately constant velocity deficit) to where a fully developed Gaussian velocity profile is attained. Thus, the near-wake length consists of the sum of the distance to the breakdown point (equation (4.9)) and the distance from this point to where a Gaussian wake deficit is first observed. Assuming that the latter is only affected by the incoming turbulence intensity, an extra term can be added to equation (4.9), which is assumed to depend directly on the natural logarithm of the turbulence intensity. The expression governing the length from the breakdown point to the start of the far wake is derived as follows:
4.10
where *C*_{2} is a proportionality constant. In the far wake, the azimuthally averaged axial velocity deficit, Δ*u*, is expected to follow a self-similar Gaussian shape [48]
4.11
where denotes the maximum deficit at a given axial position. *η*=*r*/*r*_{e} is the dimensionless radial distance normalized by *r*_{e}, the distance at which , where e denotes the Euler number. To determine how close the profile is to a Gaussian profile, a correlation coefficient, , is introduced. The correlation coefficient measures the degree of linear dependency between the computed velocity deficit and the analytical expression (equation (4.11)). It is computed for different turbulence intensities and is shown in figure 4*a* as a function of the dimensionless axial distance from the breakdown position and normalized by the natural log of the turbulence intensity. Clearly, the three simulations with a significant amount of turbulence collapse to a common line, whereas *T*_{i}=0.2% appears too small. Defining the onset of the far wake to be located where a correlation coefficient of 99% is achieved, an approximate value *C*_{2}≈−5.5 is obtained, when excluding the very low turbulence simulation. Table 2 shows the different *C*_{2}-values for a 99% correlation. The table also shows the simulated and predicted near-wake lengths. Figure 4*b* shows the mean velocity profiles for different positions inside the transitional wake region. The grey curves indicate the profiles upstream (close to the breakdown position), while the black curves show far-downstream positions. As more and more downstream profiles are considered, the curves start to collapse into a single curve in which the wake reaches its equilibrium self-similar state. From this figure, it is now possible to conclude that the 99% correlation curve gives a streamwise position in which the wake is at the border between the near wake and far wake. By adding equations (4.9) and (4.10), and introducing the computed values of *C*_{1} and *C*_{2}, the following rough-and-ready equation can be established for determining the total length of the near wake:
4.12
This expression gives a measure of the length of the near wake as a function of the intensity of the ambient turbulence level, *T*_{i}, and of parameters depending uniquely on the turbine's operational characteristics. As the constants are taken partly from the above analysis and partly from computations of the model rotor, the generality of the expression should be further validated using data from experiments and other simulations.

## 5. Dynamic wake characterization

This section employs POD to investigate the large coherent turbulent structures in the wake and how their spatial development relates to the determined near-wake length given in the previous section.

In order to simulate the behaviour of a modern wind turbine, an upscaled version of the NM80 turbine is used in the following simulations. NM80 is a three-bladed horizontal axis wind turbine with a blade length of *R*=40 m and rated to 2.75 MW at a hub velocity of *U*_{hub}=14 m *s*^{−1}. The turbine is upscaled from the original NM80 rated at 2.00 MW. The two-dimensional aerofoil data have been corrected to account for three-dimensional effects. The NM80 turbine is proprietary to Vestas Wind Systems A/S, but reference is made to the *DAN-AERO MW Experiments* [49].

The wake behind a single wind turbine has been simulated in a *x*×*y*×*z*=30*R*×30*R*×70*R* Cartesian domain. The mesh consists of 71 million grid points, yielding a rotor resolution of 25 cells per blade. Symmetric slip boundary conditions are used on the vertical and lateral sides, while the inflow is uniform and a convective boundary condition is applied on the outflow. Two different simulations are performed with turbulence intensity of *T*_{i}=7.1% and *T*_{i}=9.7% based on velocity fluctuations derived from a Mann turbulence box and introduced continously 4*R* upstream from the rotor.

POD is used to analyse the resulting flow fields. POD is a statistical method yielding an optimal linear subspace, optimal in terms of the variance, which here corresponds to the TKE content. POD can be considered as an energy filter, which reveals the large coherent turbulent structures [50]. Hence, POD ensures that the sorted modes in the optimal linear subspace contain the most energetic states. The theoretical background is omitted for brevity, but reference is made to the seminal work by Sirovich [51], Aubry *et al.* [52] and the general overview given by Berkooz *et al.* [53].

POD is employed on all three velocity components extracted over *t*=1500 s in planes behind the wind turbine, which yields a total of about 3000 POD modes. Figure 5 shows the evolution of the energy content and the spatial POD modes for the first (figure 5*b*), third (figure 5*c*) and eighth (figure 5*d*) POD modes for *T*_{i}=7.1%. The selected modes exhibit distinct spatial structures, i.e. monopole, dipole, quadrupole and hexapole structures. The first, and most energetic, POD mode clearly displays two distinct regions around the rotor perimeter, which stems from the tip vorticity. The absence of clearly defined tip vortices is mainly due to the lower resolution of this simulation. Additionally, distinct root vorticity is present. Similar structures are present 6*R* behind the rotor, although it quickly smears into two larger and larger regions as the turbulent mixing increases further downstream due to the wake breakdown initiated by the added atmospheric turbulence. Similar trends are observed for the second case with *T*_{i}=9.7% (not shown), although the distinct root structures are only present at *x*/*R*=2. It is important to note that the red and blue regions fluctuate between positive and negative according to the temporal eigenfunctions. Furthermore, the POD modes are orthogonal, so the second POD mode is comparable to the first with a phase shift of *π*/2. Therefore, only three selected modes are shown and not the orthogonal counterparts. The combination of the first two POD modes governs the very large scale motions, which translates the entire wake and is often collectively referred to as wake meandering [7]. However, it is also important to note that not all POD modes can necessarily be attributed to actual large coherent structures, but for higher POD modes the spatial distribution becomes a mere result of lumping TKE together in an optimal sense. The third POD mode shows four distinct regions along the rotor radius. This quadrupole structure is maintained up to *x*/*R*=26, where the third most energetic POD mode becomes a large central region, which pulsates back and forth. Again, a similar change is seen for *T*_{i}=9.7%, where the large central monopole appears as early as *x*/*R*=14. This change in spatial distribution of the large coherent structures in the third most energetic POD mode is due to redistribution of the TKE between different scales. A similar redistribution of energy has been reported by Johansson *et al*. [54]. The eighth POD mode contains six well-defined regions around the tip region, which smears further into the wake, but maintains six distinct regions up to *x*/*R*=44, before eventually resulting in two large regions. Once again, the well-defined regions persist for a shorter distance for *T*_{i}=9.7% as the quadrupole is only clearly visible for *x*/*R*=2, after which a couple more unstructured distributions appear followed by two large regions. Hence, the POD analysis yields specific insights into the spatial development of the large turbulent structures as the wake breaks from a stable wake with distinct tip and root vorticity over a transition region before the far wake.

This spatial development is also revealed by eigenvalues or TKE of the different POD modes. The energy contribution of the individual modes is also given as a function of distance behind the rotor and it highlights a significant change in the energy contribution. The first POD mode contains 16% of the total TKE 2*R* behind the rotor, where the tip and root vorticity dominate. Similarly, the third and eighth POD modes capture 7.0% and 2.6%, respectively, of the TKE. The energy content reaches a minimum at 16*R*, 14*R* and 12*R* behind the rotor for the first, third and eight POD modes, respectively. The energy content of the individual POD modes at the minimum constitutes only 9.7%, 3.4% and 1.8%, respectively. The minimum for the higher turbulence intensity of *T*_{i}=9.7% consistently occurs 2*R* earlier at 10*R*–14*R* for the three selected POD modes. The distance of 12*R*–16*R* is in very good agreement with the end of the near-wake region given by equation (4.12), which yields *l*_{near wake}=16.1*R* for *T*_{i}=7.1%, λ=7.2 and *C*_{T}=0.7. The near-wake length of *l*_{near wake}=14.3*R* for *T*_{i}=9.7% with the same tip speed ratio and thrust coefficient is also in very good agreement. Furthermore, POD reveals how the breakdown initiates at the smaller scales and spreads to the largest scales as the entire wake breaks down, as expected. So POD is applicable for determining the same near-wake length by examining the energy content of the POD modes. More POD modes are required to contain the TKE as the tip and root vorticities break down and to ensure turbulent mixing to obtain a Gaussian velocity profile on the transition to the far wake. The redistribution is gradually reversed in the far wake as the energy content of the highest POD modes increases again.

Hence, the spatial development of the POD modes visualizes the breakdown process as the wake transitions from the near wake into the far wake, as the large coherent turbulent structures change in shape and energy content. The derived dynamic and turbulent structures revealed by POD corroborate the analysis from the previous sections, which was based on mean values.

## 6. Conclusion

This paper demonstrates the utility of the AL technique to simulate detailed flow characteristics of wakes behind wind turbines. A short overview of state of the art of the AL technique is given and illustrated with some recent results. The advantage of the technique is that it allows detailed unsteady simulations to be carried out without having to resolve the boundary layer on the rotor blades in detail and that it can be coupled to existing aeroelastic simulation tools. Numerically obtained results are compared with measurements of a three-bladed model rotor. The computed power and thrust coefficients are generally in good agreement with measurements, with maximum deviations up to 9%. The comparison of the wake characteristics shows a very good agreement between computed and measured values of axial velocity profiles and Reynolds stresses in the near wake.

Based on the computed results, a simple expression for determining the length of the near wake has been derived. The derivation relies partly on numerical stability analysis of a perturbed flow field and partly on the computational results of the model rotor. The computations show that the first part of the wake length is inversely proportional to the number of rotor blades, tip speed ratio and thrust coefficient, and that the total length can be approximately expressed in terms of turbulence intensity.

An upscaled version of the NM80 wind turbine was used to investigate the features of the flow behind a modern wind turbine. POD was applied on numerically generated velocity components for single-wake simulations to study the dynamic characteristics of wakes. Applying POD on the computed results from the single-wake simulation confirms the main properties of the simple engineering formula for determining the length of the near wake. Thus, the energy content of the top POD modes show a clear minimum around the distance determined by the derived expression for the near-wake length, demonstrating that more POD modes are needed in order to capture all the energy as the wake breaks down into small-scale turbulence.

## Funding statement

This work has been carried out with the support of the Danish Council for Strategic Research for the project Center for Computational Wind Turbine Aerodynamics and Atmospheric Turbulence (grant no. 2104-09-067216/DSF) (COMWIND: http://www.comwind.org) and the Nordic Consortium on Optimization and Control of Wind Farms funded by the Swedish Energy Agency. Computer time was granted by the Swedish National Infrastructure for Computing (SNIC). Furthermore, proprietary data for the Vestas NM80 turbine have been used.

## Footnotes

One contribution of 17 to a theme issue ‘New perspectives in offshore wind energy’.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.