Journal topic
Wind Energ. Sci., 4, 633–644, 2019
https://doi.org/10.5194/wes-4-633-2019
Wind Energ. Sci., 4, 633–644, 2019
https://doi.org/10.5194/wes-4-633-2019

Research article 18 Nov 2019

Research article | 18 Nov 2019

# On the self-similarity of wind turbine wakes in a complex terrain using large eddy simulation

On the self-similarity of wind turbine wakes in a complex terrain using large eddy simulation
Arslan Salim Dar1,2, Jacob Berg1, Niels Troldborg1, and Edward G. Patton3 Arslan Salim Dar et al.
• 1Technical University of Denmark, DTU Wind Energy, Risø campus, Frederiksborgvej 399, 4000 Roskilde, Denmark
• 2Wind Engineering and Renewable Energy Laboratory (WIRE), École Polytechnique Fédérale de Lausanne (EPFL), 1015 Lausanne, Switzerland
• 3National Center for Atmospheric Research, Boulder, CO, USA

Correspondence: Jacob Berg (jbej@dtu.dk)

Abstract

We perform large eddy simulation of flow in a complex terrain under neutral atmospheric stratification. We study the self-similar behavior of a turbine wake as a function of varying terrain complexity and perform comparisons with a flat terrain. By plotting normalized velocity deficit profiles in different complex terrain cases, we verify that self-similarity is preserved as we move downstream from the turbine. We find that this preservation is valid for a shorter distance downstream compared to what is observed in a flat terrain. A larger spread of the profiles toward the tails due to varying levels of shear is also observed.

1 Introduction

Rotor wakes have a consequential impact on the efficiency of a wind farm, as the turbines standing in wake generally face lower wind speeds along with enhanced turbulence levels . The particular dynamics of wakes are strongly influenced by the terrain characteristics, such as surface vegetation and sloping terrain. Although a lot of emphasis has been on understanding wakes in flat terrain over the past decade , complex terrains now finally get the attention they deserve. This is partly due to a prospective shift in development of wind farms from flat to complex terrains caused by saturation of ideal flat terrains and increasing development of wind energy over the past 2 decades ; it is also partly due to the recent observational and numerical developments. Understanding wakes from turbines in complex terrains, therefore, becomes important for understanding the interaction between terrain and wakes, as well as for better resource assessment and wind farm siting. The change in topography gives rise to flow phenomenon such as speed-up across the hills, flow separation behind hills with high slopes and generation of localized turbulent structures. This makes flow prediction in complex terrains challenging and turbine behavior in such locations is far from understood.

Recent works on wake interaction in complex terrains either deal with the topic in ideal complex terrains such as Gaussian or sinusoidal hills or present site-specific studies. performed an experimental study of five turbines located on a Gaussian two-dimensional hill to understand the interaction of wind farm with the terrain. Their experimental setup was numerically regenerated in a large eddy simulation (LES) domain by in an attempt to validate their LES model. also studied wakes over periodic sinusoidal hills. Using experimental and numerical tools, they achieved good agreement between the two; however, implementing the Jensen wake model did not yield good results.

Coming to the site-specific studies, studied the impact of wakes and terrain on the performance of a wind farm located in southern Italy. used detached eddy simulation to study the wake from a single turbine in a complex site in Germany. The impact of topography on wake development was highlighted by comparison with a flat terrain case. Recently, Berg and coworkers performed a large eddy simulation study of wake from a wind turbine located at the site of Perdigão, Portugal , and highlighted some differences in wake orientation from Shamsoddin and Porté-Agel's work . From analyzing wind scanner data from Perdigão, show a strong influence of atmospheric stability on wake propagation and orientation. This further supports the hypothesis that it is difficult to generalize wind turbine wake results for complex sites.

The velocity deficit profiles as a function of downstream distance in the wake collapse on a single profile when scaled with the maximum deficit and characteristic wake-width; this property is known as self-similarity. Wakes behind turbines are known to show self-similar behavior in a flat terrain under different atmospheric conditions. showed mean velocity deficit behind turbines to be self-similar in neutral conditions. found self-similarity under various stability classes. This was further explored by , who then checked for self-similarity under different stability conditions in the presence of the Coriolis force. This self-similar behavior of mean velocity deficit in the far wake has been a fundamental assumption for many analytical wake models .

In the current work, we look for self-similarity of wind turbine wakes in a complex terrain under neutral atmospheric conditions and without the effect of the Coriolis force. For this purpose, we extend the work by and verify self-similarity of wakes under different terrain characteristics and turbine locations. If successful, this can potentially provide a basis for the development of analytical models for wakes in complex terrains . This self-similarity of rotor wakes in complex terrains has not yet been verified in the existing literature to the best of the authors' knowledge.

The rest of the paper is structured as follows: Sect. 2 details the LES code used for the study; terrain characteristics along with case configuration and data descriptions are provided in Sect. 3. A nomenclature for wake statistics necessary for the self-similarity check is presented in Sect. 4, and results from the study are shown in Sect. 5. The paper is finally concluded in Sect. 6.

2 LES framework

The large eddy simulation code used in the current study is formerly formulated in with examples of usage in complex terrains in . The governing equations are the spatially filtered incompressible Navier–Stokes and continuity equations, which for a neutrally stratified atmospheric flow read as

