Shear layer approximation of Navier-Stokes steady equations for non-axisymmetric wind turbine wakes: Description, verification and first application

The number of turbines installed in offshore wind farms has strongly increased in the last years and at the same time the need for more precise estimation of the wind farm efficiency. For this reason, the wind energy community could benefit from more accurate models for multiple wakes. Existing engineering models can only simulate single wakes, which are superimposed if they are interacting in a wind farm. This method is a practical solution, but it is not fully supported by a physical background. The limitation to single wakes is given by the assumption that the wake is axisymmetric. As alternative, we propose a new shear model which is based on the existing engineering wake models, but is extended to simulate also non- axisymmetric wakes. In this paper, we present the theoretical background of the model and two application cases. First, we proved that for axisymmetric wakes the new model is equivalent to a commonly used engineering model. Then, we evaluated the improvements of the new model for the simulation of a non-axisymmetric wake using a large eddy simulation as reference. The results encourage the further development of the model, and promise a successful application for the simulation of multiple wakes.


Introduction
When the wind passes through the wind turbine rotor, kinetic energy is extracted from the wind and is converted into electrical power. This process generates a wake which propagates downstream. Wakes can be described as shear flows with lower speed and higher turbulence fluctuations than in front of the rotor. In this sense, wakes are the main cause of power losses in wind farms [1]. Besides that, wakes hitting a turbine contribute to increase the fatigue loads of its components. For these reasons, wake modeling plays a major role in the definition of the layout of wind farms, in the evaluation of their annual energy yield and in the estimation of the life-cycle of wind turbine components. Consequently, more accurate wake models can indirectly contribute to the cost of energy reduction due to more tailored design of wind turbine and farms.
Vermeer et al. [2] give a comprehensive review about traditional wake modeling. Most of the engineering models described in their work evaluate the wind field of a single wake and combine the individual results in case of mutual interaction. More sophisticated Computational Fluid Dynamics (CFD) such as Reynolds Average Navier-Stokes (RANS) or Large Eddy Simulations (LES) can deal with wake superposition better and provide more realistic results. However, these alternatives have a much higher computational cost and therefore can become prohibitive for commercial applications.
Commercial codes for estimating wake effects in a wind farm often implement the steadystate, axisymmetric shear-layer approximation of the RANS equations, e.g. the one used in the Ainslie model [3]. Due to the axial symmetry assumption, only the wind deficit of single wakes can be simulated with such models. In the kinematic model by Katik et al. [4], the quadratic superposition of the individual wake deficits is applied to overcome this limitation and to be able to deal with multiple wakes. In a previous study, Lissaman [5] proposed a linear superposition, however this method tends to overestimate the velocity deficit and could lead to unrealistic flow reversal when many wakes merge.
Machefaux [6] compared the performance of the linear with the quadratic wake superposition approach and noticed that the former is to prefer for wakes of turbine operating at a low thrust coefficients, while the latter returns better results in the opposite case. From this observation, he developed a wake superposition model which combines the linear and quadratic superposition of single wakes using a weighted average depending on the thrust on the rotor.
As a final remark about wind farm modelling, Crespo et al. [7] declared that the classical wake superposition assumption does not rely on a physical background and, if not handled properly, could lead to unrealistic situations. This statement gives the motivation of this paper. Herewith, we aim to describe an innovative engineering model which is able to deal with the superposition of wakes and is based on physical principles without relying on a superposition method. Additionally, the contribution aims to verify the model's applicability against some test cases.
The model introduced in this paper consists of a three-dimensional steady-state shear layer (3DSL) wake model which can simulate non-axisymmetric wakes. We validated the model considering an axisymmetric wake and comparing its results with the ones calculated with the Ainslie model implemented in the code FLaP [8]. Finally, we applied the 3DSL model to a non-axisymmetric wake and compared the results with the corresponding LES.

Model description
In the following the theoretical background of the 3DSL model is provided along with the description of its numerical implementation. Moreover, it is explained how to evaluate the parameters needed to apply the model.

