The Aerodynamics of the Curled Wake: A Simpliﬁed Model in View of Flow Control

. When a wind turbine is yawed, the shape of the wake changes and a curled wake proﬁle is generated. The curled wake has drawn a lot of interest because of its aerodynamic complexity and applicability to wind farm controls. The main mechanism for the creation of the curled wake has been identiﬁed in the literature as a collection of vortices that are shed from the rotor plane when the turbine is yawed. This work extends that idea by using aerodynamic concepts to develop a control-oriented model for the curled wake based on approximations to the Navier-Stokes equations. The model is tested and 5 compared to large-eddy simulations using actuator disk and line models. The model is able to capture the curling mechanism for a turbine under uniform inﬂow and in the case of a neutral atmospheric boundary layer. The model is then tested inside the FLOw Redirection and Induction in Steady State framework and provides excellent agreement with power predictions for cases with two and three turbines in a row. not only its own wake, but the wake of other downstream turbines as well. Also, we note that the shed vortices allow for spanwise velocity components, which are vital when considering wake 25 redirection and wind farm controls. The vortices generated are not limited to only yawing, as they can also be used for tilt and combinations of tilt and yaw. This work sets a foundation for a simpliﬁed wake steering model to be used in a more general wind farm control-oriented framework. Acknowledgements. The for suggestions on turbulent modeling in the atmospheric boundary layer and Moriarty for providing wind turbine aerodynamics insights. Simulations for the NREL code SOWFA were performed using the National Renewable Energy Laboratory’s Peregrine high-performance computing system. Numerical implementation and plots were done using python with numpy, maptlotlib, and scipy libraries.


Introduction
A curled wake is a phenomenon observed in the wake of a wind turbine when the turbine is yawed relative to the free-stream velocity. When a wind turbine is yawed, the wake is not only deflected in a direction opposite to the yaw angle, but its shape changes. The mechanism behind this effect has drawn attention, not only from fluid dynamicists because of the interesting physics phenomena happening in the wake, but also from the controls community who intends to use it to control wind farm 15 flows (Fleming et al., 2017).
It has been shown that wake steering (redirection of the wake through yaw misalignment) can lead to an increase in power production of wind turbine arrays (Adaramola and Krogstad, 2011;Park et al., 2013;Gebraad et al., 2016). Previous studies have used large-eddy simulations (LES) and analytical tools to show the effects of yawing in the redirection of the wake (Jiménez et al., 2010;Bastankhah and Porté-Agel, 2016;Shapiro et al., 2018). Wind tunnel experiments have also been used to 20 study the wake of a wind turbine in yaw (Medici and Alfredsson, 2006;Bartl et al., 2018). The curled wake mechanism, in the context of wind turbine wakes, was first identified by Howland et al. (2016) during a porous disk experiment, and in LES using an actuator disk model (ADM) and actuator line model (ALM). This mechanism was described by a pair of counter-rotating vortices that are shed from the top and bottom of the rotor when the rotor is yawed. Further, these vortices move the wake to the side and create a curled wake shape. This mechanism was later confirmed and elaborated by the work of Bastankhah and Porté-Agel (2016) by performing experiments using particle image velocimetry of a scaled wind turbine and by doing a theoretical analysis using potential flow theory. Shapiro et al. show that the vorticity distribution shed from the rotor because of yaw has an elliptic shape as opposed to only a pair of counter-rotating vortices (Shapiro et al., 2018). The curled wake has 5 also been observed in LES with different atmospheric stabilities (Vollmer et al., 2016). Berdowski et al. (2018) were able to reproduce curled wake profiles by using a vortex method. Fleming et al. (2017) show that the generated vortices affect the performance of wake steering and motivate the development of engineering models (like the one in this paper), which include wake curling physics. Controllers based on such models would pursue different, and likely more effective, wind farm control strategies.

10
In this work, we describe the aerodynamics of the curled wake, and propose a new, simple and computationally efficient model for wake deficit, based on a linearized version of the Navier-Stokes equations with approximations. The model is tested and compared to LES using actuator disk and line models.

