3-D shear-layer model for the simulation of multiple wind turbine wakes: description and first assessment

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 estimations of the wind farm efficiency too. In this sense, the interaction between wakes has become a relevant aspect for the definition of a wind farm layout, for the assessment of its annual energy yield and for the evaluation of wind turbine fatigue loads. For this reason, accurate models for multiple overlapping wakes are a main concern of the wind energy community. Existing engineering models can only simulate single wakes, which are superimposed when 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 an alternative, we propose a new shear-layer model that is based on the existing engineering wake models but is extended to also simulate non-axisymmetric wakes. In this paper, we present the theoretical background of the model and four application cases. We evaluate the new model for the simulation of single and multiple wakes using large-eddy simulations as reference. In particular, we report the improvements of the new model predictions in comparison to a sum-of-squares superposition approach for the simulation of three interacting wakes. The lower deviation from the reference considering single and multiple wakes encourages the further development of the model and promises a successful application for the simulation of wind farm flows.


Introduction
When the wind passes through the wind turbine rotor, kinetic energy is extracted from the wind and is converted into electrical power.The reduced kinetic energy is revealed by a wake deficit behind the rotor, i.e. a shear flow with lower speed and higher turbulent fluctuations than in the free flow upstream and sideways.
In this sense, wakes are the main cause of power losses in wind farms (Walker et al., 2016).In addition, wakes hitting a downstream turbine contribute to the increase in the fatigue loads of its components.For these reasons, wake modelling 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 lifetime of wind turbine components.Consequently, more accurate wake models can indirectly contribute to the cost-of-energy reduction thanks to more tailored design of wind turbines and wind farms.
Despite the large progress especially in the numerical modelling, Vermeer et al. (2003) still provide a comprehensive review of traditional wake modelling.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-averaged Navier-Stokes (RANS) equations or large-eddy simulations (LESs) can provide more realistic results because the physics of the flow is resolved up to more refined length and timescales.However, these alternatives have a much higher computational cost and can therefore become prohibitive for design applications.
Engineering tools for estimating wake effects in a wind farm often implement the steady-state, axisymmetric shearlayer approximation of the RANS equations, e.g. the one used in the Ainslie model (Ainslie, 1988).Due to the axial symmetry assumption, only the wind deficit of single wakes or wakes aligned on the same axis as those illustrated in Fig. 1a can be simulated with such models.For the case of wake-turbine or wake-wake interaction, shown in Fig. 1b and c, pragmatic methods are required.In the kinematic model by Katic et al. (1986), the square addition of the individual wake deficits is applied to deal with multiple wakes.In a previous study, Lissaman (1979) proposed the linear addition of the deficits; however, this method tends to overestimate the velocity reduction and could lead to unrealistic flow reversal when many wakes merge.
Machefaux (2015) compared the performance of the linear approach with the one of the square wake addition approach and noticed that the former is to be preferred for wakes of turbines operating at low thrust coefficients, while the latter returns better results in the opposite case.From this observation, he developed a wake superposition model that combines the linear and square additions of single wakes using a weighted average depending on the thrust on the rotor.Crespo et al. (1999) declared that the classical wake superposition methods do not rely on a physical background and, if not handled properly, could lead to unrealistic results.This statement was motivation of this paper.In this regard, we aim to investigate whether a more detailed physical model could improve the simulation of multiple wakes.For this purpose, we pick up the suggestion by Ainslie (1988) to extend his model to the third dimension, dropping the hypothesis of an axisymmetric wake profile; accordingly, we develop the 3-D shear-layer (3DSL) model and test its performance in relation to Ainslie's model and the square addition approach.
For the assessment, we address four cases including a single wake, aligned wakes, and wake-turbine and wake-wake interaction; we use the wind fields extracted from LESs of the same wake conditions as reference and consider the sectionaverage wind speed and the RMSE as figures of merit.

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.