The mathematical definition
The 3DSL model implements a shear layer approximation of the steady RANS equations following the approach described by Lange et al. [8]. It is intended to model the wind turbine wake deficit u defined as using the inflow wind speed u i , the wake wind speed u w and the hub height inflow wind speed u H . The 3DSL model is valid starting from a downstream distance where the pressure gradient in the stream-wise direction is negligible. Moreover, the viscous term is not considered and no external forces are applied. Differently from other existing shear layer models, the 3DSL approach is not formulated in a polar coordinate system, but considering a Cartesian frame of reference, i.e. the stream-wise, cross-wise and vertical wind components u, v and w are defined along the x, y and z axis respectively. Considering a dimensional analysis [9] the steady RANS equation for flows with a The shear stress terms on the right hand side of the second line of equation 2 can be modelled by means of an eddy-viscosity closure introducing the eddy viscosities y , z and multiplying them by the corresponding lateral and vertical gradients of u: Further details on the eddy-viscosity model are provided in section 2.3.
At this point, the system of equation 2 is still underdeterminated. To balance the unknown variables and the equations, we assume that the wind components v and w define a conservative vector field in all the cross-sections y − z. A potential function Φ can therefore be defined such that Concerning multiple wakes, this assumption does not imply any limitation since the vector field resulting from the superposition of conservative vector fields is still conservative. However, this assumption limits the domain of possible solutions. For instance, swirling wakes in which the tangential velocity is inversely proportional to the distance from the rotation axis are accepted, while wakes rotating as a rigid body are not.
The hypothesis of a potential flow is implicit in the axial symmetry imposed by Ainslie. In its model he considers a cylindrical coordinate system defined by the radial axis r, the angular coordinate θ and the horizontal axis x. Here, the corresponding velocity vector field V (r, θ, x) = (v r , v θ , u) is conservative only if ∇ × V = 0. Considering the individual cross-section planes at a certain x coordinate, it implies that ∂v r /∂θ − ∂v θ /∂r = 0. This equation is always satisfied by the Ainslie model in which the tangential velocity v θ is neglected and radial velocity v r is the same at each radius r.
Thanks to equation 4 and considering that ∂u/∂x depends only on x and y, the conservation of mass (equation 2, first line) can be expressed as where f (y, z) = −∂u/∂x. This formulation is a second order elliptic partial differential equation of the Poisson type, which can be solved numerically.
Considering the aforementioned assumptions, the final formulation of the 3DSL model can be summarised as

The numerical implementation
The 3DSL model is implemented using finite difference schemes to obtain the numerical formulation of the physical model defined in equation 6. On the vertical cross-sections y − z, the grid points are equally spaced, while the downstream step size along the x axis is evaluated at each cross-section. This is needed to accomplish the stability constraints of the numerical solution.
The numerical problem can be solved iteratively for well defined boundary conditions. Generally u = 1 is set on the boundaries of the computational domain. Furthermore, it is possible to reproduce the ground effects by imposing ∂Φ/∂z = 0 on the lower boundary of the domain for the solution of the Poisson equation (equation 6, first line). This effect is not included in the common implementation of the Ainslie model. In the 3DSL model, the same settings can be reproduced imposing periodic boundary conditions.

Eddy viscosity model
In the 3DSL model, the eddy viscosity is evaluated as suggested by Lange et al. [8] with superposing the contribution of the wake (first addend) and of the atmosphere (second addend). The parameters describing the wake contribution to the eddy viscosity are the empirical parameter k and the filter function F (x) [3], which are included to modulate the development of turbulence generated by the shear layer within the wake deficit. Last, the parameters r y,z (x) and u a (x) are meant to represent scales of wake deficit extensions and amplitude respectively.
The parameters appearing in equation 7 to model the effect of the atmospheric conditions on the eddy viscosity are the momentum flux profile Φ m (z H /L M O ) as function of the wind turbine hub height z H and of the Monin-Obukonov length L M O [10], the Von Karman constant κ and the friction velocity u * . In a neutrally stratified atmosphere, u * is proportional to the standard deviation of the streamwise wind velocity [11]. According to experimental data [11,12], it can be approximated with where u H , u H std and T I are the inflow velocity at hub height, its standard deviation and the corresponding ambient turbulence intensity.

Wake parametrization
To generalise this model for a two-dimensional non-axisymmeric wake deficit, a representative wake deficit amplitude u a and its representative radii r y and r z in the corresponding y and z directions are regarded.
In order be consistent with Lange et al. [8] and at the same time to allow dealing with nonaxisymmetric wake deficits, the metrics representative of the wake deficit at each cross section are derived from the most representative Gaussian surface resulting from a least-square fit of to the wake deficit. The function G(y, z) is modulated by the deficit depth u a and is centered about the point g 0 = (y 0 , z 0 ). Using this function enables to identify the rotation α of the symmetry axes y and z of the Gaussian surface, along with the corresponding standard deviations σ y and σ z Two more steps are needed to calculate the representative radii r y and r z used in equation 7: First, r y and r z are defined as the distance from g 0 along y and z , at which the velocity deficit recovered to 97.13%. Then, r y and r z are calculated with a trigonometric projection onto the y and z axis of the frame of reference.
In the case of an axisymetric wake deficit, u a coincides with the center-line deficit u c−l , the rotation angle α vanishes and a single radius r can be used instead of r y and r z .
Using a Gaussian surface to estimate the wake parameters is not a limitation for the eddy viscosity model described in section section 2.3. Local irregular gradients might characterise the wind deficit of single and multiple wakes. Since G(y, z) is smooth and convex, it cannot capture these occurrences. However, it can identify the representative scales of the wind deficit amplitude and spatial extensions in the sense of integral values as required by an eddy viscosity model which prescribes a homogeneous eddy viscosity over the whole simulation domain.