Aerodynamics of the Curled Wake: A Control-Oriented Model
Here, we develop a simplified model of the wake deficit considering the aerodynamics of the curled wake. We start by writing 15 the Reynolds-averaged Navier-Stokes streamwise momentum equation for an incompressible flow: where u is the streamwise velocity, v is the spanwise velocity, and w is the wall-normal velocity, p is the pressure, ρ is the fluid density, and ν eff is the effective (turbulent and molecular) viscosity. The flow field is now decomposed into a base solution and a perturbation about this solution. The base solution includes what we consider to be the main effects that convect the wake, 20 including, but not limited to: 1. The streamwise velocity profile 2. The rotational velocity from the shed vortices caused by yawing 3. The rotational velocity due to the blade rotation.
The velocity components can be expanded as: where U , V , and W are the streamwise, spanwise, and wall-normal velocity components from the base solution, and u , v , and w are the perturbation velocities. We use a reference frame where x is the streamwise direction, y is the spanwise direction, and z is the wall-normal direction. Assuming that the inflow is fully developed, the base solution is only a function of the spanwise and wall normal directions.
For this initial test of the model, the pressure gradient is assumed to be zero, which is a valid approximation in the far wake but not necessarily near the rotor. Also, we assume that the streamwise velocity, U , is not a function of the streamwise (x) and spanwise (y) coordinates. We also assume that the viscous force from the boundary layer is balanced by its pressure gradient.

5
With these simplifications, we are left with the equation: Equation 4 describes the downstream evolution of the wake deficit, u . We note that the only velocity component being solved is the streamwise component, u . The base solution for the flow (U, V, W ) includes the spanwise effects from different features, such as rotation and curl. Because these effects are the main drivers for wake deformation, the spanwise perturbations, v and 10 w , are assumed to be zero. This assumption reduces the model to a single partial differential equation. Equation 4 will be solved numerically in the next section. The model does not use the continuity equation so mass conservation is not strictly enforced for the perturbation velocities.

Curled Wake
The curled wake effect is added to the model by adding a distribution of counter-rotating vortices to the base flow solution.
15 Figure 1 shows a schematic of the rotor and a collection of vortices being shed at the rotor plane. The superposition of all the vortices leads to a spanwise velocity distribution that creates the curled wake shape. Each vortex is described as a Lamb-Oseen vortex with a tangential velocity distribution given by: where u t is the tangential component of the velocity, r is the distance from the vortex core, Γ is the circulation strength, and σ 20 determines the width of the vortex core. The circulation related to the strength of each vortex is still unknown. We assume that the problem is symmetric and that all the circulation leaves the disk through the shed vortices. In the current implementation, the vortices are assumed not to decay as they are convected downstream; however, it would be possible to incorporate this effect into the model. It is important to note that the vorticity shed has an elliptic distribution according to the shape of the rotor (Shapiro et al., 2018). The circulation is then: Γ = ρ D U ∞ F L , where ρ is the fluid density, D is the rotor diameter, U ∞ is the 25 inflow velocity at hub height, and F L is the force perpendicular to the inflow wind. The force component perpendicular to the inflow wind generates circulation. If we use the definition of the thrust coefficient, it is now possible to redefine the strength of the circulation as a function of thrust coefficient and yaw angle (Shapiro et al., 2018) as: where C T is the thrust coefficient for the given inflow velocity, and γ is the yaw angle. The discrete elliptic distribution of shed vortices added to match the shape of the disk is:

5
where i is the index denoting each of the vortices distributed on a line between the top and bottom of the rotor diameter, N is the total number of vortices, and the coordinates y i and z i are centered at the location of the shed vortex. The size of the vortex core is set to σ = D/5. The choice of σ intends to represent a realistic vortex core size for an elliptic distribution. This value represents σ/D ∼ 0.2, which is on the order of the optimal size for flow over an airfoil (Martínez-Tossas et al., 2017). It is possible to use other values, however, in the simulations presented, this value is a good compromise between a physical width 10 and numerical stability. The strength of each vortex is associated with the elliptical distribution by: where Γ 0 = 4/π Γ is used to ensure that the total amount of circulation Γ is conserved. This means that each vortex has a unique amount of circulation. When the distribution of circulation in Equation 9 is integrated, the total circulation is obtained. The results use N =200 discrete vortices, which has shown that this has little effect on the solutions, and values around N = 20 provide the same results. Noted that, when using just the two tip vortices in the top and bottom of the rotor, the shape of the curled wake does not match the simulation results. For this reason, an elliptic distribution of shed vorticity, which matches the shape of the rotor, should be used.