$\begin{array}{}\text{(1)}& & \frac{\partial \stackrel{\mathrm{̃}}{{u}_{i}}}{\partial {x}_{i}}=\mathrm{0},\text{(2)}& & \frac{\partial \stackrel{\mathrm{̃}}{{u}_{i}}}{\partial t}+\frac{\partial \stackrel{\mathrm{̃}}{{u}_{i}}\stackrel{\mathrm{̃}}{{u}_{j}}}{\partial {x}_{j}}=-\left[\frac{\partial {\stackrel{\mathrm{̃}}{p}}^{*}}{\partial {x}_{i}}+\frac{{f}_{i}}{\mathit{\rho }}+\frac{\partial {\mathit{\tau }}_{ij}}{\partial {x}_{j}}\right]+{F}_{\mathrm{p}}{\mathit{\delta }}_{i\mathrm{1}},\end{array}$

where $\stackrel{\mathrm{̃}}{{u}_{i}}$ is the resolved velocity in x, y and z directions corresponding to i={1, 2, 3}, respectively. The pressure variable ${\stackrel{\mathrm{̃}}{p}}^{*}$ is solved using the elliptic Poisson equation in an iterative manner. fi is the body force which models the interaction of the turbine with the flow and ρ is the fluid density, which is kept constant throughout the study. The kinematic subgrid scale stresses are represented by ${\mathit{\tau }}_{ij}=\stackrel{\mathrm{̃}}{{u}_{i}{u}_{j}}-\stackrel{\mathrm{̃}}{{u}_{i}}\stackrel{\mathrm{̃}}{{u}_{j}}$. The external forcing driving the flow is represented by Fp, which in the current case is the constant pressure gradient in the streamwise direction (see later text). An important thing to note is that the Coriolis and buoyancy forces are neglected in the current study.

A terrain-following coordinate transformation is used to represent the complex geometry. The transformation that maps physical coordinates, xi=(x, y, z), onto computational coordinates, ξi=(ξ, η, ζ), xiξi, is given by the following rules:

$\begin{array}{}\text{(3)}& & \mathit{\xi }=\mathit{\xi }\left(x\right)=x,\text{(4)}& & \mathit{\eta }=\mathit{\eta }\left(y\right)=y,\text{(5)}& & \mathit{\zeta }=\mathit{\zeta }\left(x,y,z\right),\end{array}$

with the corresponding Jacobian given by

$\begin{array}{}\text{(6)}& \mathbf{J}=\mathrm{det}\left[\begin{array}{ccc}{\mathit{\xi }}_{x}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& {\mathit{\eta }}_{y}& \mathrm{0}\\ {\mathit{\zeta }}_{x}& {\mathit{\zeta }}_{y}& {\mathit{\zeta }}_{z}\end{array}\right]={\mathit{\xi }}_{x}{\mathit{\eta }}_{y}{\mathit{\zeta }}_{z}={\mathit{\zeta }}_{z},\end{array}$

where the subscript denotes partial derivative. The governing equations, Eqs. (1) and (2), can now be transformed using the chain rule and the identity (consult , for details),

$\begin{array}{}\text{(7)}& \frac{\partial }{\partial {\mathit{\xi }}_{j}}\left(\frac{\mathrm{1}}{J}\frac{\partial {\mathit{\xi }}_{j}}{\partial {x}_{i}}\right)=\mathrm{0},\end{array}$

in order to write a set of equations in a strong flux-conservation form using the volume flux variables,

$\begin{array}{}\text{(8)}& {U}_{i}=\frac{\stackrel{\mathrm{̃}}{{u}_{j}}}{J}\frac{\partial {\mathit{\xi }}_{i}}{\partial {x}_{j}},\end{array}$

so that Ui=(U, V, W) is normal to surfaces of constant ξi. In the numerical mesh, W is located on cell faces, while U and V are located in cell centers. This ensures a tight coupling to pressure defined at the cell centers. Solving the Poisson equation involves an iteration procedure which incurs additional expense compared to conventional computations in flat terrain (again we recommend the reader to consult , for details).

The code is pseudo-spectral code, with wave number representation in the horizontal directions (ξ, η), and 2nd-order finite difference in the vertical direction, ζ.

The physical (z) and computational (ζ) vertical coordinates are related to each other using a simple transformation rule, with z exponentially stretched from the surface:

$\begin{array}{}\text{(9)}& \mathit{\zeta }=H\frac{z-h}{H-h},\end{array}$

where h=h(x, y) is the local terrain height and H is the domain height.

The Deardorff SGS (subgrid scale) model as adopted in is implemented to parameterize the subgrid scale stresses.

## Turbine parameterization

An actuator disk model without any rotational effects is implemented to represent the turbines. The unavailability of a well defined free stream velocity, U, in complex terrain is compensated by employing the classical expression given by and introduced into actuator disk representation in LES by :

$\begin{array}{}\text{(10)}& {U}_{\mathrm{\infty }}=\frac{{\stackrel{\mathrm{‾}}{u}}_{\mathrm{d}}}{\mathrm{1}-a},\end{array}$

where ${\stackrel{\mathrm{‾}}{u}}_{\mathrm{d}}$ is the velocity averaged over the rotor disk and a is the axial induction factor. The model simulates turbines using disks with a thrust force given by

$\begin{array}{}\text{(11)}& {F}_{\mathrm{T}}=-\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{C}_{\mathrm{T}}^{\prime }{\stackrel{\mathrm{‾}}{u}}_{\mathrm{d}}^{\mathrm{2}}\frac{\mathit{\pi }}{\mathrm{4}}{D}^{\mathrm{2}},\end{array}$

