Journal cover Journal topic
Wind Energy Science The interactive open-access journal of the European Academy of Wind Energy
Journal topic
WES | Articles | Volume 6, issue 2
Wind Energ. Sci., 6, 555–570, 2021
Wind Energ. Sci., 6, 555–570, 2021

Research article 22 Apr 2021

Research article | 22 Apr 2021

The curled wake model: a three-dimensional and extremely fast steady-state wake solver for wind plant flows

The curled wake model: a three-dimensional and extremely fast steady-state wake solver for wind plant flows
Luis A. Martínez-Tossas1, Jennifer King1, Eliot Quon1, Christopher J. Bay1, Rafael Mudafort1, Nicholas Hamilton1, Michael F. Howland2, and Paul A. Fleming1 Luis A. Martínez-Tossas et al.
  • 1National Renewable Energy Laboratory, Golden, CO, USA
  • 2Graduate Aerospace Laboratories (GALCIT), California Institute of Technology, Pasadena, CA 91125, USA

Correspondence: Luis A. Martínez-Tossas (


Wind turbine wake models typically require approximations, such as wake superposition and deflection models, to accurately describe wake physics. However, capturing the phenomena of interest, such as the curled wake and interaction of multiple wakes, in wind power plant flows comes with an increased computational cost. To address this, we propose a new hybrid method that uses analytical solutions with an approximate form of the Reynolds-averaged Navier–Stokes equations to solve the time-averaged flow over a wind plant. We compare results from the solver to supervisory control and data acquisition data from the Lillgrund wind plant obtaining wake model predictions which are generally within 1 standard deviation of the mean power data. We perform simulations of flow over the Columbia River Gorge to demonstrate the capabilities of the model in complex terrain. We also apply the solver to a case with wake steering, which agreed well with large-eddy simulations. This new solver reduces the time – and therefore the related cost – it takes to simulate a steady-state wind plant flow (on the order of seconds using one core). Because the model is computationally efficient, it can also be used for different applications including wake steering for wind power plants and layout optimization.

1 Introduction

In this work, we present an improved formulation of the curled wake model (Martínez-Tossas et al.2019) that can be used in the context of a wind power plant without the need to use a wake superposition method. Wake superposition models are typically used because of their computational efficiency; however, they have been shown to give different results depending on the model used (Gunn et al.2016; Zong and Porté-Agel2020; Hamilton et al.2020). This inconsistency motivates the use of a more robust solver in the context of the curled wake model (Martínez-Tossas et al.2019) that does not depend on a superposition method. The new solver is developed by simplifying the Reynolds-averaged Navier–Stokes (RANS) equations to obtain a parabolic equation for the wake deficit. The equation is solved in a three-dimensional domain to obtain the wake velocity in a wind plant. This solver uses a hybrid RANS-analytical framework that aims to minimize computational cost while still preserving physics from the RANS equations.

Parabolic solvers for RANS equations are a promising tool for fast wind farm flow solvers. Ainslie (1988) developed a parabolic solver for an approximation of RANS equations in cylindrical coordinates. They proposed a mixing-length eddy-viscosity model that has a component from the ambient turbulence and another from the wake-added turbulence. Iungo et al. (2018) developed a parabolic RANS solver focused on improving the mixing-length model and used assumptions about axisymmetry in the wakes. Bradstock and Schlez (2020) have developed a parabolic wind plant RANS solver (WakeBlaster) used for commercial applications. WakeBlaster solves a simplified version of the RANS equations and has been validated using field experiments (Bradstock and Schlez2020). WakeBlaster uses a special method to solve the spanwise velocity components that does not include effects caused by yaw.

Wake steering is a promising wind plant control strategy used to maximize the power output of a wind plant (Adaramola and Krogstad2011; Park et al.2013; Gebraad et al.2016; Howland et al.2019; Fleming et al.2019; Siemens Gamesa2019). In wake steering, upstream turbines are yawed, deflecting the wakes such that downstream turbines are able to produce more power, and the wind plant as a whole can produce more power. In this work, we present a wind plant model that uses a simplified version of the RANS equations to predict the flow through a wind plant with wake steering. This tool is extremely fast (order of seconds), thereby enabling control-oriented frameworks used for wind plant operation and layout optimization.

The wake of a wind turbine in yaw has a unique shape known as the curled wake (Howland et al.2016; Bastankhah and Porté-Agel2016; Martínez-Tossas et al.2019). This shape has been observed in computational fluid dynamics simulations and in small- and large-scale experiments (Medici and Alfredsson2006; Howland et al.2016; Bastankhah and Porté-Agel2016; Vollmer et al.2016; Bartl et al.2018; Fleming et al.2018b; Schottler et al.2018). The curled wake is formed because the wake of a wind turbine in yaw introduces spanwise and vertical velocities that deform the wake and change its shape. This mechanism has been explained in the literature as a collection of vortices shed from the rotor plane (Howland et al.2016; Bastankhah and Porté-Agel2016; Shapiro et al.2018; Martínez-Tossas et al.2019). The curled wake is a unique phenomenon in wind turbine wakes because it disrupts the asymmetry of the wake. As a result of the shed vortices, the curled yawed wind turbine wake is laterally asymmetric and non-Gaussian. The curled wake is known to affect not only a turbine immediately downstream, but also subsequent turbines within a wind plant. This effect is known as secondary steering and it is important to capture it when using wake models to unravel the full potential of wake steering (Fleming et al.2018a; Bay et al.2020; King et al.2020).

The curled wake model uses a simplified version of the RANS equations to predict the wake of a wind turbine in yaw (Martínez-Tossas et al.2019). Several improvements have been proposed to the original formulation of the model, including (1) a decay model for the vortices, (2) tuning of the viscous term based on turbulence intensity, and (3) adding a pressure gradient term to account for wake expansion (Bay et al.2019, 2020; Hulsman et al.2020). Also, the new Gauss-curl hybrid model has been shown to provide a good compromise between an analytical model (Bastankhah and Porté-Agel2016) and some of the physics from the curled wake model (King et al.2020). The original formulation of the curled wake model was for a single wind turbine wake. In the case of a wind plant, the wakes are first computed individually and then superposed to obtain the flow field of the entire wind plant (Bay et al.2020). Most wake models are used in the same manner by first computing the wake of the individual turbines and using a superposition method afterward to obtain the flow over the entire domain (Annoni et al.2018). This new curled wake solver overcomes the use of a superposition method by solving the flow over the entire wind plant. This allows us to realize the benefits of the curled wake model in a much faster time frame with better scaling with domain size.

The curled wake solver presented in this work focuses on reducing computational cost and capturing wake steering effects. This is done by solving only the streamwise component of the linearized RANS equations and parametrizing the effects of the spanwise and wall-normal components using semianalytical solutions.

2 Formulation

We use the RANS equations to model the time-averaged flow field through a wind plant. The continuity equation is

(1) u x + v y + w z = 0 ,

and the RANS momentum equation for the streamwise direction is

(2) u u x + v u y + w u z = - 1 ρ p x - u u x - u v y - u w z + ν 2 u x 2 + 2 u y 2 + 2 u z 2 + 2 Ω y G z - w - 2 Ω z G y - v ,

where x, y, and z are the streamwise, spanwise, and wall-normal directions; u, v, and w are the velocity components in the respective directions (with  denoting time fluctuations and the overbar being time averaging); p is the time-averaged pressure; ρ is density; Ω is the Earth's rotational vector projected into an arbitrary, non-inertial, Earth-fixed frame of reference Ω=ω[cos (ϕ)sin (θ), cos (ϕ)cos (θ),sin (ϕ)], where ϕ is the latitude, θ is the angle measured between the domain x axis and the easting axis, and ω is the Earth's rotation rate; and G is the geostrophic wind velocity vector. This is a similar equation used in the original formulation of the curled wake model, but now we focus on a new approach to derive the equations and some generalizations used for a wind plant approach as opposed to a single wind turbine wake. We assume that the viscous effects are small (high Reynolds number limit) and are neglected in the rest of derivation. We also assume that the boundary layer is neutral and it satisfies the geostrophic balance in the free atmosphere (Allaerts and Meyers2015; van der Laan and Sørensen2017). Coriolis effects are included in the momentum balance given by Eq. (2) as a result of their influence on the wind speed and direction shear in the atmosphere (Wyngaard2010). Further, we included non-traditional Coriolis effects (Ωy) because of their potential impact in the atmospheric boundary layer in the presence of heterogeneous roughness elements such as wind turbines or terrain complexity (Howland et al.2020a).

2.1 Decomposing the velocity

The velocity is decomposed into a background flow (capital letters) and a wake deficit (Δ) by

(3) u = U + Δ u , v = V + Δ v , w = W + Δ w , p = P + Δ p .

The time-averaged fields are denoted using overbars:

(4) u = U + Δ u , v = V + Δ v , w = W + Δ w , p = P + Δ p .

The temporal fluctuations are denoted using a hash mark ():

(5) u = U + Δ u , v = V + Δ v , w = W + Δ w , p = P + Δ p .

2.1.1 Background flow

The background flow (U, V, W) is the velocity of the domain without including the wind turbines and their wakes. The background flow formulation can be obtained from an analytical formulation such as the log law or from a different time-averaged simulation. For example, you can specify uniform flow by U, V, W=U, 0, 0, or use simulation data from large-eddy simulation (LES) or experiments to define the background flow over complex terrain. For a consistent formulation of the model, the background flow should also satisfy the RANS equations.

2.1.2 Wake deficit solution

The time-averaged wake velocities are denoted by Δu, Δv, and Δw. We are interested in solving the streamwise component of the wake deficit, Δu, while the other wake velocity components, Δv and Δw, are parametrized using semianalytical models. The streamwise component of the RANS equations can be written in terms of the background flow and wake velocity as

(6) ( U + Δ u ) ( Δ u + U ) x + ( V + Δ v ) ( U + Δ u ) y + ( W + Δ w ) ( U + Δ u ) z = - 1 ρ ( P + Δ p ) x - ( U + Δ u ) ( U + Δ u ) x - ( U + Δ u ) ( V + Δ v ) y - ( U + Δ u ) ( W + Δ w ) z + 2 Ω y G z - ( W + Δ w ) - 2 Ω z G y - ( V + Δ v ) .

The background flow is defined to also satisfy the RANS equations as

(7) U U x + V U y + W U z = - 1 ρ P x - U U x - U V y - U W z + 2 Ω y G z - W - 2 Ω z G y - V .

Subtracting the background flow (Eq. 7) from the full flow (Eq. 6) leads to the equation of the curled wake model:

(8) ( U + Δ u ) Δ u x + ( V + Δ v ) Δ u y + ( W + Δ w ) Δ u z + Δ u U x + Δ v U y + Δ w U z = - 1 ρ Δ p x - ( 2 U Δ u + Δ u Δ u ) x - ( U Δ v + V Δ v + Δ u Δ v ) y - ( U Δ w + W Δ u + Δ u Δ w ) z + 2 Ω z Δ v - Ω y Δ w .

We now assume that the pressure gradient has a small effect (especially in the far wake), the Reynolds stresses are modeled using the turbulent-viscosity hypothesis (Pope2000), the second derivative of the wake deficit in the streamwise direction is neglected 2Δux2=0, the gradients of the mean flow are assumed to be small, and their influence on the convective terms is neglected:

(9) Δ u U x + Δ v U y + Δ w U z ( U + Δ u ) Δ u x + ( V + Δ v ) Δ u y + ( W + Δ w ) Δ u z ,

and the Coriolis terms in the velocity deficit equation are neglected (i.e., 2ΩzΔv(U+Δu)Δux). The Coriolis terms are included in the background flow RANS momentum balance (Eq. 7). These assumptions lead to the final form of the equation:

(10) Δ u x = - 1 U + Δ u ( V + Δ v ) Δ u y + ( W + Δ w ) Δ u z + ν eff 2 Δ u y 2 + 2 Δ u z 2 .

Equation (10) is the fundamental parabolic equation solved in the model presented. The streamwise velocity deficit, Δu, is the main unknown; all the other variables in the equation are either known a priori or parametrized at run time depending on Δu. We note that mass conservation was used in the derivation; however, it is not strictly enforced when solving Eq. (10). The equation is solved by marching in the downstream direction starting from an initial condition where the first wind turbine is (Sect. 3).

2.2 Turbulence model

The effect of turbulence in the RANS equations is described by the divergence of the Reynolds stress tensor. The streamwise component of the divergence of the Reynolds stress for the background flow solution (Eq. 7) is

(11) U U x + U V y + U W z .

The Reynolds stress term in Eq. (8) (for the wake deficit solution) is defined as

(12) ( 2 U Δ u + Δ u Δ u ) x + ( U Δ v + V Δ v + Δ u Δ v ) y + ( U Δ w + W Δ u + Δ u Δ w ) z .

The decomposition of the velocity field (background + wake, mean + fluctuation) leads to the introduction of additional stress-like terms in Eq. (8). These terms are correlations between the background flow solution and the wake deficit solution.

A mixing length model is used to represent the terms in Eq. (12). We propose using the simple model suggested in the original formulation of the curled wake model (Martínez-Tossas et al.2019) and scale the viscosity to take into account the effect from all of the extra terms in the Reynolds stresses from Eq. (12). This is the same approach suggested by Bay et al. (2019). The mixing length and eddy viscosity are defined as

(13) m = κ z ( 1 + κ z / λ ) ν eff = C m 2 d U d z ,

where m is the mixing length, νeff is the turbulent viscosity, κ is the von Kármán constant, z is the distance from the ground, and λ is the value of the mixing length in the free atmosphere. Blackadar (1962) proposed that λ0.00027G/fc, where G is the geostrophic wind speed magnitude and fc is the Coriolis parameter, resulting in a value of λ which is latitude dependent. Using typical values of G=10 m s−1 and mid-latitude ϕ=45, λ≈27 m. The precise value of λ, and more broadly m, will depend on independent parameters such as the atmospheric stability (Sun2011). In this study, we select λ=27 m and we suggest the investigation of more refined turbulence models for future work. The constant C is used to account for the additional turbulence introduced by the rotor and the wake. Tests have shown that for all of the cases tried in the manuscript, a value of C=4 has provided good agreement between the model and experiments/simulations. This value is consistent with what is suggested by Bay et al. (2019). Section A shows a sensitivity analysis of power in a wind plant simulation using the model to the constant C. We also expand the equations for the turbulence model in Sect. B. The mixing length and turbulent viscosity are difficult to approximate with constant values that depend only on height (z). A better approximation would allow turbulent viscosity to vary spatially, especially in the wake, where the local turbulence varies with the spanwise and streamwise coordinates.

The Reynolds stress model used in this study was selected because of its computational efficiency. Resolving the spatial variations in the eddy viscosity would require the solution of the full RANS momentum equations and additional transport equations for relevant parameters in the selected Reynolds stress model (van der Laan et al.2015a; Iungo et al.2018). Future work should investigate Reynolds stress models that are able to resolve the enhanced mixing and turbulence induced by the wind turbines while remaining computationally efficient for the hybrid RANS-analytical framework.

2.3 Wind turbine wakes initial condition

Wakes are initialized according to the wind speed at the rotor location in the plane closest to where the turbine is. As the solution marches downstream and new wind turbines are encountered, a new wake deficit is added to the plane (n):

(14) Δ u n = - 2 a U + Δ u n - 1 ,

where a=1-1-CTcos2α/2 is the induction from momentum theory, α is the yaw angle, CT is the thrust coefficient, and U+Δun-1 is the averaged velocity inside the disk in the plane upstream (n−1) of the rotor. The power and thrust coefficients are obtained from a lookup table based on the local velocity U+Δun-1. A Gaussian filter is used to smear the initial condition in the spanwise directions to avoid numerical instabilities described in Martínez-Tossas et al. (2019). The effects of wake curl, wake rotation, and the boundary layer are implemented using the analytical models also described in Martínez-Tossas et al. (2019). For completeness, we show the analytical formulas for the spanwise velocities from the curled wake. The effect of curl is added by modifying the spanwise velocity components according to an elliptic distribution of vorticity (Shapiro et al.2018; Martínez-Tossas et al.2019; Martínez-Tossas and Branlard2020). The spanwise velocities can be represented analytically by


where R is the turbine radius, y and z are the coordinates relative to the disk center, and Γ0=D2CTUsinαcos2α is the total circulation from yaw (Shapiro et al.2018; Martínez-Tossas et al.2019). The power and thrust coefficients are computed using the tabulated value at zero yaw angle as follows:

(17) P ( α ) = P ( α = 0 ) cos 2 ( α ) , T ( α ) = T ( α = 0 ) cos 2 ( α ) .

This relation has been used in previous work, but field experimental studies indicate that these functions are not necessarily powers of cosines and can be turbine-specific (Howland et al.2020b). The curled wake model presented here allows any function to be used to relate the power and thrust coefficient as a function of yaw angle, and future work will be focused on improving the functional relations between thrust, power, and yaw angle (e.g., model proposed by Howland et al.2020b). The solver computes the power and thrust from each turbine according to the local velocity at the rotor plane.

3 Numerical solution

Equation (10) is solved using numerical differentiation. Equation (18) shows the equation to be solved numerically with all of the terms labeled that are to be discretized:

(18) Δ u [ i + 1 , j , k ] = Δ u [ i , j , k ] - Δ x U + Δ u A ( V + Δ v ) Δ u y B + ( W + Δ w ) Δ u z C + ν eff 2 Δ u y 2 D + 2 Δ u z 2 E .

The discrete form of Eq. (18) is presented in Eq. (19).

(19) Δ u [ i + 1 , j , k ] = Δ u [ i , j , k ] - Δ x U + Δ u A V [ i , j , k ] + Δ v [ i , j , k ] Δ u [ i , j + 1 , k ] - Δ u [ i , j - 1 , k ] 2 Δ y B + W [ i , j , k ] + Δ w [ i , j , k ] Δ u [ i , j , k + 1 ] - Δ u [ i , j , k - 1 ] 2 Δ z C + ν eff Δ u [ i , j + 1 , k ] - 2 Δ u [ i , j , k ] + Δ u [ i , j - 1 , k ] Δ y 2 D + Δ u [ i , j , k + 1 ] - 2 Δ u [ i , j , k ] + Δ u [ i , j , k - 1 ] Δ z 2 E .

This numerical equation is discretized using a forward-in-time centered-in-space method with the stability criteria shown in Eq. (20) (Hoffman and Frankel2018; Martínez-Tossas et al.2019).

(20) Δ x 2 ν eff Δ u ( W + Δ w ) 2 , Δ y 2 ν eff Δ x / Δ u

We note that the model proposed is steady state and there is no time dependency. The spatial streamwise direction is treated as the “forward-in-time” part of the numerical method. The equations can be solved as a marching problem in the streamwise direction (index i) starting with an initial condition in a yz plane. The boundary conditions are set to zero wake deficit (Δu=0). Our tests have shown that the implementation has a converged and stable solution when using a grid resolution on the order of DΔy=10 in the spanwise directions (y and z) and DΔx=20 in the streamwise direction. A grid convergence study is shown in Sect. C. All the simulations and results presented were performed using uniform grid spacing.

Figure 1 shows a schematic of how the solution is computed. The main figure is a contour of streamwise velocity from a simulation with a random arrangement of turbines. The solution is marched downstream by solving Eq. (19) at each plane. Two planes are shown from the middle of the domain. The final solution includes a collection of planes for each streamwise location, which are combined to generate a full 3D solution.

Figure 1Schematic of the computational strategy used to solve Eq. (18). Dashed lines denote the location of a subset of planes, and big arrow shows the marching direction.


3.1 Computational cost

To better understand the low computational cost of the solver presented, we assess the number of floating point operations needed to obtain a solution to Eq. (16). We estimate the computational expense of the implementation by approximating the number of floating point operations (summation, subtraction, multiplication, division) in each term in Eq. (19). We assume that the total number of grid points in the computational domain is N. To solve Eq. (19), all the grid points in the domain must compute each of the terms in the equation. This leads to the following computational expense from each term: A=2N, B=4N, C=4N, D=5N, and E=5N, and, assuming one floating point operation between terms (4N), this leads to a total computational expense of ≈24N floating point operations. Assuming that we use a standard processor (1 GFLOPS), the computational time required for a simulation with N=1003 grid points based on this approximation would be 0.02 s. This can be considered an extremely fast solver for wind plant controls and layout optimization. In practice, the computational expense of the algorithm heavily depends on the implementation and software stack used. In our current implementation within the NumPy and Python frameworks (van der Walt et al.2011), the typical computational cost of a simulation is on the order of 0.1–10 s. This is 2 orders of magnitude faster than the standard curl model implementation in the FLOw Redirection and Induction in Steady State (FLORIS) framework. Figure 2 shows the time to solution of the algorithm as a function of the total number of grid points from the model presented compared to the standard FLORIS implementation with wake superposition (Bay et al.2019) compared to the linear scaling of the new solver. Also, the wind plant used for the scaling study is shown for reference. Significant speedup is expected in the presently proposed curled wake model formulation compared to the standard FLORIS implementation. The standard FLORIS implementation solves Eq. (16) for every turbine in the domain individually and then superposes the solutions. This superposition approach results in an increased computational cost, especially when more turbines are included, as well as wake superposition uncertainty. The resolutions used are finer than required for this wind plant, and the simulations lasting 0.5 s are converged and would be used for production runs. We note that this version of the model has not been optimized for performance, and future work will include code optimization and shared memory parallelization.

Figure 2Scaling of the computational algorithm (a) based on a representative wind plant composed of 36 turbines with wake steering (b).


3.2 Complex terrain capabilities

The current solver and formulations can also be used in applications with complex terrain. The complex terrain geometry can be included by specifying the boundary condition (Δu=0) along the terrain boundary. We test the model presented on a case with complex terrain over the Columbia River Gorge (Quon et al.2019). This test case is used to demonstrate the capabilities of the model in complex terrain conditions. The background flow solution is taken from a time-averaged LES (Quon et al.2019). The background flow is taken from LES, and the algorithm provides the solution for the wake deficits that would be present if turbines were there. Figure 3 shows a streamwise velocity contour for a plane in the streamwise direction. It is interesting to see how the wakes advect sideways following the background flow. Also, the combination of wakes leads to asymmetric deformation not typically observed in wakes over flat terrain. These results serve as a test case to show the applicability of the model in a case with complex terrain; further work is needed to assess the accuracy of the model under complex terrain conditions.

Figure 3Streamwise velocity contours showing a plane perpendicular to the predominant wind direction.


4 Results

We use the model proposed to compare with two different cases. The first comparison is done using supervisory control and data acquisition (SCADA) data from the Lillgrund wind plant. Second, we compare the model to a series of LES for an array of turbines with different yaw combinations.

4.1 Lillgrund wind plant

We use the model proposed to compute the flow field over the Lillgrund wind plant. Ten-minute average SCADA data are available for all turbines for different wind conditions. The SCADA data were organized by wind speed, turbulence intensity, and wind direction into bins with a width of 1 m s−1, 2 %, and 5. Three conditions from directions where the meteorological tower is not waked were chosen (185 with forty-one 10 min averages, 215 with ninety-three 10 min averages, and 255 with ninety 10 min averages). For each wind condition, we perform one simulation with the solver proposed. Figure 4 shows the layout of the Lillgrund wind plant with arrows denoting the directions for the cases that were studied. The background flow was set to a log-law streamwise velocity profile with a roughness height of z0=10-5 m.

Figure 4Layout of the Lillgrund wind plant with wind direction used in each simulation.


Figure 5 shows a comparison of power output between the SCADA data and the model proposed with a streamwise velocity contour at hub height. The data have been normalized according to the highest mean power in the experimental data. The bars in the SCADA data indicate the standard deviation of the power measurements. The mean absolute error in power is 8 %, 10 %, and 16 % for the cases with 185, 215, and 255, respectively. The agreement between the SCADA data and the model is good, with most results from the proposed solver lying within 1 standard deviation of the measurements. We can see different features of the flow, including the superposition of wakes. This allows for the solver to reach an equilibrium state in the deep array region. In this area, the power produced by the turbine flattens, providing a balance between the turbulent diffusion and the power extraction (Calaf et al.2010). Future work will focus on including wind direction uncertainty in the curled wake model (Gaumond et al.2014; van der Laan et al.2015b; Simley et al.2020).

Figure 5Comparison of turbine power versus SCADA data for the Lillgrund wind plant for cases at three different wind directions (185, 215, 255). Streamwise velocity contours at hub height are shown for all cases with the wind plant aligned with the flow direction.


Table 1List of LES cases performed for comparison study.

Download Print Version | Download XLSX

4.2 Wake steering

We now compare the model to results from LES of wakes in steering conditions. The simulations were performed using the Simulator fOr Wind Farm Applications (SOWFA) using an actuator disk model with rotation (Churchfield et al.2012). The turbine aerodynamics properties and control system are derived from the NREL 5 MW reference turbine (Jonkman et al.2009). The simulations are for cases with wind plants of four-by-three and three-by-three turbines with different offsets and yaw-angle combinations. The cases with four-by-three have spacing in the streamwise and spanwise directions of Sx=10 D and Sy=2.5 D respectively. The cases with a three-by-three array have spacing in the streamwise and spanwise directions of Sx=10 D and Sy=3 D respectively. The simulations use a precursor simulation from a neutral atmospheric boundary layer with roughness height of z0=0.15 m, wind direction of 270, and wind speed at hub height (90 m) of 8 m s−1. The simulations are time-averaged over 1600 s. Table 1 shows the main parameters for the simulations. The curled wake model uses the time-averaged LES inflow as the background flow.

Figure 6Total power output for wind plant LES with wake steering compared with the model proposed.


Figure 7Velocity at hub height normalized by average speed at hub height from the model proposed and from LES and power output for each turbine from the model proposed compared to results from LES. Simulations of a four-by-three turbine array. The bars in the LES power denote 1 standard deviation of the power. Turbine numbering is from bottom to top and left to right.


Figure 8Same as Fig. 8 but for a four-by-three array.


Figure 6 shows the total power for each case from the model proposed and from LES. There is good agreement in total power between the model and the simulations. The model proposed is able to capture the effects of yaw and general trends of power output from the different configurations. We note that the simulations still have some transient effects, and differences arise from these effects in the atmospheric boundary layer, including low-velocity streaks passing through the turbines (Munters et al.2016; Stevens et al.2018).

We select two representative cases and compare the power for all turbines and velocity at hub height. Figures 7 and 8 show power for all turbines and a velocity profile at hub height from LES and the model for cases 4 and 9. The turbine power for each turbine from the model in all cases is always within 1 standard deviation of the plots. The streaks from the precursor simulation create some of the differences in turbine power on the first row. The streaks are long structures that persist in the domain for very long periods of time (Munters et al.2016; Stevens et al.2018). To take the streaks into account, they are included as part of the background solution in the model. We notice that there are differences in the near wake between the model and the LES. These differences are present because the representation of the near wake is not well captured in the model. A better representation of the near wake and a more sophisticated turbulence model that can take into account the wake-added turbulence will be part of future work.

5 Conclusions

Fast wind power plant flow solvers are much needed for wind plant controls and layout optimization. In this work, we presented a simplified and fast solver for wind turbine wakes based on the curled wake model presented in Martínez-Tossas et al. (2019). The approach uses a hybrid RANS-analytical framework to obtain the wake velocity based on a parabolic solution for the streamwise component of the RANS equations. The computational expense of the model was shown to be on the order of seconds for a full wind plant with 36 turbines. The model was tested on three different cases: (1) SCADA data from the Lillgrund wind plant, (2) LES for flow over complex terrain, and (3) LES over flat terrain with different yaw-angle combinations. The models showed good agreement with the SCADA data from the Lillgrund wind plant. The model was also able to generate wake profiles for data in complex terrain, and future work will focus on comparing these profiles to data. Finally, the solver was able to reproduce the trends from LES with different yaw combinations. The model presented was shown to be an extremely fast solver (order of seconds) for wind turbine wakes with terrain features. This was achieved by simplifying the streamwise component of the RANS equation and making a series of assumptions. This model leverages approximations, especially with regard to the turbulence model, to improve computational speed. This trade-off provides a very computationally efficient solver at the expense of less robust turbulence modeling, compared to full three-dimensional RANS solver (van der Laan et al.2015a; Iungo et al.2018).

Some of the limitations from the different approximations of the model include a turbulence model mixing length that only depends on the vertical coordinate, a linearized solution of the vortices from curl that do not decay, a missing near-wake formulation, and no pressure term in the equations. These approximations were done to reduce the computational cost. Future work will focus on addressing the limitations and, more specifically, comparing the model with RANS, improving the turbulence model without compromising computational cost, improving the near wake, implementing a vortex decay model, using the solver for yaw-angle optimizations in a wind plant, and improving code performance to increase speed. This solver will soon be incorporated into the FLORIS framework and will be freely available.

Appendix A: The constant of the turbulence model

The turbulence model proposed uses a constant, C, to scale the turbulent viscosity. This constant is used to represent the wake-added turbulence. We performed a series of simulations with different values of C to tune the model constant. Figure A1 shows the SCADA power from the Lillgrund wind farm compared to results from the curled wake solver using a different value of the turbulent viscosity scaling, C. A value of C=4 provided best agreement between the curled wake model and Lillgrund SCADA data. We also note that this value agrees with previous observations from Bay et al. (2019).

Figure A1Comparison of turbine power versus SCADA data for the Lillgrund wind plant for cases at two different wind directions (185, 215). Different lines denote values of the constant C used to scale the viscous term in the solver.


Appendix B: Turbulence modeling

In this section, we show a different formulation for the development of the turbulence model. It is possible to invoke the eddy-viscosity hypothesis in the derivation for the base and wake deficit equations independently. When doing this, we express the Reynolds stresses as

(B1) - ( 2 U Δ u + Δ u Δ u ) x - ( U Δ v + V Δ v + Δ u Δ v ) y - ( U Δ w - W Δ u - Δ u Δ w ) z = ν f - ν b 2 U x 2 + 2 U y 2 + 2 U z 2 + ν f 2 Δ u x 2 + 2 Δ u y 2 + 2 Δ u z 2 ,

where νb is the turbulent viscosity of the base flow and νf is the turbulent viscosity of the full equation. We now choose to represent the Reynolds stresses as a function of the gradients of the wake deficit solution.

(B2) ν eff 2 Δ u x 2 + 2 Δ u y 2 + 2 Δ u z 2 = ν f - ν b 2 U x 2 + 2 U y 2 + 2 U z 2 + ν f 2 Δ u x 2 + 2 Δ u y 2 + 2 Δ u z 2 ,

where νeff is an effective viscosity if we only use the gradients of the wake deficit. Rearranging the equation, the effective turbulent viscosity can be defined as

(B3) ν eff = ν f - ν b 2 U x 2 + 2 U y 2 + 2 U z 2 2 Δ u x 2 + 2 Δ u y 2 + 2 Δ u z 2 + ν f .

We can define the effective viscosity as

(B4) ν eff = C ν b ,

and subtracting from Eq. (B3) leads to

(B5) C = ν f / ν b - 1 2 U x 2 + 2 U y 2 + 2 U z 2 2 Δ u x 2 + 2 Δ u y 2 + 2 Δ u z 2 + ν f / ν b .

The value of C is a function of space. However, in this work, we have chosen a constant value for C that minimizes the error between the observations and the model results.

Appendix C: Grid refinement study of the curled wake model

Here, we show a convergence study for the curled wake model based on one of the simulations for the Lillgrund wind farm in Sect. 4.1. We evaluate grid convergence using power output from turbines. The power is computed by taking the velocity average in the rotor area for each turbine and using a lookup table. Figure C1 shows the power for all turbines, and each line represents a different resolution in the spanwise directions. The results converge when using D/Δy=9. At this point, the average error percentage in power for the finest resolution is below 3 %. We also refined the streamwise direction for all cases studied, but we noticed that the error in power from refining the grid in the streamwise direction is always less than 1 % as long as the grid meets the stability criteria presented in Sect. 3.

Figure C1Comparison of turbine power for all turbines using different number of grid points across the turbine diameter for the Lillgrund wind plant for cases at wind direction of 185.


Code availability

Code will be available upon request by contacting the correspondence author.

Author contributions

LAMT led the model development and wrote the article. All authors provided input to this paper.

Competing interests

The authors declare that they have no conflict of interest.


The views expressed in the article do not necessarily represent the views of the DOE or the US Government. The US Government retains and the publisher, by accepting the article for publication, acknowledges that the US Government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this work, or allow others to do so, for US Government purposes.


The authors would like to thank Paula Doubrawa for help in allocating funding for this work and Sheri Anstedt from NREL for their help in editing the manuscript. The authors would also like to acknowledge the comments from Paul van der Laan and the anonymous reviewer. Data were furnished to the authors under an agreement between the National Renewable Energy Laboratory, Siemens Gamesa Renewable Energy A/S, and Vattenfall. Data and results used herein do not reflect findings by Siemens Gamesa Renewable Energy A/S and Vattenfall.

Financial support

A portion of the research was performed using computational resources sponsored by the US Department of Energy's Office of Energy Efficiency and Renewable Energy and located at the National Renewable Energy Laboratory. This work was authored by the National Renewable Energy Laboratory, operated by Alliance for Sustainable Energy, LLC, for the US Department of Energy (DOE) under contract no. DE-AC36-08GO28308. Funding provided by the US Department of Energy Office of Energy Efficiency and Renewable Energy Wind Energy Technologies Office.

Review statement

This paper was edited by Katherine Dykes and reviewed by Paul van der Laan and one anonymous referee.


Adaramola, M. and Krogstad, P.-A.: Experimental investigation of wake effects on wind turbine performance, Renew. Energy, 36, 2078–2086, 2011. a

Ainslie, J.: Calculating the flowfield in the wake of wind turbines, J. Wind Eng. Indust. Aerodynam., 27, 213–224,, 1988.  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,, 2015. a

Annoni, J., Fleming, P., Scholbrock, A., Roadman, J., Dana, S., Adcock, C., Porte-Agel, F., Raach, S., Haizmann, F., and Schlipf, D.: Analysis of control-oriented wake modeling tools using lidar field results, Wind Energ. Sci., 3, 819–831,, 2018. a

Bartl, J., Mühle, F., Schottler, J., Saetran, L., Peinke, J., Adaramola, M., and Hölling, M.: Wind tunnel experiments on wind turbine wakes in yaw: effects of inflow turbulence and shear, Wind Ener. Sci., 3, 329–343,, 2018. a

Bastankhah, M. and Porté-Agel, F.: Experimental and theoretical study of wind turbine wakes in yawed conditions, J. Fluid Mech., 806, 506–541,, 2016. a, b, c, d

Bay, C. J., Annoni, J., Martínez-Tossas, L. A., Pao, L. Y., and Johnson, K. E.: Flow Control Leveraging Downwind Rotors for Improved Wind Power Plant Operation, in: 2019 American Control Conference (ACC), 10–12 July 2019, Philadelphia, PA, USA, 2843–2848, 2019. a, b, c, d, e

Bay, C. J., King, J., Martínez-Tossas, L. A., Mudafort, R., Hulsman, P., Kühn, M., and Fleming, P.: Towards flow control: an assessment of the curled wake model in the FLORIS framework, in: The Science of Making Torque from Wind (TORQUE 2020), 28 September–2 October 2020, the Netherlands, 2020. a, b, c

Bradstock, P. and Schlez, W.: Theory and verification of a new 3D RANS wake model, Wind Energ. Sci., 5, 1425–1434,, 2020. a, b

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

Churchfield, M., Lee, S., Moriarty, P., Martinez, L., Leonardi, S., Vijayakumar, G., and Brasseur, J.: A Large-Eddy Simulation of Wind-Plant Aerodynamics, in: 50th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, 9–12 January 2012,Nashville, Tennessee,, 2012. a

Fleming, P., Annoni, J., Churchfield, M., Martinez-Tossas, L. A., Gruchalla, K., Lawson, M., and Moriarty, P.: A simulation study demonstrating the importance of large-scale trailing vortices in wake steering, Wind Energ. Sci., 3, 243–255,, 2018a. a

Fleming, P., Annoni, J., Martínez-Tossas, L. A., Raach, S., Gruchalla, K., Scholbrock, A., Churchfield, M., and Roadman, J.: Investigation into the shape of a wake of a yawed full-scale turbine, J. Phys.: Conf. Ser., 1037, 032010,, 2018b. a

Fleming, P., King, J., Dykes, K., Simley, E., Roadman, J., Scholbrock, A., Murphy, P., Lundquist, J. K., Moriarty, P., Fleming, K., van Dam, J., Bay, C., Mudafort, R., Lopez, H., Skopek, J., Scott, M., Ryan, B., Guernsey, C., and Brake, D.: Initial results from a field campaign of wake steering applied at a commercial wind farm – Part 1, Wind Energ. Sci., 4, 273–285,, 2019. a

Gaumond, M., Réthoré, P.-E., Ott, S., Peña, A., Bechmann, A., and Hansen, K. S.: Evaluation of the wind direction uncertainty and its impact on wake modeling at the Horns Rev offshore wind farm, Wind Energy, 17, 1169–1178,, 2014. a

Gebraad, P. M. O., Teeuwisse, F. W., van Wingerden, J. W., Fleming, P. A., Ruben, S. D., Marden, J. R., and Pao, L. Y.: Wind plant power optimization through yaw control using a parametric model for wake effects?a CFD simulation study, Wind Energy, 19, 95–114, 2016. a

Gunn, K., Stock-Williams, C., Burke, M., Willden, R., Vogel, C., Hunter, W., Stallard, T., Robinson, N., and Schmidt, S. R.: Limitations to the validity of single wake superposition in wind farm yield assessment, J. Phys.: Conf. Ser., 749, 012003,, 2016. a

Hamilton, N., Bay, C. J., Fleming, P., King, J., and Martínez-Tossas, L. A.: Comparison of modular analytical wake models to the Lillgrund wind plant, J. Renew. Sustain. Energ., 12, 053311,, 2020. a

Hoffman, J. D. and Frankel, S.: Numerical Methods for Engineers and Scientists, 2nd Edn., CRC Press, Boca Raton, 840 pp.,, 2018. a

Howland, M. F., Bossuyt, J., Martínez-Tossas, L. A., Meyers, J., and Meneveau, C.: Wake structure in actuator disk models of wind turbines in yaw under uniform inflow conditions, J. Renew. Sustain. Energ., 8, 043301,, 2016. a, b, c

Howland, M. F., Lele, S. K., and Dabiri, J. O.: Wind farm power optimization through wake steering, P. Natl. Acad. Sci. USA, 116, 14495–14500,, 2019. a

Howland, M. F., Ghate, A. S., and Lele, S. K.: Influence of the geostrophic wind direction on the atmospheric boundary layer flow, J. Fluid Mech., 883, A39,, 2020a. a

Howland, M. F., González, C. M., Martínez, J. J. P., Quesada, J. B., Larrañaga, F. P., Yadav, N. K., Chawla, J. S., and Dabiri, J. O.: Influence of atmospheric conditions on the power production of utility-scale wind turbines in yaw misalignment, J. Renew. Sustain. Energ., 12, 063307,, 2020b. a, b

Hulsman, P., Martinez-Tossas, L. A., Hamilton, N., and Kuhn, M.: Modelling and assessing the near-wake representation and turbulence behaviour of control-oriented wake models, in: The Science of Making Torque from Wind (TORQUE 2020), 28 September–2 October 2020, the Netherlands, 2020. a

Iungo, G. V., Santhanagopalan, V., Ciri, U., Viola, F., Zhan, L., Rotea, M. A., and Leonardi, S.: Parabolic RANS solver for low-computational-cost simulations of wind turbine wakes, Wind Energy, 21, 184–197,, 2018. a, b, c

Jonkman, J., Butterfield, S., Musial, W., and Scott, G.: Definition of a 5-MW reference wind turbine for offshore system development, Tech. rep., NREL – National Renewable Energy Lab., Golden, CO, USA, 2009. a

King, J., Fleming, P., King, R., Martínez-Tossas, L. A., Bay, C. J., Mudafort, R., and Simley, E.: Controls-Oriented Model for Secondary Effects of Wake Steering, Wind Energ. Sci. Discuss. [preprint],, in review, 2020. a, b

Martínez-Tossas, L. A. and Branlard, E.: The curled wake model: equivalence of shed vorticity models, J. Phys.: Conf. Ser., 1452, 012069,, 2020. a

Martínez-Tossas, L. A., Annoni, J., Fleming, P. A., and Churchfield, M. J.: The aerodynamics of the curled wake: a simplified model in view of flow control, Wind Energy Sci., 4, 127–138,, 2019. a, b, c, d, e, f, g, h, i, j, k, l

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

Munters, W., Meneveau, C., and Meyers, J.: Shifted periodic boundary conditions for simulations of wall-bounded turbulent flows, Phys. Fluids, 28, 025112,, 2016. a, b

Park, J., Kwon, S., and Law, K. H.: Wind farm power maximization based on a cooperative static game approach, Proc. SPIE, 8688, 86880R, 2013. a

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

Quon, E. W., Doubrawa, P., Annoni, J., Hamilton, N., and Churchfield, M. J.: Validation of Wind Power Plant Modeling Approaches in Complex Terrain, in: AIAA Scitech 2019 Forum, 7–11 January 2019, San Diego, California,, 2019. a, b

Schottler, J., Bartl, J., Mühle, F., Sætran, L., Peinke, J., and Hölling, M.: Wind tunnel experiments on wind turbine wakes in yaw: redefining the wake width, Wind Energ. Sci., 3, 257–273,, 2018. a

Shapiro, C. R., Gayme, D. F., and Meneveau, C.: Modelling yawed wind turbine wakes: a lifting line approach, J. Fluid Mech., 841, R1,, 2018. a, b, c

Siemens Gamesa: Siemens Gamesa now able to actively dictate wind flow at offshore wind locations, Press release, available at: (last access: 13 April 2021), 2019. a

Simley, E., Fleming, P., and King, J.: Design and analysis of a wake steering controller with wind direction variability, Wind Energ. Sci., 5, 451–468,, 2020. a

Stevens, R. J., Martínez-Tossas, L. A., and Meneveau, C.: Comparison of wind farm large eddy simulations using actuator disk and actuator line models with wind tunnel experiments, Renew. Energy, 116, 470–478, 2018. a, b

Sun, J.: Vertical variations of mixing lengths under neutral and stable conditions during CASES-99, J. Appl. Meteorol. Clim., 50, 2030–2041, 2011. a

van der Laan, M. P. and Sørensen, N. N.: Why the Coriolis force turns a wind farm wake clockwise in the Northern Hemisphere, Wind Energ. Sci., 2, 285–294,, 2017. a

van der Laan, M. P., Sorensen, N. N., Rethore, P.-E., Mann, J., Kelly, M. C., Troldborg, N., Schepers, J. G., and Machefaux, E.: An improved k-ϵ model applied to a wind turbine wake in atmospheric turbulence, Wind Energy, 18, 889–907,, 2015a. a, b

van der Laan, M. P., Sørensen, N. N., Réthoré, P.-E., Mann, J., Kelly, M. C., Troldborg, N., Hansen, K. S., and Murcia, J. P.: The k-ϵ-fP model applied to wind farms, Wind Energy, 18, 2065–2084,, 2015b. a

van der Walt, S., Colbert, S., and Varoquaux, G.: The NumPy Array: A Structure for Efficient Numerical Computation, Comput. Sci. Eng., 13, 22–30,, 2011. a

Vollmer, L., Steinfeld, G., Heinemann, D., and Kühn, M.: Estimating the wake deflection downstream of a wind turbine in different atmospheric stabilities: an LES study, Wind Energ. Sci., 1, 129–141,, 2016.  a

Wyngaard, J. C.: Turbulence in the Atmosphere, Cambridge University Press, Cambridge, 2010. a

Zong, H. and Porté-Agel, F.: A momentum-conserving wake superposition method for wind farm power prediction, J. Fluid Mech., 889, A8,, 2020. a

Publications Copernicus
Short summary
In this paper a three-dimensional steady-state solver for flow through a wind farm is developed and validated. The computational cost of the solver is on the order of seconds for large wind farms. The model is validated using high-fidelity simulations and SCADA.