5
It is important to include wake rotation in the model, because the rotation will move the wake in a preferred direction. Wake rotation is taken into account by adding a tangential velocity distribution that is caused by the rotation inside the rotor area.
The tangential induction factor is defined as: where a is the induction factor based on the thrust coefficient from standard actuator disk theory, R is the rotor radius, r is the 10 radial distance from the center of the rotor, and λ is the tip speed ratio (Burton et al., 2002). From this equation, the tangential component of the velocity u t is a singular vortex with 1/r behavior: A Lamb-Oseen vortex is used to desingularize the behavior near the center of the rotor. The circulation strength for the wake rotation vortex based on Equation 10 is now: We assume that the wake rotation vortex does not decay or deform as it moves downtstream. This is not necessarily true, as turbulence mixing will decrease the wake rotation. However, the present model does not diffuse the spanwise velocity components and some of the errors in the model can be attributed to this. The current implementation uses a Lamb-Oseen vortex with a core size σ = D/5 to eliminate numerical instabilities caused by high velocities near the vortex core, but other 20 values could be used.

Atmospheric Boundary Layer
The atmospheric boundary layer can be specified as part of the background flow. A profile including streamwise and spanwise velocity components can be specified. The streamwise profile is described by using a power law:

25
where U h is the velocity at hub height, z h is the hub height, and α is the shear exponent. Also, it is possible to add a profile for the spanwise velocity component to take veer into account.

Turbulence Modeling
The turbulent viscosity in Equation 4 is determined by using a mixing length model. The turbulent viscosity in the atmospheric boundary layer is dependent on a mixing length, m , and a velocity gradient (Pope, 2001). The mixing length for flows in the atmospheric boundary layer is defined by: where κ is the von-Kàrmàn constant, z is the distance from the wall, and λ = 15m is the value reached by m in the free atmosphere (Blackadar, 1962;Sun, 2011). Now the turbulent viscosity is given by: where du dz is the streamwise velocity gradient in the wall-normal direction. Figure 2 shows a profile of turbulent viscosity as a function of height based on a power law velocity profile with shear exponent α = 0.2.

Ground Effect
The presence of the ground will have an effect on the shed vortices. The ground effect is incorporated by applying a symmetry boundary condition at the ground (Bastankhah and Porté-Agel, 2016). This is done by using Equations 7 and 8, with the y-and z-coordinates placed below the ground and inverting the sign of the circulation. This condition has a more dominant effect on the vortices close to the ground as they interact with the boundary.

Superposition of Solutions
Superposing all the effects mentioned earlier leads to a base flow that includes all the features presented. The linearized equation allows us to add features by superposing them onto the velocity components. Notice that in this implementation of the model, the base solutions is a function of only the spanwise directions, y and z, and there is no dependency on the streamwise coordinate, x. However, dependency on the streamwise direction could also be included in the base solution. The vortices do 5 induce motion on each other and the spanwise component of the momentum equation should be used. However, these motions are smaller than the streamwise motions and solving one equation as opposed to three reduces computational cost significantly.

Initial and Boundary Conditions
The initial condition for the perturbation velocity, u , is specified as the yawed disk projected onto a plane normal to the streamwise direction. This shape represents the shape of the wake downstream right after the rotor. The exact shape of the 10 wake near the rotor is much more complicated, and is not taken into account in the present model. The initial profile is set as a uniform distribution of wake deficit (u = −2 a U ∞ ) inside the rotor projected area, where a is the induction. This step function is smoothened using a Gaussian filter to avoid numerical oscillations in the spanwise directions. The lateral boundaries are set to zero perturbation (u = 0) because there is no wake in that region.