where D is the rotor diameter and ${C}_{\mathrm{T}}^{\prime }$ is related to the thrust coefficient CT using one-dimensional momentum theory:

$\begin{array}{}\text{(12)}& {C}_{\mathrm{T}}^{\prime }=\frac{{C}_{\mathrm{T}}}{\left(\mathrm{1}-a{\right)}^{\mathrm{2}}}=\frac{\mathrm{4}a}{\left(\mathrm{1}-a\right)}.\end{array}$

For the current work, we choose an induction, $a=\mathrm{1}/\mathrm{4}$, which corresponds to ${C}_{\mathrm{T}}^{\prime }=\mathrm{4}/\mathrm{3}$. In the numerical code, the thrust force is then distributed proportionally to the fractional area of each cell covered by the rotor. The force per unit volume in the longitudinal direction is then given by

$\begin{array}{}\text{(13)}& {f}_{x}\left(i,j,k,t\right)=-\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{C}_{\mathrm{T}}^{\prime }{\stackrel{\mathrm{‾}}{u}}_{\mathrm{d}}\left(t{\right)}^{\mathrm{2}}\frac{\mathrm{\Lambda }\left(j,k\right)}{\mathrm{\Delta }x},\end{array}$

where Λ is the fractional area and Δx is length of a given cell. In the implementation, we have added a 10 min exponential time filter to ${\stackrel{\mathrm{‾}}{u}}_{\mathrm{d}}$.

3 Case configurations

Simulations are performed on a domain spanning over $\mathrm{5120}\phantom{\rule{0.125em}{0ex}}\mathrm{m}×\mathrm{2560}\phantom{\rule{0.125em}{0ex}}\mathrm{m}×\mathrm{3000}$ m, which represents a double-ridge-configured site named Perdigão in Portugal . The site has been the focus of an intensive field measurement campaign under the Perdigão experiment 2017 and the New European Wind Atlas (NEWA) project . Fig. 1a shows the overview of the site, where the wind direction is along the main transect (235 relative to north). With its double ridge configuration and challenging slopes, the terrain is ideal for a comprehensive study of interaction between wind turbine wake and topography. In Fig. 1b we show the numerical terrain. In order to make the terrain periodic to comply with the pseudo-spectral LES model, a little smoothing had to be performed. The atmosphere is assumed to be neutrally stratified and boundary conditions used are described below. A roughness length of z0=0.5 m is chosen, keeping the rugged, forested terrain in mind, giving ${z}_{\mathrm{0}}/H=\mathrm{1.67}×{\mathrm{10}}^{-\mathrm{4}}$. The simulations are initiated from a random incompressible velocity field and run for a time of TS=100TE until stationarity was achieved. Here, ${T}_{\mathrm{E}}=H/{u}_{*}$ (u* is defined as friction velocity) is the timescale corresponding to the size of the largest eddy that can possibly fit in the domain and is based on a simple momentum transfer argument borrowed from the flat terrain. Figure 2 shows an instantaneous snapshot of streamwise velocity along the main transect passing through the turbine located on top of the first ridge with respect to flow direction. The turbine is assumed to have a rotor diameter and hub height of 80 m each. This is done to somewhat match the dimensions of an ENERCON 2 MW turbine actually deployed at the upwind ridge at Perdigão.

Figure 1Terrain elevation variation at Perdigão. (a) Shaded rectangle represents the simulated area, with a yellow line passing through the main transect of the turbine location (source: ); (b) shaded region from (a) with intersection of red lines representing the turbine locations.

Figure 2Instantaneous velocity in the streamwise direction along the main transect with turbine (in red).

The iteration procedure performed when solving the Poisson equation is expensive due to the strong coupling between pressure and vertical wind component and the fact that, in contrast to flat terrain, pressure gradients are quite high in complex terrains and thus require more time to fully converge. Slight smoothing of the terrain can help in reducing the computation time by reducing the number of pressure iterations required to attain convergence. This is done by applying an exponential filter in the wave number (x and y) space. This filtering process is given as follows:

$\begin{array}{}\text{(14)}& \begin{array}{rl}{h}^{\mathrm{model}}\left(x,y\right)=& \int \int {h}^{\mathrm{true}}\left({k}_{x},{k}_{y}\right)\stackrel{\mathrm{̃}}{F}\left({k}_{x},{k}_{y}\right)\\ & {e}^{i\left({k}_{x}x+{k}_{y}y\right)}\mathrm{d}{k}_{x}\mathrm{d}{k}_{y},\end{array}\end{array}$

where $\stackrel{\mathrm{̃}}{F}\left({k}_{x}$, ky) is the exponential filter given by

$\begin{array}{}\text{(15)}& \stackrel{\mathrm{̃}}{F}\left({k}_{x},{k}_{y}\right)={e}^{-\mathit{\alpha }\left[{\left({k}_{x}\mathrm{\Delta }x\right)}^{\mathrm{2}}+{\left({k}_{y}\mathrm{\Delta }y\right)}^{\mathrm{2}}\right]}.\end{array}$

Here, α is the control parameter which determines the level of smoothing. Two different values of α used in the current study are 4 and 0.5, which result in maximum terrain slopes of 0.57 and 0.77 across the main transect, respectively. The steep terrain nearly matches the original terrain. The α value of 0.5 instead of 0 is chosen for this terrain to avoid the Gibbs phenomenon. Figure 3 shows a comparison of two terrains at a fixed lateral position passing through the turbine location. From here onward, a terrain with maximum slope of 0.57 will be referred to as “smooth”, whereas one with maximum slope of 0.77 will be referred to as “steep”.