Model application
An evaluation of the 3DSL model is presented in this section with regards to the Ainslie model applied to an axisymmetric wake first, and then to a non-axisymmetric wake simulated with LES.

Verification against the Ainslie model ( FLaP)
As explained in section 2, the 3DSL model using periodic boundary conditions is supposed to be equivalent to the Ainslie model for the simulation of axisymmetric wakes. In this regard, we considered six test cases to verify the 3DSL model against the the Ainslie model implemented in the wind farm layout software FLaP [8].
The test cases deal with axisymmetric wakes of an NREL offshore 5-MW baseline wind turbine (hub height z H and rotor diameter D of 80 m and 126 m respectively) [13] operating in different atmospheric conditions. They correspond to a neutral atmospheric stratification, i.e.
We evaluated the wake deficit on a 10D long domain with a 5Dx5D cross-section. We used a fixed downstream incremental step in the FLaP simulations, providing 17 downstream positions x. Differently, the 3DSL model implements a dynamic step size to ensure the numerical convergence of the solution. The resulting number of downstream positions computed for each test case is listed in table 1. The radius r and the center-line deficit u c−l were chosen as basis for the evaluation and are addressed in figure 1(a) and (b) respectively. The two models provided very similar results with a maximal discrepancy in the radius lower than 0.04% and about 0.2% for what concern the center-line deficit. In all test cases, the 3DSL model estimated a larger wake radius and, consistently with the mass conservation, a lower center-line deficit. For both parameters the difference between the two models accumulates fast until about 3D downstream. Afterwards, in the case of the radius it reduces, while it seems to converge to a constant value in the case of the center-line deficit. The oscillations of the curves is related to the rougher discretisation of the x axis used in FLaP.

Benchmark against Large Eddy Simulations
To evaluate the advantages provided by the 3DSL model in comparison to the Ainslie model, we addresed a non-axisymmetric wake with a near neutrally stratified atmosphere in a benchmark. The wake comes from a LES coupled with an enhanced actuator disc model based on blade element lift and drag that is able to create heterogeneous thrust of the rotor and rotation of the wake. The case analysed here reproduces a measured wake situation at the alpha ventus wind farm in the German North See. The LES were setup as described in detail by Vollmer et al. [14] and were re-run modelling the Adwen wind turbine installed in alpha ventus instead of the NREL offshore 5-MW baseline wind turbine [13]. This change led to slight differences in the wake with an even closer replication of the measured wake situation.
The simulated wind turbine has a hub height z H of 90 m and a rotor diameter D of 116 m. A top view of the horizontal wind speed on the LES domain at 87.5 m is given in figure 2(a), where a solid black segment indicates the rotor position, the red box and the dash-dot black lines identify the area considered for the benchmark and the corresponding cross-sections (2D upstream, and 2D, 3D, 5D, 8D downstream the rotor).
The inflow conditions 2D upstream the rotor are described in figure 2(b). The vertical profile of the horizontal wind speed and direction are plotted in the left and right panel respectively. The latter reveals a wind direction change of about 10 o across the rotor which explains the asymmetry of the wake: The high shear in the cross-stream component of the wind stretches the wake deficit. In the simulations, we defined the stream-wise direction as the wind direction at hub height. The initial wake deficit 2D downstream was calculated with equation 10 where u i and u w are the stream-wise wind speeds from the LES on the cross-sections 2D upstream (a) Wind speed on the hub-height horizontal plane and on the 2D downstream crosssection.
(b) Vertical profiles of the inflow wind speed and direction sampled 2D upstream the wind turbine rotor. and downstream the rotor respectively. For the normalisation, u H was chosen as the maximum value of the stream-wise wind speed. The wind field that defines the initial condition for this test case is not suitable for FLaP because it is not axisymmetric. Thus, the formula defining u 0 in equation 10 was fitted to the actual initial deficit from the LES and used to start the simulation with FLaP. This solution is meant to provide the axisymmetric wake which can best approximate the LES initial deficit.
The results are compared visually in figure 3(a). From the sequence of downstream crosssections, it is possible to notice, on one side, how the 3DSL model differently from FLaP tends to maintain the same form of the deficit given by the LES; on the other side, how the implementation of the diffusion effects in both models cannot reproduce the downstream development of the LES deficit. Consequently, an overall slower recovery rate is observed. To provide a quantitative understanding of these results, the same parameters used in equation 9 served for the evaluation of the models under analysis. In this sense, G(y, z) was fitted to the wind deficit on the downstream cross-sections and the results were collected in the plots of figure 4(a). From the top panel, it is evident how the 3DSL model is able to catch the asymmetry of the LES wake deficit, but it cannot reproduce its downstream stretching. In fact r y and r z from the 3DSL model increase at a low rate and with a very slight convergence, while the corresponding values from the LES diverge at higher rates. From this perspective, it seems that the 3DSL model diffuses the wake almost homogeneously. Considering FLaP, the development of the two radii collapse to a single line which increases between r y and r z resulting from the 3DSL model.
The rotation angle in the middle panel reveals some differences between the LES and the 3DSL model. The one of the former is stable, while the one of latter follows it until about 4D then it start to increase. These effect could be a consequence of the diverse expansion resulting from the two simulations. No rotation was to considered in FLaP since it assumes that the wake deficit is axisymmetric.
Concerning the amplitude of the wake deficit, the two engineering models provided an almost identical evolution. They could follow fairly good the amplitude depth evaluated for the LES and until 5D, then they overestimated it.   The previous observations are reflected also in the regression analysis presented generally in figure 3(b), which includes the scatter plots of the wake deficit evaluated at corresponding points with the two engineering models and by means of LES. The solid red line in each frame is the result of the regression analysis. At the first downstream cross-section considered (3D), it is very close to the dashed black line representing a perfect correlation. Moving downstream, these lines separate progressively. The exact results of the regression analysis including additional crosssections are summarized in figure 4(b). A part from the coefficient of determination R 2 of the regression line, its slope and offset are very similar for both engineering models. The increase of the offset downstream means that the wake deficit recovery is slower in the engineering models than in LES. This attests the results given by figure 4(a) for the wake deficit depth. While the offset increases, the slope decreases suggesting that the compared deficits could converge to different values.
Although the slope and the offset of the two models in question are very similar, this does not mean that they have the same performances. In fact, the coefficients of determination R 2 evaluated with the regression analysis for FLaP and the 3DSL model are clearly in favour of the latter. This result is in agreement with the impression suggested by figure 3(b), that the points identifying the deficit evaluated with the 3DSL model are less scattered around the regression line than the ones calculated with FLaP.