15
It is now possible to discretize Equation 4 and solve it numerically. Because the time derivative and pressure-gradient terms were dropped, the equation is parabolic, and it can be solved as a marching problem. The equation is solved by starting from an initial condition at the rotor plane and marching downstream. This is done by using a first-order upwind discretization for the streamwise derivative and second-order finite differencing for the spanwise derivatives. The discrete equation is now: 20 where i, j, and k are the indices denoting the grid points in the x, y, and z coordinates, ∇ 2 are the wall-normal and spanwise components in the Laplacian operator and ν eff is the effective viscosity. The resolution used is on the order of 30-40 grid points per diameter. The computational expense of the algorithm without any optimization is small (∼1-3 s) and can be used to generate curled wake profiles quickly.

25
The proposed numerical method uses a forward-time, centered-space method (Hoffman and Frankel, 2001 for stability from the forward-time, centered-space method (Hoffman and Frankel, 2001). Equation 17 shows a guideline for the stability requirement of the algorithm: This stability criterion is based on a two-dimensional equation and the equation we are solving is three dimensional. However, after testing various conditions, this criterion served as a good guideline for the three-dimensional version of the equation.

5
After some algebraic manipulation, it is possible to show that the maximum grid spacing in the streamwise direction, ∆x, is independent of the spanwise grid resolution. Equation 17 can be written as: After testing the model with several grid spacings, it was found that a resolution on the order of D/∆ ∼ 30 − 40 provided converged results for the model without numerical oscillations.  are based on the cross-stream components of velocity. The overall shape of the curled wake is well-captured by the model.
The overall streamline shapes between the LES and proposed wake model are similar. The pair of counter-rotating vortices are clearly visible in both cases; however, the LES computes streamlines with a more complex shape than the simpler model is capable of capturing. The resulting effect is that, in both cases, the wake deficit cross sections are deformed and curled in a similar fashion. The model does not contain the tower, which is present in the simulation; however, this has a small effect on 25 the wake. Figure 4 shows downstream velocity contours for the case of a LES using an actuator line model under uniform inflow (Howland et al., 2016) and the proposed model including curl and rotation. Again, the streamlines are similar in both cases, but the LES produces more complex patterns. Further, the resultant wake deficit deformation is also similar in both cases, with more deficit remaining at the top of the wake. In this case, asymmetry is observed with respect to the centerline across the y-axis. This asymmetry is caused by the wake rotation induced by the torque applied to the fluid by the rotor. Interestingly, the combination of curl and rotation pushes most of the deficit in a preferred direction. In this case, it is pushed upward in the positive z and negative y directions. This insight can be used to steer the wake accordingly. We also observe that in the LES, the wake diffuses faster than in the model. This is because the model does not take yet into account the turbulence generated by the wake. 5 Figure 5 shows the axial velocity along a horizontal line at hub height for different downstream locations from the simulations in Figures 3 and 4. There is good agreement between the simulations and the model (in general), although some differences can be observed. Near the edges of the wake there is an acceleration in the LES, which is caused by the blockage effect, and the model is not able to capture this. In the case of the ALM, the main difference can be attributed to the different initial condition.
The model assumes a step function, which is different from the wake resolved by the LES. However, the general shape of the wake and its deviation to the sides, is well-captured by the model.

Large-Eddy Simulations using the Actuator Line Model in the Atmospheric Boundary Layer
The framework presented can easily be extended by adding more features. As an example, we present a comparison of the model with a LES of a wind turbine inside a neutral atmospheric boundary layer with a yaw angle, γ = 20 o . The simulation was performed using the Simulator fOr Wind Farm Applications (SOWFA) from the National Renewable Energy Laboratory (Churchfield and Lee, 2012). To add the effects of the atmosphere to the model, a vertical profile of velocity in the streamwise 5 direction is added to the base solution. Also, a linear spanwise velocity component is added to the base solution to take veer into account. The veer profile was chosen as a linear profile that matched the inflow from the LES results. Figure 6 shows the mean streamwise velocity contours for the LES of a wind turbine inside a neutral atmospheric boundary layer and results from the model. In general, there is good agreement between the model and the simulation. We can see that the main difference comes from the wake in the LES diffusing more than in the model. This is expected because the turbulence model does not 10 take into account the turbulence generated by the turbine wake, only the turbulence caused by the velocity gradients in the atmospheric boundary layer. There are also differences resulting from the lateral motion of the vortices. The present model does not take into account the convection of the vortices. This is shown in Figure 6, where the top and bottom vortices stay at the same place when using the model. In reality, these vortices are convected to the sides, as shown in the LES. Figure 7 presents velocity along a line at hub height comparing the model and large-eddy simulations. We observe good agreement in terms of the wake shape and how it moves and curls. There are differences present, which we attribute to the simplifications of the proposed model.

5
We now test the proposed model into the FLOw Redirection and Induction in Steady State (FLORIS) framework (Gebraad et al., 2016;Annoni et al., 2018). We compare the new curled wake model to the two-dimensional Gaussian steering model from Bastankhah and Porté-Agel (2016). In the Gaussian steering model (Bastankhah and Porté-Agel, 2016), the wake deflection due to yaw misalignment of turbines is defined by doing budget analysis on the Reynolds-averaged Navier-Stokes equations.
In the curled wake model, the wake steering is computed by solving a linearized version of the Navier-Stokes equations. by Fleming et al. (2017), that the curled wake mechanism does affect the wake of the second turbine, but current yaw steering models are not able to take this effect into account. Now, we present results for three turbines aligned with the upstream turbine yawed 25 o . The curled wake model predicts deflections up to the third turbine's wake. This approach addresses previous concerns about models not being able to capture wake deflection from downstream turbines (Fleming et al., 2017). We notice that the power predictions from the curled wake 5 overpredict results from LES. This outcome is expected because the vortices resulting from curl do not decay as they travel  downstream. Table 1 shows a comparison between the performance of the Gaussian and curled wake models and simulations performed in SOWFA. The curled wake model provides improvements in predicting power gains for more than two turbines in a row. This outcome is because the vortices from the first turbine are propagated downstream. However, because the vortices do not decay in time, the power may be overpredicted.

Possible Improvements for the Model
The key differences between the model and simulations can be summarized as follows: 1. The vortices caused by the curl effect in the model do not change their position and do not decay. In reality, these vortices induce motion on each other and are advected by the free-stream flow, which may have a lateral component.
2. The turbulence model does not take into account the wind turbine wake. It can only take into account the turbulence 5 from the atmospheric boundary layer background flow. This is why the wake decays faster in the large-eddy simulations compared to the model. The model can be further improved by taking some of these factors into account. However, the present model is able to capture the main dynamics of the curled wake with a reduced computational cost. Further improvements are part of future work.

15
A new model has been proposed to study the aerodynamics of the curled wake. The model solves a linearized version of the Navier-Stokes momentum equation with the curl effect added as a collection of vortices with an elliptic distribution shed from the rotor plane. The model has the ability to include several features of the wake including effects due to yaw ('curl'), wake rotation, a boundary layer profile, and turbulence modeling. The model has been implemented and tested to reproduce curled wake profiles. Good agreement is observed when comparing the model to large-eddy simulations of flow past a yawed 20 turbine using an actuator disk/line model. The model was tested into the FLORIS framework. Good agreement was observed in predicting power extraction by yawing the first turbine in a row of two and three turbines. We observe that the effects of the vortices shed by a yawed turbine propagate for downstream distances longer than the separation between two turbines. This means that a yawed turbine can be used to redirect, not only its own wake, but the wake of other downstream turbines as well. Also, we note that the shed vortices allow for spanwise velocity components, which are vital when considering wake 25 redirection and wind farm controls. The vortices generated are not limited to only yawing, as they can also be used for tilt and combinations of tilt and yaw. This work sets a foundation for a simplified wake steering model to be used in a more general wind farm control-oriented framework.