Figure 3Comparison of terrains with two different slopes along the main transect.

Table 1 gives an overview of various cases analyzed in the current study. It is to be noted that the number of ensembles used is fairly high to guarantee stationarity and an apparent difference in the number of ensembles for different cases does not affect the results. The aspect ratio of 1:4 in the case of 20 m resolution lies toward the limit of the allowed aspect ratio for the chosen subgrid scale model (rule of thumb).

Table 1Description of cases analyzed in the study.

## 3.1 Boundary conditions

The boundaries in the horizontal directions are fully periodic, where a constant pressure gradient given by dPd$x=-{u}_{*}^{+\mathrm{2}}/H$ is applied in the x direction to drive the flow. Here, ${u}_{*}^{+}$ is the friction velocity chosen to be 0.6 m s−1 and H is the height of domain equal to 3000 m. It is worth mentioning here that due to complexity of the terrain, the value used for friction velocity is not the effective friction velocity in the domain (discussed later). The periodicity in the horizontal directions means that the terrain is repeating itself, i.e., the inflow is affected by the terrain in the domain itself. We have made tests with extended buffer regions (not shown) which indicated that the turbulence level in the incoming flow was not affected by the terrain itself, and we therefore do not consider the periodicity to be a shortcoming of the later conclusions.

The lower surface of the domain is modeled assuming a log-law behavior in the first cell as a pointwise implementation . A no-stress condition is applied at the upper boundary of the domain. Under this condition, the gradients of horizontal velocity components are set to zero, and the vertical velocity component itself is set to zero. Moreover, the subgrid scale energy (e) is also set to zero at upper boundary.

## 3.2 Data description

We base our analysis on averaged LES data, a choice associated with the need of stationary flow characteristics required for evaluation of self-similar wake statistics. We save 30 min averaged LES fields and by averaging a number of these fields together a total time necessary to guarantee stationarity is obtained. The rule for this averaging is defined as follows:

$\begin{array}{}\text{(16)}& 〈u〉=\frac{\mathrm{1}}{N}\sum _{n=\mathrm{1}}^{N}{u}_{\left(\mathrm{30},n\right)},\end{array}$

where u is the three-dimensional ensemble averaged velocity field, u(30,n) is the nth 30 min averaged field and N is the total number of 30 min averages used, as given in Table 1 for each case. To simplify the notations, mean streamwise velocity computed using the above rule will be denoted by u=〈u.

4 Wake flow nomenclature

The quantities required to characterize the wake flow for the self-similarity check are defined in the current section. We first define the normalized mean velocity deficit as follows:

$\begin{array}{}\text{(17)}& \mathrm{\Delta }U=\frac{{u}_{\mathrm{1}}-{u}_{\mathrm{0}}}{{u}_{\mathrm{ref}}},\end{array}$

where u1 is the averaged streamwise velocity in the simulation including the turbine and u0 is the same in the absence of the turbine. The mean difference is then normalized with a reference velocity, which corresponds to the velocity at the hub height of the first turbine in the absence of the turbine. Table 2 lists the reference velocities used for different cases.

Table 2Reference velocities at hub height used for different cases.

Wake half-width is defined as the distance from wake center to the point where the velocity deficit is reduced to 50 % of the maximum velocity deficit. As wind turbine wakes can be asymmetric, the wake half-width is computed for two sides of the wake individually. The nomenclature used when describing the wake characteristics is defined in Fig. 4. The wake centerline is traced by identifying the point of maximum velocity deficit at each downstream location.

Figure 4Nomenclature for wake half-width.

5 Results

## 5.1 Inflow velocity

To characterize the incoming flow velocity, normalized mean velocity profiles at 20 rotor diameters (20 D) upstream of the first ridge for different cases along with a logarithmic profile are shown in Fig. 5. Here, uh refers to the velocity at the hub height approximately 20 rotor diameters upstream in the streamwise direction. It can be observed that the simulated velocity profiles show some deviation from the logarithmic profile, especially in the case of complex terrain. These deviations are somewhat expected, as the logarithmic profile is not completely valid in heterogeneous terrains. Moreover, the deviations observed are due to differences in the slopes of various profiles, whereas the quantitative values lie very close to the logarithmic profile. The periodic boundary conditions also play a role in determining the inflow velocity and therefore could be responsible for some of the deviations from the logarithmic profile.

Figure 5Comparison of normalized mean streamwise velocity at x=20 m. Smooth and flat cases have a resolution of 20 m, whereas the steep case has a resolution of 10 m.

## 5.2 Impact of grid resolution

The steep slope case is simulated at two different horizontal grid resolutions of 20 and 10 m. It should be mentioned that our procedure includes contributions from two effects: the grid resolution itself but also terrain effects, since more terrain is resolved in the 10 m case compared to the 20 m case. Figure 6 shows vertical profiles of streamwise velocity across the turbine located on top of the first ridge for the two grid resolutions. Comparing the velocity profiles in the two cases, good agreement is observed at respective downstream locations. However, it is important to note that the agreement between velocity profiles for the two cases is sensitive to the chosen reference velocity. Figure 7 shows a comparison of streamwise velocity along the main transect at the hub height, as well as local slopes, for the two grid resolutions. Whereas the qualitative trend of the velocity profiles matches well for the two resolutions, the quantitative values differ by up to 15 % between the two cases. This apparent discrepancy can be somewhat justified by the change in terrain characteristics caused by the change in resolution. As the resolution gets finer, the terrain becomes more detailed and the slopes are slightly changed. This can be verified by comparing the relative vertical distance between local surface and rotor positions up- and downstream of the ridge.

