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

This work focuses on minimizing the computational cost of steady-state wind power plant flow simulations that take into account wake steering physics. We present a simple wake solver with a computational cost on the order of seconds for large wind plants. The solver uses a simplified form of the Reynolds-averaged Navier-Stokes equations to obtain a parabolic equation for the wake deficit of a wind plant. We compare results from the model to supervisory control and data acquisition (SCADA) from the Lillgrund wind plant; good agreement is obtained. Results for the solver in complex terrain are also shown. Finally, the 5 solver is demonstrated for a case with wake steering showing good agreement with power reported by large-eddy simulations. This new solver minimizes the time–and therefore the related cost–it takes to conduct a steady-state wind plant flow simulation to about a second on a personal laptop. This solver can be used for different applications including wake steering for wind power plants and layout optimization, and it will soon be available within the FLOw Redirection and Induction in Steady State (FLORIS) framework. 10


Introduction
In this work, we present an improved formulation of the curled wake model  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é-Agel, 2020). This inconsistency motivates the use of a more robust solver 15 in the context of the curled wake model ) 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.
Parabolic solvers for RANS equations are a promising tool for fast wind farm flow solvers. Other researchers have developed parabolic solvers for wind plant applications. Iungo et al. (2018) developed a parabolic RANS solver focused on improving the 20 mixing-length model and used assumptions about axisymmetry in the wakes. Bradstock and Schlez (2019a, b) 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. WakeBlaster uses a special method to solve the spanwise velocity components that does not include effects caused by yaw. The curled wake solver presented in this work focuses on minimizing 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.
Wake steering is a promising wind plant control strategy used to maximize the power output of a wind plant (Adaramola and Krogstad, 2011;Park et al., 2013;Gebraad et al., 2016;Howland et al., 2019;Fleming et al., 2019;Siemens Gamesa, 2019). In wake steering, upstream turbines are yawed, deflecting the wakes such that downstream turbines are able to produce 30 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 controls-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é-Agel, 2016;Martínez-Tossas et al., 2019). This shape has been observed in computational fluid dynamics simulations 35 and in small-and large-scale experiments (Medici and Alfredsson, 2006;Howland et al., 2016;Bastankhah and Porté-Agel, 2016;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é-Agel, 2016;Shapiro et al., 2018;Martínez-Tossas et al., 2019). The curled wake is a unique phenomenon 40 in wind turbine wakes because it disrupts the asymmetry of the wake. The curled wake cannot be characterized by a symmetric profile such as a Gaussian distribution and requires a different modeling approach. 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).

45
The curled wake model uses a simplified version of the RANS equations to predict the wake of a wind turbine in yaw . 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., 2020;Hulsman et al., 2020). Also, the new Gauss-Curl Hybrid model has shown to provide a good compromise between an analytical model (Bastankhah and Porté-Agel, 2016) and some 50 of the physics from the curled wake model . 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, then superposed to obtain the flow field of the entire wind plant . 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 55 wind plant. This allows us to realize the benefits of the curled wake model in a much faster time frame with better scaling.
We use the RANS equations to model the time-averaged flow field through a wind plant. The RANS equation for the streamwise direction is: 60 where x is the streamwise direction; y is the spanwise direction; z is the wall-normal direction; u, v and w are the velocity components in the respective directions (with denoting time fluctuations and the overbar is time averaging); p is the timeaveraged pressure; and ρ is density. This is the same 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.

Decomposing the velocity
The velocity is decomposed into a background flow (capital letters) and a wake deficit (∆) by: The time-averaged fields are denoted using overbars:

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

Background flow
The background flow (U , V , W ) is the velocity of the domain where the wind turbines would be without including the wind turbines and their wakes. The background flow formulation can be obtained from an analytical formulation such as the log-75 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 LES or experiments to define the background flow over complex terrain.

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.

80
The streamwise component of the RANS equations can be written in terms of the background flow and wake velocity as: The background flow is defined to also satisfy the RANS equations as Subtracting the background flow (Equation 6) from the full flow (Equation 5) leads to the equation of the curled wake model: We now assume that the pressure gradient has a small effect (especially in the far wake). The Reynolds stresses are modeled as a viscous term using a mixing length model and are dominated by spanwise gradients (Pope, 2000). The gradients of the mean 90 flow are assumed to be small, and their influence on the convective terms is neglected: This leads to the final form of the equation: The equation is solved by marching in the downstream direction starting from an initial condition where the first wind turbine is (section 3).

Turbulence model
The effect of turbulence in the RANS equations is described by the divergence of the Reynolds stress tensor. The streamwise 100 component of the divergence of the Reynolds stress for the background flow solution (Equation 6) is: The Reynolds stress term in Equation 7 (for the wake deficit solution) is defined as: The decomposition of the velocity field (background + wake, mean + fluctuation) leads to the introduction of additional stress-105 like terms in Equation 7. These terms are correlations between the background flow solution and the wake deficit solution.
These terms need to be taken into account when solving Equation 7. A mixing length model is used to represent the terms in Equation 11. We propose using the simple model suggested in the original formulation of the curled wake model  and scale the viscosity to take into account the (2019). The mixing length and turbulent viscosity are defined as: 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, λ = 15m is the value of the mixing length in the free atmosphere. and C is a scaling constant (Blackadar, 1962;Sun, 2011).

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 with a diameter including the expansion of the wake: where a = (1 − √ 1 − C T )/2 is the induction from momentum theory, C T is the thrust coefficient, and < U + ∆u > is the averaged velocity inside the disk. A Gaussian filter is used to smear the initial condition in the spanwise directions to avoid

