Brief communication: A double Gaussian wake model

In this paper, an analytical wake model with a double Gaussian velocity distribution is presented, improving on a similar formulation by Keane et al. The choice of a double Gaussian shape function is motivated by the behavior of the near wake region, observed in numerical simulations and experimental measurements. The method is based on the conservation of momentum principle, while stream-tube theory is used to determine the wake expansion at the tube outlet. The model is calibrated and validated using large eddy simulations replicating scaled wind turbine experiments. Results show that the tuned 5 double Gaussian model is superior to a single Gaussian formulation in the near wake region.


Introduction
Analytical engineering wind farm models are low-fidelity approximations used to simulate the performance of wind power systems. A wind farm model includes both a model of the wind turbines and a model of the modifications to the ambient flow induced by their wakes, together with their mutual interactions. Analytical wake models, as opposed to high-fidelity 10 computational fluid dynamics (CFD) models, are simple, easy to implement and computationally inexpensive. In fact, they only simulate macroscopic average effects of wakes and not their small scales and turbulent fluctuations. Engineering wake models find applicability in all those cases that do not need to resolve small spatial and fast temporal scales, such as the calculation of the power production of a wind plant over a sufficiently long time horizon. Such models are also extremely useful in optimization problems, where a large number of simulations might be required before a solution is reached, or where 15 calculations need to be performed on the fly in real-time. Analytical wake models are thus often utilized in wind farm layout planning and in the emerging field of wind farm and wake control (Scholbrock, 2011;Churchfield, 2013;Boersma et al., 2017).
Because of their indisputable usefulness, engineering wake models have been extensively studied in the literature. The Jensen (PARK) formulation is one of the most widely used wake models, to the extent that it is sometimes considered as the industry standard (Keane et al., 2016). The model was first introduced by Jensen (1983), and later further developed by Katic below 3D. This raises the necessity of developing models that accurately represent the wake not only far away from the rotor disk, but also in the near and mid-wake regions. Keane et al. (2016) developed a wake model featuring a double Gaussian velocity deficit distribution, in an attempt to formulate a model that closely resembles observed speed distributions in both the near and far wake regions. In fact, while a single Gaussian function is considered to be a good approximation of the wake velocity distribution in the far wake (Bastankhah 30 and Porté-Agel, 2014, the near wake is better approximated using a double Gaussian distribution. This is due to the presence of two peaks in the speed profiles close to the rotor disk, as also observed in experimental measurements and highfidelity CFD simulations (Wang et al., 2017). The double Gaussian model by Keane et al. (2016), which is referred to as the Keane model in this paper, was developed in a similar fashion to the EPFL Gaussian model (Bastankhah and Porté-Agel, 2014), and was intended to respect the principles of mass and momentum conservation.

35
In this short note, a double Gaussian wake model, based on Keane's model and with emphasis on near wake flow behavior, is derived, calibrated and validated. The present formulation addresses and resolves some issues found in the original implementation of Keane et al. (2016), primarily concerning momentum conservation. In addition, the wake expansion function is defined such that mass flow deficit conservation is achieved at the stream-tube outlet. This paper is organized as follows. The derivation of the double Gaussian wake model is detailed in Section 2, along 40 with the formulation of the wake expansion function. In Section 3, the model is tuned and validated, using both experimental measurements obtained with scaled models in a boundary layer wind tunnel and by numerical results of high-fidelity large eddy simulations (LES). Additionally, the performance of the double Gaussian model is compared to a standard single Gaussian formulation. Finally, concluding remarks and future work recommendations are given in Section 4. Appendix A derives some integrals appearing in the formulation. 45 2 Wake model description

Double Gaussian velocity deficit
The double Gaussian wake model is derived in a similar way to the Frandsen (Frandsen et al., 2006) and EPFL single Gaussian models (Bastankhah and Porté-Agel, 2014). Following their approach, the conservation of momentum principle is applied on an ansatz velocity deficit distribution, which includes an amplitude function. Thereby, an expression for the amplitude is 50 obtained that assures conservation of momentum.
At the downstream distance x from the wind turbine rotor and at the radial distance r from the wake centerline, the wake velocity deficit U ∞ − U (x, r) is modeled as the product of the normalized double Gaussian function g(r, σ(x)), which dictates the spatial shape of the deficit, with the amplitude function C(σ(x)). This yields

55
where U ∞ represents the ambient wind speed and U (x, r) the local flow velocity in the wake. The double Gaussian wake shape function, which is symmetric with respect to the wake center, is defined as where r 0 is the radial position of the Gaussian extrema. The standard deviation of the Gaussian function, noted σ(x), represents the width (cross-section) of each of the two single Gaussian profiles. The wake expands with downstream distance x, causing 60 the transformation of the initial double Gaussian profile in the near wake, through a flat-peak transition region, into a nearly single Gaussian profile in the far wake. The wake expansion function is discussed in further detail in §2.2. x and radial distance r from the wake centerline;ṁ is the mass flow rate through the stream tube; AW is a planar cross-sectional area large enough to contain the wake deficit, and A0 is the rotor disk area; T is the thrust force (by the principle of action and reaction, an equal and opposite force is applied by the rotor onto the flow).
The conservation of momentum principle is now applied on the ansatz velocity deficit distribution, using the amplitude function C(σ(x)) as a degree of freedom. Accordingly, the axial thrust force T is related to the rate of change of momentum p of the flow throughout the stream tube (see Fig. 1), i.e.
whereṁ is the mass flow rate through the stream tube, ∆Ũ an effective wake velocity deficit, ρ the air density and A W a planar cross-section at least large enough to contain the wake deficit. Equation (3) is only valid if there is an equal pressure and negligible flow acceleration at the inlet and outlet sections of the the stream tube and, additionally, if shear forces on the control volume can be neglected. The thrust force T is customarily expressed through the non-dimensional thrust coefficient where A 0 is the rotor swept area.
If the wake velocity, defined in Eqs. (1) and (2), is substituted into the Eq.
(3), one obtains Note that, as the double Gaussian wake expands all the way to infinity, the integral boundary is set accordingly. The integration of Eq. (5), whose details are provided in Appendix A, yields where By substituting the thrust given by Eq. (4) into Eq. (6), and solving the resulting quadratic equation for the amplitude function where d 0 = 4 A 0 /π is the rotor diameter. Both solutions of the amplitude function C(σ) would theoretically lead to the con-85 servation of momentum at all downstream distances. However, the velocity profiles obtained by using C + (σ) are characterized by a negative speed (i.e., in the direction opposite to the ambient flow), and thus C + (σ) is deemed to be a nonphysical solution.
Therefore, the true solution for the amplitude function is C − (σ). In addition, a momentum-conserving solution exists only if M 2 − 1/2 N C T d 2 0 ≥ 0, which might not always be the case for large values of C T . The derived expressions for M and N presented in this paper differ from the results reported in the original publication 90 by Keane et al. (2016), even though all assumptions are identical. The expressions reported in the original paper were also evaluated numerically, yielding nonphysical results that violate the conservation of mass and momentum underlying the formulation.

Wake expansion function
In the previous section, following the conservation of momentum, the shape of the double Gaussian wake deficit has been 95 defined as a function of the Gaussian parameter σ. In this section, a wake expansion function σ(x) is introduced, which is linear with respect to the downstream distance x. By mass conservation, the wake expansion at the position of the stream tube outlet is therefore identified.
In previous work by Frandsen et al. (2006) and Bastankhah and Porté-Agel (2014), stream tube theory was employed to derive an equation for the initial wake width at the turbine rotor. Thereby, the number of tunable parameters of the wake expansion function is reduced, facilitating model calibration. However, this approach includes the assumption that the stream tube outlet is located exactly at the turbine rotor itself, which is hardly true. Results indicate that the derived initial wake width is too large to fit experimental measurements, which in turn requires a model re-tuning (Bastankhah and Porté-Agel, 2014).
In the present work, the stream tube outlet is not assumed to be located at the turbine rotor (x = 0), but at the unknown downstream position x 0 . Therefore, the expansion function is defined as where parameter k * controls the rate of expansion, while represents the wake expansion at x 0 . The wake expansion function is assumed to be linear as in Bastankhah and Porté-Agel (2014).
To derive , mass conservation between the Betz stream tube and the wake model is enforced. Starting from Eq. (3) where U ST is the uniform cross-sectional wake velocity. In this expression, β is the ratio between the stream tube outlet area A ST and the rotor disk area A 0 , which can be expressed as a function of the thrust coefficient C T as

115
At the Betz stream tube outlet (x = x 0 ), the mass flow deficit rate of the double Gaussian (noted DG) wake model writeṡ By equating both mass flow deficits (i.e.,ṁ DG =ṁ ST ), the initial wake expansion can be derived. The solution was computed numerically as a function of the thrust coefficient C T and the spanwise location of the Gaussian extrema r 0 . The resulting surface is presented in Fig. 2. Note that the solution of stream tube theory is defined only in the range 0 ≤ C T < 1, and that 120 tends to infinity as the thrust coefficient approaches the value of 1, due to mass conservation.
The remaining parameters x 0 and k * in the linear wake expansion function expressed by Eq. (9) are not explicitly modeled, and should be tuned based on experimental measurements or high-fidelity simulations, as shown in the next section. Note that the underlying momentum conservation statement expressed by Eq.
(3) has only been defined for ambient pressure. Therefore, the formulation is, strictly speaking, only valid for x ≥ x 0 . However, as pressure recovers rapidly immediately downstream 125 of the rotor, resonable approximations can also be expected for x < x 0 . Finally, k * is expected to depend on atmospheric conditions (Peña et al., 2016) and turbine thrust (Bottasso et al., 2019).

Parameter identification and results
The double Gaussian model proposed in this work has three tunable parameters: k * , x 0 , and k r . Parameters k * and x 0 are used 150 to describe the wake expansion downstream of the turbine rotor, as expressed by Eq. (9). The third parameter, k r , is defined as r 0 = k r /2 and describes the position of the Gaussian extrema. When k r = 1 the curve extrema are located at the tip of the rotor blades, while for k r = 0 the two Gaussian functions coincide at the wake center, leading to a single Gaussian wake profile.
In the original formulation by Keane, k r was fixed at 75% blade-span, as it was argued that most lift is extracted from the flow at this location. In the present work the parameter is tuned based on measurements, as the assumed 75% blade-span 155 position did not lead to satisfactory results. The goal of model calibration is to ensure that the wind velocity profiles match the reference data set as closely as possible.
To this end, the squared error between the wake model and CFD-computed wake profiles is minimized with respect to the free parameters. This estimation problem was solved using the Nelder-Mead simplex algorithm implemented in the MATLAB function fminsearch (Lagarias et al., 1998). To ensure the generality of the results, only a subset of the reference data 160 was used for parameter estimation (namely, the downstream distances 1.7 D, 3 D and 6 D), while the others were used for verification purposes.
The identified parameters are presented in Table 1. The Gaussian extrema were found to be at approximately 53.5% of blade-span (k r = 0.535), while the wake width at x 0 is = 0.23 D. Model calibration also resulted in the positioning of the stream tube outlet at x 0 = 4.55 D, which appears to be a realistic value for the investigated turbine.  The performance of the model clearly depends on the data used for its calibration. Using reference data close to the turbine rotor is important for accurately gauging the positions of the velocity profile extrema, while a wider rage of distances leads to an improved expansion behavior. In the present case, more data from the near wake region (1.7 D and 3 D) was considered in the tuning process than in the far wake (6 D). This leads to a slight overestimation of the velocity deficit at 9 D, which could 175 be attributed to an underestimation of the expansion slope k * . However, tuning the model using the entire set of reference data points leads to only small differences in the identified parameters, which in turn produce wake profiles that are fairly similar to the ones presented here. On the other hand, identifying the model using only data points from the far wake resulted in better fitting results at 9 D, but with either very small values of r 0 -which lead to nearly single-Gaussian profiles-, or with high values of the k * expansion slope -which lead to non-physical solutions of Eq. (8) for the amplitude function in the near wake 180 region, due to excessively small Gaussian widths. Figure 4 depicts with a solid blue line and * symbols the root mean square error (RMSE) between the DG wake model (based on the parameters reported above, obtained from measurements at 1.7 D, 3 D and 6 D) and the reference CFD data as a function of downstream distance. To identify a lower error bound, the DG wake model parameters were also tuned separately at each downstream distance, obtaining seven different local parameterizations. The corresponding RMSE with respect to the CFD solution is reported in the same figure using a dashed blue line. The small difference between the two curves shows that the single (global) parameterization computed using only three distances is only marginally sub-optimal, in the sense that it is very close to the best possible fitting that a double Gaussian shape function can achieve. The plot shows also a slight increase of the difference between the two curves in the far wake region, which can again be attributed to the fact that only one large-distance (6 D) measurement was used in the global tuning.

190
As a comparison, Fig. 4 also shows the results obtained with the EPFL single Gaussian (SG) model (Bastankhah and Porté-Agel, 2014). The SG model was identified with the same procedure and measurements used for the DG model, obtaining SG = 0.3177 and k * SG = 0.0082; the corresponding RMSE with respect to the CFD results is reported in the figure using a solid red line with circles. The lower error bound for the SG model, here again obtained by tuning the parameters independently at each one of the seven available downstream distances, is shown using a dashed red line. As expected, the SG wake model 195 shows a significantly larger RMSE in the near wake region than in the DG case. In fact, here the wake profile is indeed characterized by two peaks, so that the double Gaussian shape function allows for a more precise representation of the actual flow characteristics. Here again it should be noted that the difference between the global and local parameterizations is quite small, which strengthens the conclusion that improved results are due to the ansatz velocity deficit distribution and not to the specific parameterizations.

Conclusions
This short paper presented an analytical double Gaussian wake model. The proposed formulation corrects and improves a previously published model proposed by Keane et al. The shape of the velocity deficit distribution in the wake is described by two Gaussian functions, which are symmetric with respect to the wake center, while the amplitude of the velocity deficit is derived using the principle of momentum conservation. A linear expansion of the width of the Gaussian profiles was assumed, 205 and stream tube theory was used to estimate the conditions at the stream tube outlet.
The model was calibrated and validated using a set of time-averaged CFD simulation results, which replicate wind tunnel experiments performed with a scaled wind turbine in a boundary layer wind tunnel. Results show that the model fits the reference data with good accuracy, especially in the near wake region where a single Gaussian wake is unable to describe the typically observed bimodal velocity profiles. In the far wake, a slight overestimation of the wake deficit could be observed. It 210 is speculated that this might be due to the wake expansion gradient being slightly different in the near and far wake regions.
This claim, however, would need additional work to be substantiated. Future work could extend the wake model to include wake deflection, which could be done in a rather straightforward manner by following Bastankhah and Porté-Agel (2016). In this case, a non-symmetric double Gaussian shape function could be used to model the kidney shape of a deflected wake (Bartl et al., 2018).

215
Code and data availability. A MATLAB implementation of the wake model and the data contained in this article can be obtained by contacting the authors.
Appendix A: Integration of the momentum flux conservation formula Equation (5) can be written as In the following, integrals M and N are solved to obtain Eq. (7a) and (7b).
Competing interests. The authors declare that they have no conflict of interest.