Figure 6Vertical profiles of normalized mean streamwise velocity around first ridge for two different grid resolutions. Blue and yellow vertical lines show velocity profiles without and with turbine for 20 m resolution, whereas red and purple vertical lines show profiles without and with turbine for 10 m resolution. Black horizontal lines trace rotor top and bottom for the coarse grid resolution, whereas cyan horizontal lines trace the same for finer-grid resolution. The error bars correspond to $\mathit{\sigma }/\sqrt{N}$, where σ is the standard deviation and N is the number of 30 min ensembles used.

Figure 7Comparison of streamwise velocity at the hub height without turbine and local slopes along the main transect for the two grid resolutions.

Figure 8(a) Streamwise velocity profile for smooth and steep terrain; (b) wake development in the rotor plane at different downstream locations. First row: first turbine, smooth terrain – 20 m resolution; second row: first turbine, steep terrain – 10 m resolution; third row: second turbine, steep terrain – 10 m resolution. Rotor position is highlighted by the red circle.

Recall that in Sect. 3.1 the flow was defined to be pressure-gradient-driven with a friction velocity of 0.6 m s−1. This balance between pressure gradient and friction velocity would hold true in flat terrain, whereas in complex terrain there is an additional contribution from the form drag generated by the two ridges. The immediate impact of this additional contribution is a reduced friction velocity in the terrain. As the terrain gets more detailed in the finer grid, the friction velocity is further decreased by virtue of an increased contribution from form drag. We believe that this difference in friction velocity in the two grid resolutions can explain the difference in the velocity profiles observed here.

A comparison of the LES flow characteristics in different terrains and resolutions with field measurements, although without any turbine in the terrain, was recently presented by . Good agreement between the measurements and simulations was obtained, keeping in mind the simple flow conditions defined in the simulations.

Figure 9Lateral profiles of streamwise velocity deficit for different cases.

Figure 10Vertical profiles of streamwise velocity deficit for different cases.

## 5.3 Wake characteristics

A fundamental question can be formulated around the interaction of wake from the turbine with the local terrain. Understanding this kind of interaction is important in order to know how the wake behaves within a complex terrain and how the terrain-generated phenomena impact the behavior of the rotor wake. Figure 8a shows the wakes from turbines in two different terrains and two different locations (top of the respective ridges). An important impact of local topography on wakes is the change in the orientation of wakes behind the two turbines. The strong recirculation regions developed behind the ridges due to flow separation generate their own shear layers, which do not allow the turbine wakes to mix with them. It is for this reason that we observe an upward orientation for the turbines located on top of the first ridge in the two cases. For the smooth terrain, the recirculation starts on top of the ridge and the height of this region extends higher than the ridge top. This causes an initially horizontal wake to shift in an upward direction. For the steep case, however, the high slopes result in an upward inclination of the wake from the rotor position. For the turbine located on the second ridge, the wake is oriented in a slightly downward, somewhat terrain-following, orientation. This is due to the fact that the recirculation region behind the second ridge is smaller than the ridge top and thus allows the wake to move in a downward direction.

Figure 8b shows wake development in the rotor plane. The sideways movement of wake can be observed in the three cases, as the wake moves downstream. Although the most probable reason for this sideways movement is the terrain complexity and variation in pressure distributions across the terrain, any conclusive statement cannot be made due to the limitations of grid resolution employed in the current study.

An important feature of the wake is the relatively faster recovery as compared to the flat terrain. High levels of turbulent mixing in the atmosphere due to complexity of the terrain provide a catalyst effect and thus promote quick recovery of the wake. It is thereby observed that the wake is almost recovered at a distance of around 5 rotor diameters in the downstream. This faster wake recovery in the complex terrain has been previously reported by various other studies .

### Velocity deficit profiles

Figures 9 and 10 show normalized mean velocity deficit profiles for different cases (flat, smooth, steep – 10 m and first turbine; steep – 10 m and second turbine) at different downstream locations. Comparing the lateral profiles, it is observed that the finer-grid resolution in the lateral direction captures the wake structure much better than the coarser resolutions. Moreover, the velocity deficit decays at a much slower rate in the flat case than the other (complex) cases. The velocity deficit for the turbine located on the second ridge is observed to be the smallest with the fastest recovery. This can be attributed to the highest ambient turbulence in the wake of the particular turbine among the considered cases. To support this argument, we plot normalized profiles of turbulent kinetic energy, TKE =$\frac{〈{\mathit{\sigma }}_{u}^{\mathrm{2}}+{\mathit{\sigma }}_{v}^{\mathrm{2}}+{\mathit{\sigma }}_{w}^{\mathrm{2}}〉}{\mathrm{2}}$ (see Fig. 11). As turbulent kinetic energy is responsible for transport of energy in the domain, it can be used to quantify the wake recovery. From the figure it is clearly observed that the turbulent kinetic energy for complex terrain cases is significantly higher than the flat terrain case. The turbine located on top of the second ridge in the steep case shows the highest levels of TKE, thus causing fastest wake recovery.