Mathematical definition
The 3DSL model is meant to add the third dimension to the shear-layer approximation of the steady RANS equations for wind turbine wake simulations first described by Ainslie (1988), maintaining all his assumptions but the one of an axisymmetric wake profile.The 3DSL model is intended to simulate the development in the wake of the normalised wind velocity u D , which can be defined as using a representative vertical profile of the inflow wind speed u i , the corresponding hub height value u H and the wind speed u at the desired position.For sake of brevity we will refer to u D simply as wake velocity in the following.The 3DSL model is generally valid starting from a downstream distance at which the pressure gradient in the streamwise direction is negligible.Moreover, the viscous term is not considered and no external forces are applied.Different from other existing shear-layer models, our 3DSL approach is not formulated in a polar coordinate system, but considering a Cartesian frame of reference, i.e. the stream-wise wake velocity u D , the cross-stream, and the vertical wind components v and w are defined along the downstream x, lateral y and upward z axes respectively.In the same way as u D , the two latter wind components are also normalised by u H .
Considering a dimensional analysis (Cebeci and Cousteix, 2005), the steady RANS equation for flows with a shear layer along the cross-stream and vertical component can be simplified to The shear stress terms on the right-hand side of the second line of Eq. ( 2) can be modelled by means of an eddy viscosity closure introducing the eddy viscosities y and z and multiplying them by the corresponding cross-stream and vertical gradients of u D : Further details on the eddy viscosity model are provided in Sect.2.4.
At this point, the system of Eq. ( 2) is still underdetermined.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 a 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 his model, he considered a cylindrical coordinate system defined by the radial coordinate r, the angular coordinate θ and the axial coordinate x.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 the radial velocity v r does not vary along the angular coordinate θ when a constant radial distance r is considered.
The explanation above shows that, like the 3DSL model, the Ainslie model also assumes a potential flow and therefore no vorticity on the cross sections y −z.In the vortex cylinder model of the actuator disc (Burton et al., 2011), the flow field of a wind turbine wake is conservative everywhere but on the surface of the vortex cylinder that encloses the wake, along the root vortex and on the bound vortex sheet swept by the rotor blades.Accordingly, our approximation to a potential flow is reasonable for most of the simulation domain and, even if the real flow is not strictly conservative, the 3DSL model enables us to find one of the solutions for the underdetermined, three-dimensional shear-layer problem that respects the conservation of mass and the momentum balance in the stream-wise direction.
Thanks to Eq. ( 4) and considering that, at each individual vertical cross section, ∂u D /∂x depends only on y and z, the conservation of mass (Eq.2, first line) can be expressed as where g(y, z) = ∂u D /∂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

Numerical implementation
The 3DSL model is implemented using finite difference schemes to obtain the numerical formulation of the physical model defined in Eq. ( 6).Stream-wise gradients are approximated with a forward finite difference scheme, while a central one is used for the gradient in the other directions.The solution of the wind field on each consecutive cross section is accomplished with the following steps: 3. correction of v and w on the previous cross section with the values derived from the definition of (Eq.6, third and forth lines), 4. reiteration of the cycle from step 2 until sufficient convergence of v and w is reached, 5. evaluation of u D on the current section by means of numerical integration of Eq. ( 6), second line.
For the initial condition on the first cross section, a disc actuator model can be applied to estimate u D , while v and w are set to zero.The vertical cross sections y − z are defined by a regular grid with spacing y = z = h; the resolution x along the x axis is evaluated at each cross section.This is needed to accomplish the stability constraints of the numerical solution.In fact, the stream-wise momentum balance (Eq.6, fifth line) is similar to the much simpler problem The solution of this problem with a so-called forward-time central-space finite difference scheme is numerically stable only if µ t/h 2 ≤ 1 4 , where t and h = y = z are the time and space discretisation increments respectively and µ is the diffusive parameter of the problem (Press et al., 2007, chap. 20.5).Inspired by this constraint, we conservatively define the downstream step size at each cross section as The boundary conditions are assigned in two different ways: periodic conditions are applied to solve the Poisson equation (Eq.6, first line), while, for the solution of the stream-wise momentum balance (Eq.6, fifth line), u D is set as in the initial conditions on the boundaries.