135
where R is the turbine radius, y and z are the coordinates relative to the disk center, and Γ 0 = D 2 C T U ∞ sin α cos 2 α is the total circulation from yaw (Shapiro et al., 2018;Martínez-Tossas et al., 2019).
It is now possible to write Equation 1 using numerical differentiation. Equation 16 shows the equation to be solved numerically with all of the terms labeled that are to be discretized: The terms in Equation 16 can be defined discretely as: This numerical equation is discretized using a forward-in-time, centered-in-space formulation with the stability criteria shown 145 in Equation 18 (Hoffman and Frankel, 2018;Martínez-Tossas et al., 2019). We note that the model proposed is steady stated 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): ∆x ≤ 2ν eff ∆u (W + ∆w) 2 , ∆y ≥ 2ν eff ∆x/∆u.

150
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-20 in the spanwise directions (y and z) and D ∆x 30-40 in the streamwise direction.

Computational cost
To better understand the low computational cost of the solver presented, we asses the number of floating point operations needed to obtain a solution to 15. We estimate the computational expense of the implementation by approximating the number of floating point operations (summation, subtraction, multiplication, division) in each term in Equation 16. We assume that   Figure 2 shows the time to solution of the algorithm as a function of total 170 number of grid points from the model presented compared to the standard FLORIS implementation with wake superposition  compared to the linear scaling of the new solver. Also, the wind plant used for the scaling study is shown for reference. The resolutions used are finer than required for this wind plant, and the simulations lasting 0.5 seconds are converged and would be used for production runs.

175
We use the model proposed to compare with three different cases. The first comparison is done using SCADA from the Lillgrund wind plant. Second, we showcase the use of the solver in complex terrain. Finally, we compare the model to a series of LES for an array of turbines with different yaw combinations.

Lillgrund Wind Plant
We use the model proposed to compute the flow field over the Lillgrund Wind Plant. Ten-minute average SCADA is available 180 for all turbines for different wind conditions. Three conditions from directions where the meteorological tower is not waked were chosen (185 o with 41 10-minute averages, 215 o with 93 10-minute averages, and 255 o with 90 10-minute averages). For each wind condition, we perform one simulation with the solver proposed. Figure 3 shows the layout of the Lillgrund Wind  according to the local velocity. 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; it is a balance between the turbulent diffusion and the power extraction.

Complex terrain: Columbia River Gorge
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 195 from a time-averaged LES (Quon et al., 2019). Figure 5 shows a volume rendering of streamwise velocity from a simulation using the proposed model. We can see the three-dimensionality of the solution and how the wakes conform to the terrain. 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 6 shows streamwise velocity contours for planes in all directions. It is interesting to see how the wakes advect sideways following the background flow. Also, the combination of wakes leads to asymmetric deformation not   10 https://doi.org/10.5194/wes-2020-86 Preprint. Discussion started: 9 July 2020 c Author(s) 2020. CC BY 4.0 License.

Wake steering
We now compare the model to results from LES of wakes in steering conditions. The simulations were performed using the Simulator fOr Applications (SOWFA) using an actuator disk model with rotation (Churchfield et al., 2012).   We select two cases and show the power for all turbines and a plane at hub height with wake profiles from the model proposed. Figure 8 shows power for all turbines and a velocity profile at hub height for cases 2 and 6. There is good agreement 215 between the model and the LES. The unsteadiness in the atmospheric boundary layer inflow creates some of the differences in turbine power. These differences are expected to diminish with longer time averaging. Also, we observed consistent differences in the power output of turbines further donwstream. These differences are caused by the lack of a more sophisticated turbulence model that can take into account the wake-added turbulence. We are currently working on improving the turbulence model and incorporating some of the findings from .

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).
This solver is based on a parabolic equation for the streamwise component of the Reynolds-averaged Navier-Stokes equation.
The computational expense was shown to be on the order of seconds for a full wind plant with 36 turbines. The model was tested 225 on three different cases: 1) SCADA 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 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.

230
This was achieved by simplifying the streamwise component of the RANS equation and making a series of assumptions. This model uses many simplifications, especially with regard to the turbulence model,to improve computational speed. This tradeoff provides a very computationally efficient solver at the expense of less robust turbulence modeling. This solver will soon be incorporated into the FLORIS framework and will be freely available. Future work will consist of implementing a vortex decay model and using the solver for yaw-angle optimizations in a wind plant.

235
Code availability. The code will soon be available withing the FLORIS framework.
Author contributions. LAMT led the model development and writing the article. All authors provided input to this manuscript.
Competing interests. No competing interest.

Acknowledgments
The authors would like to thank Paula Doubrawa for help in allocating funding for this work and Sheri Anstedt from NREL for The U.S. Government retains and the publisher, by accepting the article for publication, acknowledges that the U.S. Government retains a nonexclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this work, 250 or allow others to do so, for U.S. Government purposes.