Figure 11Vertical profiles of normalized turbulent kinetic energy for different cases.

## 5.4 Wake self-similarity

Figure 12 shows how the maximum velocity deficit develops as the wake moves in the downstream. An increasing–decreasing behavior can be observed, with the smooth case showing the highest velocity deficit. At a downstream distance of 5 rotor diameters, the three complex cases show very close numbers for maximum velocity deficit. The flat case (as expected) shows slowest recovery in wake.

Figure 12Maximum normalized mean velocity deficit as a function of downstream distance for different cases.

To further characterize the wake, we plot right and left wake half-widths for the lateral and vertical profiles in Fig. 13. The wake half-widths generally show an increasing trend with the downstream distance. For the lateral profiles the values are much closer near the turbine location and spread differently for different cases, whereas for vertical profiles the values have a wider spread due to varying levels of shear. The flat terrain is observed to have the smallest wake half-widths, whereas the second turbine in steep slope shows highest values. This can be attributed to the rate of wake recovery in the respective cases. Finally, the difference in the left and right wake widths for a specific location highlights the wake symmetry or asymmetry.

Figure 13Right and left wake half-width for the (a) lateral and (b) vertical velocity deficit profiles.

We now plot normalized velocity deficit profiles against normalized distance from wake center. The velocity deficit profiles are normalized with the maximum velocity deficit for the respective profiles, whereas the centered distance is normalized with wake half-width on either side of the wake center.

Figure 14a shows lateral self-similar velocity deficit profiles for the flat case. The profiles follow the Gaussian profile and last from a downstream distance of 1 to 9 rotor diameters. The same profiles for three complex cases shown in Fig. 14b also collapse on the Gaussian profile; however, two important things should be noted. First the self-similarity holds for a much shorter distance (from 1 to 3 rotor diameters). This can be attributed to the faster wake recovery in complex cases and to the fact that the localized terrain changes disrupt the wake structure and thus deviations from self-similarity occur. The second important feature is the wide spread in the tails of the profiles for three complex cases. Whereas the profiles for a single case (denoted by same color) are much closer to each other, wider spread is seen while comparing profiles for two different cases. This can be due to the difference in the shear toward the edges of the wakes in different cases and to the differences arising due to differences in terrains, flow separation characteristics and levels of turbulence. Now, Fig. 14b shows velocity deficit profiles extracted from a horizontal plane passing through the hub height of the respective turbine; however, from Fig. 8, we observe significant deflection in the vertical direction. In order to account for this wake deflection, we track the maximum velocity deficit point in the vertical direction at each downstream location and extract the corresponding lateral velocity profiles. Figure 15 shows the resulting self-similar profiles of streamwise velocity deficit from 1 to 6 D. Computing self-similar profiles beyond 6 D is not possible, as the wake is almost recovered (the max velocity deficit drops below 0.04 for >6 D). Comparing with Fig. 14b, the profiles are observed to be self-similar for a longer distance; however, beyond 3.5 D, the left-side tail of the profiles depart from the Gaussian profile with increasing downstream distance.

Figure 14Self-similar lateral profiles of streamwise velocity deficit (a) flat case (b) complex cases for a downstream distance of 1 D to 3 D with intervals of 0.5 D; inset shows normalized velocity deficit profiles from 3.5 D to 5 D.

Figure 15Self-similar lateral profiles of streamwise velocity deficit for complex cases accounting for wake deflection in the vertical direction for a downstream distance of (a) 1 D to 3.5 D and (b) 4 D to 6 D.

Figure 16Self-similar vertical velocity deficit profiles (a) flat case (b) complex cases for a downstream distance of 1 D to 5 D with intervals of 0.5 D. Arrows indicate direction of increasing downstream distance.

Figure 16a shows vertical self-similar velocity deficit profiles for the flat case. Self-similarity is observed from a distance of 3 to 10 rotor diameters. In the complex cases (Fig. 16b), self-similarity in the vertical profiles is observed from 1 to 3 rotor diameters. As indicated by arrows in Fig. 16b, the spread in tails of profiles increases with increasing downstream distance. Beyond 3 D, this spread is significantly high. These profiles also show a much wider spread when compared to the lateral profiles. Moreover, they exhibit more asymmetric behavior than the lateral profiles. This is due to the impact of surface and vertical wind shear. In addition, as the terrain is complex, varying levels of slope lead to varying distances between the projected rotor area and the surface, which eventually gives rise to high variation among profiles at different downstream locations.

6 Conclusion

Terrain-generated phenomena can have a significant impact on the wake characteristics of a turbine located in a complex terrain. In this study, we attempted to verify whether the self-similarity, which is one the fundamental characteristics of wakes, still holds in complex terrains or not. In this context, we performed large eddy simulation of wind turbine(s) located in a flat and complex terrain. By varying turbine location, as well as terrain complexity, we simulated different flow scenarios in complex terrain. By plotting normalized velocity deficit profiles in the lateral, as well as vertical, direction, we looked for self-similarity in the simulated cases. We observed that wakes in complex terrain preserve self-similarity in both directions. The region of this preservation is, however, over much shorter distances downstream than the flat terrain counterpart. This is attributed to the deformation of wakes in the far wake region by varying heights of the terrain and to the faster wake recovery in the complex terrain. Finally, vertical profiles show a wider spread and higher asymmetry than the lateral ones.

Data availability
Data availability.