Model initialisation
To run simulations with the 3DSL model it is necessary to initialise it with the wind field at the downstream outlet of the induction zone of the rotor, i.e. the region where the pressure field is influenced by the operation of the wind turbine.In fact, as explained in Sect.2.1, the 3DSL model is not valid in the near field behind the rotor where the pressure gradients have a major influence on the flow.Werle (2015) and Madsen et al. (2010) suggested possible methodologies suitable for this purpose.Here, we apply a classic disc actuator approach (Burton et al., 2011) to estimate the initial wake velocity u D,o at the outlet of the induction zone.
We consider the stream tube depicted in Fig. 2 and defined by the cross sections at the inlet, at the rotor and at the outlet of the induction zone.We indicate the corresponding diameters as D i , D r and D o respectively.We use the same subscripts for the section-averaged wind speed U SA and for the stream-wise wind component u.Following the disc actuator theory, we assume that -U SA is homogeneous on each cross section of the stream tube, the induction factor a defined by the thrust coefficient C T as in Eq. ( 16) regulates the evolution of U SA through the stream tube such that According to the conservation of mass of an incompressible flow across the stream tube (see Fig. 2), we can combine with Eq. ( 9) to calculate the inlet and the outlet cross-section diameters of the stream tube: The initial conditions u D,o for the 3DSL model are calculated in three steps: first, we estimate the wind speed u o at the outlet as applying the induction factor a to the inflow wind speed u i on the inlet cross section of the stream tube homogeneously.
Then, the wind field is expanded according to Eq. ( 11).Finally, the initial wake velocity u D,o is given, replacing u with u o in Eq. (1).
To calculate the stream tube cross sections and the corresponding average wind speeds, this method needs to be applied iteratively until convergence.In fact, the induction factor a has to be known.Usually, it can be derived from the thrust coefficient C T associated with the undisturbed wind speed at the inlet of the stream tube according to the wind turbine specification.In the case described here, the undisturbed wind speed is defined as average over the inlet cross section by U SA,i , which in turn is dependent on the induction factor a (see Eq. 11).For this reason, an iterative process is applied starting with the rotor diameter D r as a first guess to approximate the diameter D i of the inlet cross section.
As already mentioned, shear-layer wake models are valid only outside the induction zone.However, Madsen et al. (2010) noticed that the turbulent mixing influences the wake velocity profile already within this region.Therefore, they simulated wakes with their shear-layer model starting from the rotor position.To compensate for the effect of pressure gradients not included in their model but actually present in reality until two to three rotor diameters downstream of the turbine, they applied a linear filter to the ambient eddy viscosity within this range.In the same way, the 3DSL model also first evaluates the wake velocity outside the induction zone to initialise the wake simulation, which in turn starts directly behind the rotor.Then, it applies the linear filter to the ambient eddy viscosity to mimic the effects of the pressure gradients within the near wake.

Eddy viscosity model
In the 3DSL model, the eddy viscosity is evaluated following the approach suggested by Ainslie (1988), who combined the contribution of the wake and of the atmosphere.Experimentally, he found that the proportionality coefficient k = 0.015 links the wake contribution to r i and u a i , which are the characteristic length and velocity scales of turbulent fluctuations within a wake in the cross-wise and vertical directions for i = y, z respectively.Furthermore, he introduced the filter function1 to properly modulate the development of the turbulence generated by the shear layer within the wake.
To model the effect of the atmospheric conditions on the eddy viscosity, Ainslie used the momentum flux profile m (z H /L MO ) as a function of the wind turbine hub height z H and of the Monin-Obukhov length L MO (Dyer, 1974), the von Kármán constant κ and the friction velocity u * .
Based on the definitions above, in the complete eddy viscosity model the first and second addends represent the wake and atmospheric contributions respectively.As explained in Sect.2.3, the filter function F 2 was added following the example by Madsen et al. (2010) to compensate for the pressure effect within the near wake when the 3DSL model is initialised at the rotor position.In Eq. ( 15) we indicate the spatial dependence of the parameters because we want to stress the fact that, thanks to the three dimensions resolved by the model, the eddy viscosity does also not need to be axisymmetric anymore and can be defined locally.For instance, it can vary linearly over the height z or depend on the local strain rates of the wind field, as it will be explained in Sect.2.4.1.

