3.2.1. Wind Turbine Wakes and Turbulence Modeling
Wake and turbulence modeling [
57,
58,
59,
281] are two of the most challenging research topics in the field of fluid mechanics and are considered unsolved problems in classic physics. However, a significant amount of contributions have been performed in the field during the last four decades resulting in two different types of WT wake modeling approaches. The first modeling approach is based on Engineering Wake Models (EWMs) [
18,
19,
21,
23,
282,
283,
284], which, in the present context, are analytical functions derived from empirical regressions, correlations, or simplified Newtonian-based laws of mass flow and momentum conservation, which describe in a simplified way the WT energy conversion and the WT wake evolution processes (e.g., the evolution of the streamwise wind speed at any radial position measured from the WT hub height) while making certain assumptions regarding the near and far wake expansion, the wake speed recovery, and the wake merging process, which in turn are a function of the WT aerodynamics and operation characteristics and the local atmospheric condition. The second modeling approach consist of advanced wake models that are based on the Reynolds-Averaged Navier-Stokes (RANS) equations or Large Eddy Simulations (LES), both complemented with turbulence models and subgrid-scale models respectively, and are combined with the concepts of actuator line and actuator disk approaches for the WT rotor modeling [
25,
26,
57,
285,
286,
287,
288].
The second approach (based on either the RANS or the LES approaches) has shown improved accuracy when describing the evolution and the side effects of WT wakes (in both steady and unsteady conditions) but requires a significant amount of computational resources, thus is commonly not adopted by wind farm designers nor by most of the authors of the reviewed works. The first approach (based on EWMs), although simple, physically limited and computationally inexpensive, has proven to describe the wind speed deficit in the far wake zone with reasonable accuracy [
60,
65,
66,
67,
68]. However, EWMs typically overpredict wake effects [
58,
60,
66,
281,
289]. Moreover, some commercial software packages (e.g., WindFarmer, WindFarm, WindPro, OpenWind) complement EWMs with Engineering Turbulence Models (ETMs) [
288,
290,
291,
292,
293,
294,
295], which describe the net TI evolution within the WT wake. Nonetheless, none of the reviewed works considered the inclusion of ETMs, since TI effects on the energy conversion process are often neglected.
The basic mathematical structure of the EWMs is shown in Equation (15):
where
is the horizontal wind speed in the wake at some axial (
[m]) and radial (
[m]) position downwind from the WT hub, and
is the velocity deficit, which represents the fraction of the free-stream wind speed that was affected during the energy conversion process. The wind speed in the wake exhibits a natural recovery due to the turbulent and viscous interactions of the wake with the unaffected neighboring free stream. The wake recovery is commonly modeled through the use of wake decay coefficients, which are assumed to be exclusively dependent on the ambient turbulence intensity levels [
19,
23,
282,
283].
The differences between common EWMs are essentially the way they approximate
. The first wake model, categorized as a Kinematic wake model [
59], was proposed by Lissaman in 1977 [
19,
296,
297], and later was updated by him in 1982 [
24]. He assumed that the wake of a single WT can be decomposed into two different main regions: the near and the far wake regions. For each region, self-similar wind speed profiles were assumed. Between these two regions is a transition zone, in which the flow profile changes from one self-similar profile to another. The wake is assumed to grow by turbulent entrainment at its edge, which introduces free stream momentum as well as mass into the wake, so these quantities are not conserved in the wake. Nonetheless, the momentum deficit in the wake, considering an idealistic constant pressure field, is conserved. The momentum deficit is defined by the integral shown in Equation (16) across the wake:
where the free-stream wind speed (
) was assumed constant and homogeneous. Equation (16), along with the wake growth (
), and the wake profile at a radius
r [m], measured from the wake centerline and given by a function of the form
(where
is the effective outer radius of the wake), are the most important mechanisms modeling the wake. The momentum deficit of the wake can be computed from the initial state of the wake, which can be directly related to the thrust force, or the thrust coefficient (
), of the WT rotor and its output power. The thrust coefficient of the WT represents the normalized axial force produced by the wind on the WT rotor, or conversely, of the rotor on the wind. The thrust coefficient can be viewed as a basic parameter representing the aero-structural loading of the WT and is dependent on the same variables as the power coefficient (e.g.,
,
, TI and
), as discussed in
Section 3.1.1. Similarly to the power coefficient, the thrust coefficient is strongly influenced by the adopted control strategy and is commonly reported by WT manufactures as an empirical function of the wind speed at hub height, while assuming standardized inflow conditions and disregarding any other dependence on the previously mentioned atmospheric variables.
Lissaman assumed that turbulence near the edge of the wake controls the lateral growth. This turbulence is a functional combination of: (1) the ambient turbulence intensity; (2) the turbulence generated by the shear due to the wind speed gradients within the wake itself (or wake shear); and (3) the mechanically added turbulence generated by the viscous interactions of the free-stream flow and the WT rotor (i.e., the mechanically generated turbulence). These terms have a different relative importance in the different regions of the wake. The ambient turbulence (assumed isotropic and uniform in the free stream) controls the wake growth after the mechanically added turbulence and the turbulence generated by the wake shear have decayed. Thus in the near wake region, the three turbulent components drive the wake growth, and in the far region only the ambient component drives the wake growth.
At the beginning of the near wake region, the radius of the wake is given by
and at the end of the near wake region is given by
, where
R [m] is the WT rotor radius and
is the inverse of the normalized wake speed at the wake core (
). In the near wake region, the transversal wake profile has a potential core of uniform velocity (
) and radius
. Outside the wake core (e.g.,
), the transversal wake profile is defined by a characteristic shear layer, which in turn defines the wake velocity deficit, as shown in Equation (17):
where
, so that
is the width of the turbulent mixed zone in the transversal wake profile. The shear layer expands outward (
i.e.,
) and inward (
i.e.,
) from the edge of the wake as it develops in the streamwise direction (
[m]). It is considered that the near wake region terminates when the outer shear layer meet the centerline of the wake, so the center velocity deficit starts to decrease with increasing distance from the WT. The downstream extent (
[m]) of the near wake can be computed using Equation (18):
where
,
, and
are the wake growth contributions due to ambient turbulence, due to the turbulence generated by the wake shear, and due to the mechanically added turbulence from the WT rotor. Moreover,
is the rotor solidity,
[%] is the ambient turbulence intensity,
is the thrust (or drag) coefficient of the blade profile, and
is the tip speed ratio. The wake growth rate is assumed linear and given by
.
In the transition region, the velocity deficit is given by a linear combination of profiles at the extremes of the transition region. At the beginning of the transition region, the radius of the wake is
and at the end of the transition region is
. The extent of the transition region can be calculated using Equation (19):
The wake growth rate at the transition region is the same as in the near wake region, given by .
In the far wake, the velocity deficit is given by Equation (20):
where
is the velocity deficit at the wake core. The wake growth rate in the far region is governed only by the ambient turbulence intensity, resulting in
. Thus the radius of the far wake at any point downstream (
) can be calculated as
.
The ground has the effect of reducing vertical downward wake growth, both by suppression of the turbulence near the ground and by the direct presence of the ground. The effects of the ground were modeled by classical reflection techniques (or imaging techniques) [
281], in which an imaginary WT is added below the ground (in a symmetrical fashion), so the wakes of both WTs expand beyond the ground and then their velocity deficits are added, thus the momentum deficit of the original wake is conserved. However, as pointed out by Crespo
et al. [
281], this approach is not completely valid as the wake velocity deficit is not actually conserved because of the friction (viscous) effects at the surface. To overcome this issue, Crespo
et al. [
281] considered an antisymmetric WT wake, so that velocity deficits are subtracted and result in a zero perturbation at the ground. However, although this alternative procedure eliminates the previously mentioned issue, it is not clear if it is valid throughout the wind field, where surface effects might not be important. Finally, the tower wake was introduced by treating the tower as a drag producing cylinder, which can be regarded as an approximately two-dimensional problem and has appropriate two-dimensional growth rates (
i.e., the wake of the tower expands only in the plane parallel to the ground).
Another interesting wake model, categorized as a parabolic field model [
59], was developed by Ainslie [
26] in 1988. The model considers that the wake is axisymmetric, fully turbulent, having nil circumferential velocities, and stationary with time. In addition, the pressure gradients in the co-flowing fluid outside the wake were assumed negligible. The model is based on the Reynolds-Averaged Navier-Stokes equations, which are complemented with an eddy viscosity model as turbulence closure. The eddy viscosity is predicted by a simple analytical form of the Prandtl’s free shear layer model and includes contributions by the ambient turbulence. The set of equations were simplified by adopting the thin shear layer approximation and by disregarding the viscous terms. The complete Ainslie far wake model is shown in Equations (21) and (22):
where
[m/s] is the axial mean velocity,
[m/s] is the radial mean velocity,
x [m/s] is the axial distance coordinate (downwind from the WT position),
r [m] is the radial distance coordinate (measured from the wake centerline),
[m
2/s
2] is the Reynolds stress cross-correlation,
[m
2/s] is the eddy viscosity,
and
are suitable length and velocity scales describing the wake shear layer, and
[m
2/s] is the ambient turbulence contribution to the eddy viscosity.
The length and velocity scales were taken to be proportional to the wake width [m] and the velocity difference [m/s] across the wake shear layer (where is the centerline wake velocity), respectively. is given by the atmospheric eddy diffusivity of momentum , where k is the Von Karman constant (), [m/s] is the friction velocity, z [m] is the height above the ground, L [m] is the Monin-Obukhov length, and is a similarity function reflecting the influence of the thermal stability condition on the turbulent mixing process (e.g., for neutral conditions). The value of the friction velocity depends on the atmospheric surface shear layer, which can be approximated with the use of the logarithmic laws. For neutral conditions, the logarithmic law states that the axial wind speed ( [m/s]) at any height within the atmospheric surface layer (e.g., up to 100 m) is governed by the roughness length ( [m]) and the friction velocity as given by .
A modification is required in order to use the eddy viscosity equation in the near wake region (which length was taken to be approximately five rotor diameters), since a lack of equilibrium between the mean velocity field and the turbulence field was observed. A filter function
, shown in Equation (23), was proposed to modulate the build-up of turbulence in the shear layer of the wake:
The resultant eddy viscosity closure, incorporating the previous approximation, results in Equation (24):
where
is a constant, which is a property of the shear layer and largely independent of the ambient turbulence and the WT characteristics. Ainslie conducted wind tunnel experiments for different inflow conditions at low ambient turbulence intensities and found that the value
k1 = 0.015 was adequate.
The solution of the Ainslie model should be calculated considering downstream distances of more than two rotor diameters, where pressure gradients no longer dominate the flow. The wake transversal profile was defined using a Gaussian profile of the form , where , [%] is the ambient turbulence intensity, is the thrust coefficient of the WT, and .
To complement the wake model, Ainslie provided a simple model for the treatment of the wake meandering [
59,
281], a phenomenon in which the wake meanders relative to an observer fixed on the ground due to fluctuations in wind direction. The wake meandering results in a reduced centerline velocity deficit measured by the fixed observer when compared with the centerline velocity deficit of a stationary case, because the measurement point sweeps across a region of the wake profile during the averaging period of the measurement. When meandering occurs (e.g., in neutral and unstable atmospheres), the deficit measured by a fixed observed can be determined using Equation (25):
where
is the uncorrected centerline velocity deficit and
is the standard deviation of wind direction fluctuations, measured over the same averaging time as the wake measurements (but not including very short time timescale fluctuations). Finally, the Ainslie model do not account for ground effects, tower effects, and variations in turbulent properties across the wake in its basic formulation.
Another relevant wake modeling approach was proposed by Larsen, who derived two similar wake models in 1988 [
25], which were updated by him later in 2009 [
287]. The models are based on the turbulent boundary equations (a simplified version of the Reynolds-Averaged Navier-Stokes equations), from which closed-from analytical solutions were obtained for the wake width as well as for the mean wind speed profile within the WT wake. Similarity assumptions and the Prandtl’s mixing length theory for the description of the turbulent stresses were considered. Both models have different levels of sophistication; the first model only takes into account the dominant terms of the turbulent boundary equations, whereas the second model takes into account the full system of equations. The free-stream airflow was considered turbulent, homogeneous, incompressible and stationary. Larsen simplified the original set of parabolic differential equations, expressed in cylindrical coordinates, by assuming an axisymmetric wake, by neglecting some terms (e.g., the pressure terms and distance terms of the form
), by conducting analyses of order of magnitude (e.g., the thin shear approximation), and by using similarity assumptions for the description of the wake profile evolution, resulting in the first order model. The analytical solution to the first order model, together with the semi-analytical calibration of the involved constants, is shown in Equations (26)–(28):
where
[m/s] is the axial wind speed magnitude at some axial (
x [m]) and radial (
r [m]) position within the wake,
[m/s] is the radial wind speed magnitude at some axial and radial position within the wake,
is the thrust coefficient of the WT,
A [m
2] is the rotor swept area,
,
D [m] is the rotor diameter,
,
,
is a measured axial velocity that depends on
[m], which represents a distance measured from the WT rotor, and
. The second order model of Larsen included the previously neglected distance terms of the form
and its analytical solution, together with additional assumptions, resulted in a very large model that is not described here for the sake of brevity.
The updated version of the Larsen model [
287] included two major improvements: an approach to merge wakes, which is necessary to quantify the performance of wind farms, and a proper calibration of the model with full-scale experiments. The updated first order model is shown in Equations (29) and (30):
where
,
,
,
is the wake radius at 9.6 rotor diameters downstream from the WT, which was calculated empirically form the Vindeby offshore wind farm data and expressed by
, where the
coefficients are resumed in [
287].
In order to calculate the effective wind speed of merged wakes in wind farms, Larsen considered the calculation of the mean effective wind speed reaching the rotor of a WT, given by either
or
, where the wind shear was assumed to be described by the logarithmic laws. The effective wind speed within the wind farm was calculated using a linear superposition of wakes, as shown in Equation (31):
where
is the number of upstream WT. A second approach was suggested by using the effective kinetic energy reaching the WTs, as shown in Equation (32):
Finally, like the Ainslie wake model, the Larsen models do not account for the ground effects or the tower effects in its basic formulations.
Another wake modeling approach was suggested by Frandsen
et al. in 2006 [
283]. For the description of single WT wakes, they considered expressions derived from the Betz-Lanchester theory, which link the thrust and power coefficient of the WT to the flow speed deficit in the wake. A cylindrical control volume with constant cross-sectional area equal to the wake area, with horizontal axis parallel to the mean wind velocity vector, and with no flow across its surface was used to perform a momentum balance between the free-stream flow and the WT wake, considering a simplified integral version of the momentum equation, to determine the momentum deficit in the wake. The simplified momentum equation disregarded pressure gradients (e.g., the control volume extends upwind,
, and downwind,
, far enough for the pressure to become equal to the free-stream pressure), the gravity force, viscous forces (e.g., shear forces in the cylinder surface were neglected), rotational effects, and tower effects. The free-stream wind speed was considered constant, homogeneous, stationary, and incompressible. The wake was considered axisymmetric and its transversal profile constant (
i.e., a top-hat profile was considered). The model defines the wake wind speed (
) as shown in Equation (33):
where
[m
2] is the swept area of the WT rotor,
is the area of the wake, which is assumed to be a function of the downstream distance (
x [m]), and “+” in Equation (33) applies when
and “−” when
. The wake area after it initially expands (assumed to happen at the rotor position, so
) is given by
, where
.
The evolution of the wake diameter downstream the WT was defined as shown in Equation (34):
where
is a wake decay factor,
,
,
, and the initial wake initial wake diameter is
. Nonetheless,
must be determined depending on the study case; for small
and large
, the decay factor
is of the order 10
.
The model was used to predict the performance of regular arrays of WTs. Three different scenarios were considered; in the first scenario, a row of WTs, with each WT exposed to multiple wakes, was studied and an analytical link between the wake expansion and the asymptotic flow speed deficit was derived. The wake-affected wind speed reaching a WT was derived from the momentum balance within the control volume, resulting in the expression given in Equation (35):
where
is the wake speed reaching a WT
,
is the wake area at the position of the WT
,
,
, and
is a constant separation distance among the WTs. The effects from boundaries, such as the ground, were included implicitly through the wake area growth
. For an infinitely large wind farm, it was assumed that there is an asymptotic, non-zero wake flow speed, thus
when
. In such cases, the wake speed and the wake area are given by Equations (36) and (37):
The right hand side of Equation (37) is a constant value, thus the wake cross-sectional area expands linearly with increasing downstream distance. The wake decay factor when is given by .
In the second scenario, the wakes from neighboring rows of WTs merge and the resultant wake can only expand vertically after that point. Since the area must, asymptotically, expand linearly (e.g., Equation (37)), the height of the wake must increase linearly. Therefore, the growth (
h) of the internal boundary layer in the limit for
is given by Equation (38):
where
is the dimensionless distance between neighboring rows,
and
are the axial distance and wake height when the neighboring wakes start to merge.
The third scenario is when the wind farm is infinitely large and the wake-affected flow is in balance with the atmospheric boundary layer. The interested reader is referred to [
283,
298] for a complete description of the model, which links the surface wind speed to the geostrophic wind speed by defining two flow layers divided at the WT hub height and by introducing the geostrophic drag law.
The most widely used EWM is the Jensen-Katic model [
21,
23,
282], as shown in
Figure 3. The Jensen-Katic wake model defines the velocity deficit as shown in Equation (39):
where
[%] is the ambient turbulence intensity, representing the effect of turbulent forces on the WT wake recovery and wake expansion,
is the thrust coefficient, and
D [m] is the WT rotor diameter.
is sometimes written as
, where
denotes the momentum entrainment constant (or wake decay coefficient) modulating both the wind speed recovery and the wake expansion, as first proposed in [
23,
298] and as shown in Equation (40):
where
is the height above ground [m], and
is the roughness length [m], representing the height above ground level at which a zero wind speed is found. Equation (40) implies that the wake recovery and expansion is lower with increasing height, since lower ambient turbulence intensities are found far from the surface. The wake expansion is commonly assumed to be proportional to the downwind distance (
x [m]) from the WT hub such that
, where
[m] is the radius of the wake and
[m] is the rotor radius.
Figure 3.
Timeline of Wake Models used in the WFDO Problem.
Figure 3.
Timeline of Wake Models used in the WFDO Problem.
Like Frandsen’s wake model, the Jensen-Katic wake model disregards pressure gradients, the gravity force, viscous forces (e.g., shear forces in the edge of the wake were neglected), rotational effects, ground effects, and tower effects. The free-stream wind speed was considered constant, homogeneous, stationary, and incompressible. Moreover, the wake was considered axisymmetric and its transversal profile constant (i.e., a top-hat profile was considered). The Jensen-Katic wake model relies on momentum conservation in a fully developed wake, thus its validity is limited to the far wake region. However, conflicting views about which is the axial length at which the wake can be considered to be in the far region arise in the literature. It is often considered that the far wake starts at an axial distance of about three to four rotor diameters downstream from the WT. A formal definition for such a distance, however, should be based on the necessary distance required for the reestablishment of the axial and radial pressure field, turbulent kinetic energy and angular momentum values, which in turn are a function of the WT operation mode, as previously described.
It has been observed that the total velocity deficit, composed of several merged wakes, depends mostly on the closest WT [
59,
299,
300]. However, the interaction of multiple wakes is not fully understood and is subject of many studies in the fluid mechanics field. Different methods for calculating the total effect of merged WT wakes have been proposed in the literature [
23,
269,
300]. The wake merging process is based on the concept of constructive superposition, where several wakes can be combined to produce a stronger wake. One effective and recurring way to calculate the effective value of the velocity deficit (
) in merged wakes is to sum the kinetic energy deficits of all WTs contributing to the wake (
i.e., the sum of the square of the velocity deficits) as shown in Equation (41):
where
is the velocity deficit contribution of the upstream WT
i, and
is the number of upstream WTs. The effect on the energy conversion of a WT operating under partial wake conditions is commonly modeled through the use of a geometric factor [
156,
178,
301] that affects the value of the effective wind speed deficit (
i.e., Equations (39) and (41)). The geometric factor is commonly defined as the ratio of the wake affected rotor area to the WT rotor area (
). Therefore, the modeled partial wake operation is less harmful, in terms of power conversion, than full wake operation. The assumption, however, does not account for the increased fatigue loading the WT will undergo while operating in partial wake conditions, which is important as denoted by the experimental findings of Seifert
et al. [
69]. Moreover, the tower and nacelle contributions to the wake are usually not considered explicitly by EWMs. Nevertheless, its effects can be included indirectly through the use of empirical values of
, which may contain contributions from all passive and active parts of the WT.
New semi-empirical approximations that complement the Jensen-Katic wake model have been proposed to describe: (1) the wake expansion and recovery as a function of the net TI (composed of both the ambient and the WT mechanically added turbulence intensities); (2) the effects of yaw misalignment and unsteady conditions in wake development; and (3) the effects of operating WTs under partial wake conditions [
282].