The full data set of the LES simulations is on the order of hundreds of gigabytes. By contacting the authors a smaller subset can be made available.

Author contributions
Author contributions.

ASD performed the research work and prepared the article. JB did the simulations, conceived the research plan and supervised the research work and the article preparation. NT supervised the research work. EGP took part in developing the LES code and guided aspects of the analysis. All authors contributed to the improvement and preparation of the final article.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Special issue statement
Special issue statement.

Acknowledgements
Acknowledgements.

The authors would like to thank Niels N. Sørensen, DTU Wind Energy, for help in preparing the surface meshes and Peter P. Sullivan, National Center for Atmospheric Research, Boulder, CO, US, for sharing the LES code.

Review statement
Review statement.

This paper was edited by Raúl Bayoán Cal and reviewed by Dries Allaerts and one anonymous referee.

References

Abkar, M. and Porté-Agel, F.: Influence of atmospheric stability on wind-turbine wakes: A large-eddy simulation study, Phys. Fluids, 27, 035104, https://doi.org/10.1063/1.4913695, 2015. a, b

Abkar, M., Sharifi, A., and Porté-Agel, F.: Wake flow in a wind farm during a diurnal cycle, J. Turbulence, 17, 420–441, 2016. a

Alfredsson, P. H. and Segalini, A.: Introduction Wind farms in complex terrains: an introduction, Philos. T. Roy. Soc. A , 375, 20160096, https://doi.org/10.1098/rsta.2016.0096, 2017. a

Allaerts, D. and Meyers, J.: Large eddy simulation of a large wind-turbine array in a conventionally neutral atmospheric boundary layer, Phys. Fluids, 27, 065108, https://doi.org/10.1063/1.4922339, 2015. a

Astolfi, D., Castellani, F., and Terzi, L.: A Study of Wind Turbine Wakes in Complex Terrain Through RANS Simulation and SCADA Data, J. Sol. Energy Eng., 140, 031001, https://doi.org/10.1115/1.4039093, 2018. a

Barthelmie, R. J., Frandsen, S. T., Nielsen, M. N., Pryor, S. C., Rethore, P. E., and Jørgensen, H. E.: Modelling and measurements of power losses and turbulence intensity in wind turbine wakes at Middelgrunden offshore wind farm, Wind Energy, 10, 517–528, 2007. a

Bastankhah, M. and Porté-Agel, F.: A new analytical model for wind-turbine wakes, Renew. Energy, 70, 116–123, 2014. a

Berg, J., Troldborg, N., Sørensen, N. N., Patton, E. G., and Sullivan, P. P.: Large-Eddy Simulation of turbine wake in complex terrain, J. Phys.: Conf. Ser., 854, 012003, https://doi.org/10.1088/1742-6596/854/1/012003, 2017. a, b, c

Berg, J., Troldborg, N., Menke, R., Patton, E. G., Sullivan, P. P., Mann, J., and Sørensen, N. N.: Flow in complex terrain-a Large Eddy Simulation comparison study, J. Phys.: Conf. Ser., 1037, 072015, https://doi.org/10.1088/1742-6596/1037/7/072015, 2018. a

Bou-Zeid, E., Meneveau, C., and Parlange, M. B.: Large-eddy simulation of neutral atmospheric boundary layer flow over heterogeneous surfaces: Blending height and effective surface roughness, Water Resour. Res., 40, W02505, https://doi.org/10.1029/2003WR002475, 2004. a

Calaf, M., Meneveau, C., and Meyers, J.: Large eddy simulation study of fully developed wind-turbine array boundary layers, Phys. Fluids, 22, 1–16, 2010. a, b

Castellani, F., Astolfi, D., Mana, M., Piccioni, E., Becchetti, M., and Terzi, L.: Investigation of terrain and wake effects on the performance of wind farms in complex terrain using numerical and experimental data, Wind Energy, 20, 1277–1289, 2017. a

Chamorro, L. P. and Porté-Agel, F.: A wind-tunnel investigation of wind-turbine wakes: Boundary-Layer turbulence effects, Bound.-Lay. Meteorol., 132, 129–149, 2009. a

Feng, J., Shen, W. Z., Hansen, K. S., Vignaroli, A., Bechmann, A., Zhu, W. J., Larsen, G. C., Ott, S., Nielsen, M., and M. M. Jogararu: Wind farm design in complex terrain: the FarmOpt methodology, Paper presented at China Wind Power 2017, October 2017, Beijing, China, 2017. a

Hansen, M. O. L.: Aerodynamics of wind turbines, Routledge, New York, 2015. a

Hyvarinen, A. and Segalini, A.: Qualitative analysis of wind-turbine wakes over hilly terrain, J. Phys.: Conf. Ser., 854, 012023, https://doi.org/10.1088/1742-6596/854/1/012023, 2017. a

Iungo, G. V.: Experimental characterization of wind turbine wakes: Wind tunnel tests and wind LiDAR measurements, J. Wind Eng. Indust. Aerodynam., 149, 35–39, 2016. a

Iungo, G. V., Wu, Y. T., and Porté-Agel, F.: Field measurements of wind turbine wakes with lidars, J. Atmos. Ocean. Tech., 30, 274–287, 2013. a

Jimenez, A., Crespo, A., Migoya, E., and Garcia, J.: Advances in large-eddy simulation of a wind turbine wake, J. Phys.: Conf. Ser., 75, 012041, https://doi.org/10.1088/1742-6596/75/1/012041, 2007. a