Characteristic scales of turbulent fluctuations within wakes
In the 3DSL model, the characteristic turbulence length scales r y and r z are both approximated with a representative wake radius r(x) derived as a function of the normalised downstream distance x and the thrust coefficient C T using the analytical wake model by Frandsen et al. (2006) and revised by Rathmann et al. (2006) as In case of multiple wakes, only the turbine closest to the cross section considered is regarded in the evaluation of r(x) within the overlapping area.
On each cross section, we define the local characteristic turbulence velocity scale u a i as a function of the position P = (x, y, z).For this purpose, the local characteristic velocity scale is derived with the classic turbulence mixing length theory (Pope, 2000), similarly as in the model by Keck et al. (2012).Accordingly, the turbulent velocity scales are modelled by means of the local strain rates of the wake velocity u D y (P ) = ∂u D ∂y P and u D z (P ) = ∂u D ∂z P together with the turbulence length scale approximated with r(x) in the direction considered.

Multiple wakes
The 3DSL model is suited for simulation of multiple wakes and does not require the addition of individual wakes to resolve the wind field where wakes from different turbines are overlapping.Still, for simulations of multiple wakes it has to deal with the definition of the inflow wind field of a wind turbine hit by other wakes.This is a delicate matter because it generates a sort of conflict between the actuator disc model used for the initialisation of the 3DSL model and the recovery of the wake within the upstream induction zone of the downstream turbine.
The induction zone, that is, the region directly affected by the pressure gradients across the rotor, already begins in the inflow.For instance, the IEC 61400-12-1 standard for power performance measurements suggests measuring the wind speed of the free inflow at least two rotor diameters upstream from the wind turbine.Power performance measurements exclude the case of wind turbines operating in www.wind-energ-sci.net/2/569/2017/Wind Energ.Sci., 2, 569-586, 2017 wakes.We could have followed this indication anyway, but we would have disregarded the recovery of the wake.When a wind turbine operates within a wake, the 3DSL model uses the wind field on the rotor cross section as the inflow in the evaluation of the section-average wind speed U SA,i .Doing this it neglects the effect of the induction zone upstream of the wind turbine, but this is necessary in order to consider the recovery of the wake.Recent studies that investigate how to model the induction zone upstream of the wind turbine rotor (Meyer Forsting et al., 2016) could provide tools to improve this pragmatic approach, but it is out of the scope of the present work.

Wake simulations
In this section we consider single and multiple wind turbine wakes from LES wind fields as a reference to evaluate and compare results from simulations carried out with the 3DSL model and with the Ainslie model as implemented in the wind farm layout software FLaP (Lange et al., 2003).In the latter case we apply the square addition approach to multiple wakes.Accordingly, the total wake velocity resulting from the overlapping of the consecutive wakes is assumed as where u D,sw i is the wake velocity of the ith single wake.The comparison includes three cases of multiple wakes (namely aligned wakes and wake-turbine and wake-wake interactions), preceded by a single-wake simulation.

Test cases and reference wind fields
All the test cases are simulated with the same atmospheric conditions and as a reference consider wakes generated with the LES model implemented in PALM (Raasch and Schröter, 2001), whose solver is coupled with an actuator disc model (Calaf et al., 2010).These LES wind fields deal with wakes from the Siemens SWT-3.6-120wind turbine (120 m rotor diameter D, 90 m hub height z H ). In the test cases two, three and four, the turbines are placed with a consecutive downstream displacement of 6 D and a cumulative separation in the cross-stream direction of 0.0 D, 0.5 D and 1.5 D respectively.These layouts lead to the hub height maps of the wake velocity displayed in Fig. 3.The wind field is evaluated on a uniform grid with a spacial resolution of 10 m (0.083 D) and a total domain size of approximately 20, 5 and 3.5 km along the stream-wise, cross-stream and vertical axes respectively.The reference wind field results from the temporal average of 45 min simulations with a time step close to 1 s.With a roughness length z 0 = 0.002 m and a vertically constant potential temperature, the wind conditions should resemble a typical offshore boundary layer in neutral stratification ( m (z H /L MO ) = 1).The friction velocity u * evaluated fitting the logarithmic profile u = (u * /κ) ln(z/z 0 ) to the average vertical profile of the wind speed on the inflow section is about 0.3 m s −1 .Un-  der these conditions, the hub height wind speed 3.3 D upstream of the first rotor is 8 m s −1 with 5 % turbulence intensity TI.According to this inflow wind speed, the wind turbines are operating in partial load with a thrust coefficient2 C T = 0.858.