Discussion
In the comparison with FLaP, the 3DSL model performed very well simulating an axisymmetric and non-axisymmetric wakes. In contrast, the deformation and the recovery of the wake resulting from the 3DSL model and from the LES are only in partial agreement. A different modeling of the transport-diffusion processes is one of the responsible for the different propagation of the wake downstream.
The vertical veer -i.e. a high shear in the cross-stream wind speed-influences the deformation of the wake downstream in the LES. The current implementation of the 3DSL model gives only  the possibility to set an initial condition on the stream-wise wind speed. This means that no veer can be imposed on the simulated flow to improve the results. Future developments could investigate the possibility to initialize v and w with a conservative wind field compliant with the assumption implied by the shear layer formulation.
A more sophisticated diffusion model could also partly compensate for the effect of the vertical veer. In the 3DSL, the global eddy-viscosity which controls the diffusion process at each step is evaluated on the basis of the representative amplitude of the wake deficit and the corresponding radius. Alternatively, the eddy-viscosity linked to the shear in the wake and to the ambient turbulence could be evaluated locally on the basis the local gradients of the wake deficit [15] and of the vertical profile of the horizontal wind speed [16] respectively.

Conclusions and outlook
This paper presents a new wake shear layer model (3DSL) that can deal with non-axisymmetric flows such as multiple wakes partially overlapping. It demonstrates that the 3DSL model is equivalent to the commonly used Ainslie model [3] implemented for instance in the wind farm layout software FLaP [8]. Furthermore, it provides an evaluation of the two models using a wake deficit from LES as reference wind field.
To allow the simulation of overlapping wakes, the 3DSL model abandons the assumptions of an axisymmetric wake implemented for example in the Ainslie model and add a third dimension to the simulation domain. In order to do this, it assumes a potential flow on the vertical cross-sections.
The validation against the Ainslie model considered a wind turbine operating at high and low thrust, with different turbulence conditions. The wake radius and center-line deficit resulting from the simulation of the two codes are in agreement with a maximal variation of about 0.04% and 0.21% respectively.
The benchmark against the LES proved that the 3DSL model can simulate non-axisymmetric wakes. However, it was not possible to fully catch the wake expansion and deformation found in the reference wind field. This does not mean that the 3DSL failed. It just makes evident the different level of simplification used in the two simulations confronted.
The proposed model is just at the beginning of its development. It positively passed the first tests and the results indicated the direction to follow for further improvements such that in the near future the simulations of wake interaction within a wind farm could benefit from the purely physical approach adopted in the 3DSL model.