Katic, I., Højstrup, J., and Jensen, N. O.: A Simple Model for Cluster Efficiency, in: European Wind Energy Association Conference and Exhibition, 7–9 October 1986, Rome, Italy, 407–410, 1986. a, b

Lutz, T., Schulz, C., Letzgus, P., and Rettenmeier, A.: Impact of Complex Orography on Wake Development: Simulation Results for the Planned WindForS Test Site, J. Phys.: Conf. Ser., 854, 012029, https://doi.org/10.1088/1742-6596/854/1/012029, 2017. a

Luzzatto-Fegiz, P.: A one-parameter model for turbine wakes from the entrainment hypothesis, J. Phys.: Conf. Ser., 1037, 072019, https://doi.org/10.1088/1742-6596/1037/7/072019, 2018. a

Mann, J., Angelou, N., Arnqvist, J., Callies, D., Cantero, E., Chavez Arroyo, R., Courtney, M., Cuxart, J., Dellwik, E., Gottschall, J., Ivanell, S., Kühn, P., Lea, G., Matos, J. C., Palma, J. M. L. M., Pauscher, L., Peña, A., Sanz Rodrigo, J., Söderberg, S., Vasiljević, N., and Veiga Rodrigues, C.: Complex terrain experiments in the New European Wind Atlas, Philos. T. Roy. Soc. A, 375, 20160101, https://doi.org/10.1098/rsta.2016.0101, 2017. a

Medici, D. and Alfredsson, P. H.: Measurements on a wind turbine wake: 3D effects and bluff body vortex shedding, Wind Energy, 9, 219–236, 2006. a

Menke, R., Vasiljević, N., Hansen, K. S., Hahmann, A. N., and Mann, J.: Does the wind turbine wake follow the topography? A multi-lidar study in complex terrain, Wind Energ. Sci., 3, 681–691, https://doi.org/10.5194/wes-3-681-2018, 2018.  a, b

Moeng, C.-H.: A large-eddy-simulation model for the study of planetary boundary-layer turbulence, J. Atmos. Sci., 41, 2052–2062, 1984. a

Politis, E. S., Prospathopoulos, J., Cabezon, D., Hansen, K. S., Chaviaropoulos, P. K., and Barthelmie, R. J.: Modeling wake effects in large wind farms in complex terrain: The problem, the methods and the issues, Wind Energy, 15, 161–182, 2012. a

Porté-Agel, F., Yu-Ting Wu, Lu, H., and Conzemius, R. J.: Large-eddy simulation of atmospheric boundary layer flow through wind turbines and wind farms, J. Wind Eng. Indust. Aerodynam., 99, 154–168, 2011. a

Shamsoddin, S. and Porté-Agel, F.: Large-Eddy Simulation of Atmospheric Boundary-Layer Flow Through a Wind Farm Sited on Topography, Bound.-Lay. Meteorol., 163, 1–17, https://doi.org/10.1007/s10546-016-0216-z, 2017. a, b

Sullivan, P. P., McWilliams, J. C., and Moeng, C.-H.: A subgrid-scale model for large-eddy simulation of planetary boundary-layer flows, Bound.-Lay. Meteorol., 71, 247–276, 1994. a

Sullivan, P. P., Patton, E. G., and Ayotte, K. W.: Turbulent flow over sinusoidal bumps and hills derived from large eddy simulations, in: 19th Symposium on Boundary Layers and Turbulence, August 2010, Keystone, CO, USA, 2010. a

Sullivan, P. P., McWilliams, J. C., and Patton, E. G.: Large-Eddy Simulation of Marine Atmospheric Boundary Layers above a Spectrum of Moving Waves, J. Atmos. Sci., 71, 4001–4027, 2014. a, b, c

Tabib, M., Rasheed, A., and Fuchs, F.: Analyzing complex wake-terrain interactions and its implications on wind-farm performance, J. Phys.: Conf. Ser., 753, 032063, https://doi.org/10.1088/1742-6596/753/3/032063, 2016. a

Tian, W., Ozbay, A., Yuan, W., Sarakar, P., Hu, H., and Yuan, W.: An experimental study on the performances of wind turbines over complex terrain, in: 51st AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition 2013, 7–10 January 2013, Grapevine, Texas, USA, 2013. a

Vasiljević, N. L. M. Palma, J. M., Angelou, N., Carlos Matos, J., Menke, R., Lea, G., Mann, J., Courtney, M., Frölen Ribeiro, L., and Gomes, V. M. M. G. C.: Perdigão 2015: methodology for atmospheric multi-Doppler lidar experiments, Atmos. Meas. Tech., 10, 3463–3483, https://doi.org/10.5194/amt-10-3463-2017, 2017. a

Wildmann, N., Vasiljevic, N., and Gerz, T.: Wind turbine wake measurements with automatically adjusting scanning trajectories in a multi-Doppler lidar setup, Atmos. Meas. Tech., 11, 3801–3814, https://doi.org/10.5194/amt-11-3801-2018, 2018. a

Xie, S. and Archer, C.: Self-similarity and turbulence characteristics of wind turbine wakes via large-eddy simulation, Wind Energy, 17, 657–669, 2014. a, b

Xie, S. and Archer, C. L.: A Numerical Study of Wind-Turbine Wakes for Three Atmospheric Stability Conditions, Bound.-Lay. Meteorol., 165, 87–112, 2017. a