Simulations with the shear-layer models
The simulation domains of the 3DSL and of the Ainslie model are different.In the first case, the cross sections are resolved with 111 and 81 points in the lateral and vertical directions respectively, extend from y = −7 to y = +3 D and are 8 D high.The adaptive step in the downstream direction leads to 2291 points from x = 0 to x = 20 D. With these settings, the simulation of three wind turbine wakes takes about 11 s.
In FLaP, we impose the initial condition taking into account the turbulence intensity, the thrust coefficient and the tip speed ratio of the turbine according to Lange et al. (2003).Additionally, for test case 2 and test case 3, we consider the turbulence added by the wake following the empirical formula suggested by Hassan (1992) as reported in Burton et al. (2011).
For the simulation of a single wake with FLaP, 181 points are considered along the downstream direction from x = 2 to x = 20 D; the radial coordinate counts 20 000 points in the range from 0 to 7 D. The enormous number of points in the radial direction is dictated to achieve a convergent result with a downstream step close the one of the LES wind field.This simulation set-up requires a computational time of about 6 s for a single wake.
The computational times reported above refer to simulations performed on one core of a 2.7 GHz standard processor with 16 GB of RAM available, using MATLAB R2016a for the 3DSL model and a compiled Fortran implementation of the Ainslie model for FLaP.

Results
For a quantitative assessment of the results, we sample the wind fields using several virtual turbines of the same type as the one used for the simulations; their rotors are centred on the black dots printed in the wind fields of Fig. 3.An illustrative sketch of a row of the virtual turbine rotors is given in Fig. 4.
With regard to the virtual turbine rotor j , to the corresponding N j grid points and in relation to the reference stream-wise wind component u ref , we analyse the relative deviation of rotor-average wind speed (RAWS) On the one hand, the first two figures of merit are individually considered for each virtual turbine.On the other hand, we calculate the overall values RAWS and E RMS , averaging the absolute values RAWS, j for the former and considering all virtual turbines at once in the calculation of the RMSE for the latter.These overall values are collected in Table 1.
The three methods of evaluation are related, but each has its own specific character.The RAWS is often used as parameter to evaluate the operational state of a wind turbine.In this sense, it is very close to the application field.However, it cannot give precise information about the accuracy of the simulated wind field because inaccurate previsions of the wake velocity could cancel out in the averaging process.The RMSE does not suffer from this problem and can express the accuracy of the simulations with more confidence.Last, we also included the regression analysis in our study because in this way we could see how well the models are correlated to the reference in terms of the coefficient of determination R 2 , and of the corresponding regression line slope A and intercept B. These statistical parameters are included in Table 1 too.
To provide further information on the intermediate results of the simulations, we include figures describing the development of the horizontal and vertical profiles of the wake velocity at different cross sections in Appendix A.

Test case 1: single wake
In the first test case, we address a single wake to assess the general accuracy of the two shear-layer wake models and at the same time to have a term of comparison for the simulation of multiple wakes.
Looking at the results in Fig. 5, the 3DSL model and FLaP tend to have fair and very similar results with values of RAWS (top panels) and E RMS /u H (bottom panels) below 10 % after 6 D downstream.Higher errors occur in the preceding region, especially around the centre of the rotor (y = 0 D) where the RAWS is overestimated.Here, the 3DSL model seems to perform slightly better, in particular from the graphics of E RMS .In the far wake, starting from 12 D the profiles of RAWS and E RMS do not vary much moving downstream.
The difference between the results of the two models perceived in RAWS and E RMS /u H is not found in the overall RAWS RAWS and in the average root mean square error E RMS /u H . Similarly, the results of the regression analysis are essentially the same for the two models.The corresponding scatter plots in Fig. 6 and intercept B suggest that, in general, the two models tend to overestimate the wind speed values, i.e. to underestimate the wake effects.

Test case 2: multiple aligned wakes
Even though the simulation of consecutive aligned wakes with the Ainslie wake model does not require the square addition approach because the wake velocity profiles remain axisymmetric, we apply this approach to be consistent with the other test cases.
The main results are collected in Fig. 7, the top panels of which show that FLaP overestimates the RAWS relative deviation RAWS , particularly around the axis of the real turbine rotor (y = 0 D) where the maximum of the deviation is reached.Moving sideways, the deviation decreases gradually.
Differently, the 3SDL model underestimates the RAWS around the axis of the real turbine rotor and overestimates it around the boundaries of the wakes (y = ±1 D).Also in this case, the highest absolute deviation from the reference is around the axis of the real turbine rotors.
In general, the results give the impression that the 3DSL model simulations are a little more accurate in terms of RAWS.The same conclusion is not evident in the values of the RMSE drawn in Fig. 7c, d.Since in both figures of merit the two models have a very similar behaviour, it is hard to draw a clear conclusion from the comparison.
Contrary to the previous test case, the overall statistics RAWS and E RMS /u H sustain the impression suggested by Fig. 7: the former indicates that 3DSL simulations have a deviation from the reference on average six percentile points lower than FLaP.In contrast, the latter suggests that the two models have the same accuracy in terms of overall RMSE.
The slope and the intercept from the regression analysis (Table 1) show that the 3DSL model approaches an almost perfect regression line.FLaP does not have such good results in these terms, but it is characterised by a lower spread of the data as indicated by the higher coefficient of determination R 2 .This outcome can be explained with the different distribution of the deviation from the reference of the two models (see Fig. 8): on the one hand, the 3DSL model tends to underestimate the lower values of the wind speed, i.e. in the near wake especially in the region around the axis of the real turbine rotor.On the other hand, it tends to overestimate the higher wind speed values around the boundaries of the wake at a further downstream distance.The resulting uneven distribution leads to an almost perfect regression line.Differently, FLaP mainly overestimates the wind speeds in the whole domain, causing a higher intercept and a lower slope of the regression line.The same arguments explain the results of RAWS described before.
Considering all these results, we conclude that the two models simulate the wake of this test case differently, but they have very similar overall performance.FLaP (Ainslie)   The wakes predicted with FLaP and the square addition rule overestimate the rotor equivalent wind speed values almost everywhere and are higher than in the case of the sim-ulation with the 3DSL model.In the results from FLaP, we also notice that the maximal deviation of the RAWS at each cross section is higher than in test case 2 in which the aligned wakes are supposed to be axisymmetric.Furthermore, it increases passing through the third turbine.not observe such behaviour in the 3DSL model, in which the maximum peaks of RAWS have a similar level as in test case 2 on all cross sections.This difference between the two models might be due to the three-dimensional, nonaxisymmetric character of the multiple wakes simulated in this test case, which is better reproducible by the 3DSL model.
The results of the RMSE (Fig. 9c, d) are derived from a different point of view; however, they lead to the same observations.
The overall statistics provide a quantitative summary of the results mentioned above; in particular, the 3DSL model achieves a deviation from the reference RAWS ( RAWS ) 21 percentile points lower than FLaP.Considering the overall RMSE, the spread between the two models is even more acute: in the simulations with the 3DSL model, E RMS is almost 20 percentile points lower than in FLaP simulations.
Here, the regression analysis (see Fig. 10) replicates the results of test case 2, with the difference that, for the 3DSL model simulations, the slope A and the intercept B of the regression line are not so close to their ideal values 1 and 0 respectively.In turn, the coefficient of determination R 2 is a little higher indicating less scatter of the data.For the simulations with FLaP we observe a remarkable increase in the intercept, which indicates a larger overestimation of the wind speed, meaning a larger underestimation of the wake effects.
4.4 Test case 4: multiple wakes with 1.5 D lateral separation Due to the increased cross-stream separation between the three turbines considered in this test case, the flow seems composed by single wakes.The results presented in Fig. 11 are therefore comparable with those of test case 1, but with an amplified difference between the performance of the two models.In fact, with regard to the reference, both the deviation of the RAWS and the RMSE evaluated for FLaP are clearly higher than the ones evaluated for the 3DSL model.
The corresponding overall values give a measure of this difference: both the deviation RAWS of FLaP and the average root mean square error E RMS /u H are more than 10 percentile points larger than the ones of the 3DSL model.
The regression analysis displayed in Fig. 12 provides results close to those of test case 1 for both models, apart from the intercept, which in test case 4 is lower for the simulation with the 3DSL model, while it is higher for the FLaP simulation.

Discussion
In Sect.4, we compared two shear-layer wake models with a different level of detail in the physical description of the flow.The results are not always easy to interpret because in some cases one model was accurate where the other was not and vice-versa.We dealt with this problem by estimating different figures of merit that are generally in agreement.This temporarily solves the conflict between the applicative point of view of the RAWS and the more wind-field-oriented approach of the RMSE.The object of comparison was the performance of the two models in the simulation of multiple wakes.In this regard, the figures of merit are generally in favour of the 3DSL model.This is a positive outcome of our research and encourages the further development of this new model.Nonetheless, the two models provided similar results for axisymmetric wakes (test case 1 and test case 2).This points out the advantage given by the third dimension resolved by the 3DSL model.In fact, in the other test cases, i.e. when multiple wakes have a lateral separation, the additional details in the physical description of the flow implemented in the 3DSL model seem to improve the results.
Despite the different performances, we found similar deficiencies in the two models.This particularly regards the flow of single wakes near the rotor cross section as indicated by the results of test case 1 and in test case 4. Furthermore, the regression analysis and the scatter plot indicate the tendency to overestimate the wind velocity in the same cases.More in detail, it is possible to notice that the lowest wind speed in single wakes near the rotor was underestimated, while further downstream there was a general overestimation of the wind speed in the wake.This indicates that single wakes might have recovered in the transition between near and far wake too fast.The analysis of the individual wake profiles could help to understand this interpretation and at the same time could provide hints about how to deal with this issue.In many cases, a possible solution could be provided by different eddy viscosity models.In this sense, the three-dimensional domain of the 3DSL model offers the possibility to develop proper spatial distributions of these quantities, while the axisymmetric two-dimensional models would have more limits in the accomplishment of this task.

Conclusions
This paper investigated the possibility of improving the simulation of multiple wakes with engineering wake models such as the commonly used Ainslie model (Ainslie, 1988), implemented for instance in the wind farm layout software FLaP (Lange et al., 2003).In this regard, the paper presented a new wake shear-layer model (3DSL) that can deal with nonaxisymmetric flows and is therefore suitable to simulating multiple wakes at once.Differently, when the Ainslie model is applied in a wind farm, the flow of multiple wakes is evaluated, superimposing the deficit of the individual wakes according to a linear or square addition approach.To allow the simulation of multiple wakes without the superposition of the individual wakes, the 3DSL model abandons the hypothesis of an axisymmetric wake assumed by Ainslie (1988) and adds a third dimension to the simulation domain.In order to do this, it assumes a potential flow on the vertical cross sections.
In a benchmark with LESs as a reference case, we considered four test cases and compared wake simulations performed with FLaP and with the 3DSL model.The assessment was based on the average wind speed on the rotor of www.wind-energ-sci.net/2/569/2017/Wind Energ.Sci., 2, 569-586, 2017 several fictive turbines and on the corresponding RMSE.The two models provided similar results when they simulated axisymmetric wakes, but the 3DSL model performed better in the test cases including non-axisymmetric wakes.In part, this might be one of the advantages of the third dimension included in the 3DSL model.Since only a few test cases using wakes simulated within LESs were addressed here, these results cannot be generalised.For the same reason, we cannot make any statement about how these results could affect the estimation of the annual energy yield of a wind farm.Nevertheless, we are confident that the additional details in the physical description of the wake flow implemented in the 3DSL model can in general offer new possibilities for improving the simulation of single and multiple wakes at an affordable computational cost.

Appendix A: Normalised wake wind velocity profiles
The comparison of the 3DLS model and of FLaP with the reference wind field from large-eddy simulations (LESs) reported in this paper deals with figures of merit representative for integral results and therefore cannot provide detailed results about the output of the two models at specific positions in the simulation domain.For this reason we show the downstream development of the wake velocity for the test cases analysed in this paper in this appendix.

Figure 2 .
Figure 2. Sketch of the stream tube used to describe the disc actuator approach.The dashed lines represent the inflow, rotor and outlet cross sections, which are indicated with the subscripts i, r and o in the definition of the diameter D and the section-average wind speed U SA .

Figure 3 .
Figure 3. Colour map of the hub height wake velocity u D evaluated for the test cases from the large-eddy simulations (LESs).The black dots indicate the position of the virtual turbine rotors used to compare the simulation results.

Figure 4 .
Figure 4. Illustrative sketch of the rotors of the virtual turbines considered to assess the performance of the engineering models in relation to the large-eddy simulations.
et al.: 3-D shear-layer simulation model the linear regression of the stream-wise wind component values on the grid points within the rotor area.

Figure 5 .
Figure 5. Test case 1 (single wake).Deviation of the rotor-average wind speed RAWS (a, b) and of the root mean square error E RMS (c, d) evaluated in relation to the large-eddy simulation wind field for the simulation with the 3DSL model and with FLaP.

Figure 6 .Figure 7 .
Figure 6.Test case 1 (single wake).Scatter plot and corresponding regression line of the wind speed derived from the 3DSL model (a) and from the FLaP wake simulations in relation to the reference wind field calculated with large-eddy simulations (LESs).

Figure 8 .
Figure 8. Test case 2 (multiple aligned wakes).Scatter plot and corresponding regression line of the wind speed derived from the 3DSL model (a) and from the FLaP wake simulations in relation to the reference wind field calculated with large-eddy simulations (LESs).

Figure 9 .
Figure 9. Test case 3 (multiple wakes with 0.5 D lateral separation).Deviation of the rotor-average wind speed RAWS (a, b) and of the root mean square error E RMS (c, d) evaluated in relation to the large-eddy simulation wind field for the simulation with the 3DSL model and with FLaP.

Figure 10 .
Figure 10.Test case 3 (multiple wakes with a 0.5 D lateral separation).Scatter plot and corresponding regression line of the wind speed derived from the 3DSL model (a) and from the FLaP wake simulations in relation to the reference wind field calculated with large-eddy simulations (LESs).

Figure 11 .
Figure 11.Test case 4 (multiple wakes with 1.5 D lateral separation).Deviation of the rotor-average wind speed RAWS (a, b) and of the root mean square error E RMS (c, d) evaluated in relation to the large-eddy simulation wind field for the simulation with the 3DSL model and with FLaP.

Figure 12 .
Figure 12.Test case 4 (multiple wakes with a 1.5 D lateral separation).Scatter plot and corresponding regression line of the wind speed derived from the 3DSL model (a) and from the FLaP wake simulations in relation to the reference wind field calculated with large-eddy simulations (LESs).

Figure A1 .Figure A2 .
Figure A1.Test case 1 (single wake).Downstream development of the vertical (a) and horizontal (b) profiles of the wake velocity evaluated along the wind turbine rotor axis from the 3DSL model and from the FLaP simulations and from the reference large-eddy simulation (LES) wind field.

Figure A3 .
Figure A3.Test case 3 (multiple wakes with 0.5 D lateral separation).Downstream development of the vertical (a) and horizontal (b) profiles of the wake velocity evaluated along the common axis of the wind turbines from the 3DSL model and FLaP simulations and from the reference large-eddy simulation (LES) wind field.The position of the profiles considered is illustrated in the top right corner of (b).

Figure A4 .
Figure A4.Test case 4 (multiple wakes with 1.5 D lateral separation).Downstream development of the vertical (a) and horizontal (b) profiles of the wake velocity evaluated along the common axis of the wind turbines from the 3DSL model and FLaP simulations and from the reference large-eddy simulation (LES) wind field.The position of the profiles considered is illustrated in the top right corner of (b).

Table 1 .
Overall performance of the 3DSL model and FLaP (Ainslie model) in relation to the reference large-eddy simulation wind field.Namely, the average deviation RAWS of the rotor-average wind speed, the total root mean square error E RMS , the coefficient of determination R 2 ), the corresponding regression line slope A and intercept B are included.