Brief communication: A momentum-conserving superposition method applied to the super-Gaussian wind turbine wake model

Accurate wind farm flow predictions based on analytical wake models are crucial for wind farm design and layout optimization. In this regard, wake superposition methods play a key role and remain a substantial source of uncertainty. Recently, new models based on mass and momentum conservation have been proposed in the literature. In the present work, such methods are extended to the superposition of super-Gaussian-type velocity deficit models, allowing the full wake velocity deficit estimation and design of closely packed wind farms.

1 Introduction

Wind farm design and layout optimization rely on analytical flow models due to a large number of configurations to be evaluated and the computational efficiency of such numerical methods. A typical wind farm flow solver consists of a combination of several sub-models, including a minima a velocity deficit model; a wake-added-turbulence (WAT) model; and possibly a wake deflection model, a blockage model, and a coupled wake–atmospheric-boundary-layer model. The velocity deficit and WAT models usually apply to a single wind turbine: wake superposition methods accumulate the wakes and estimate a wind farm power production for given environmental conditions. Concerning the superposition of velocity deficits, the available methods lacked theoretical justification, until the recent work of Zong and Porté-Agel (2020) and Bastankhah et al. (2021). In these studies, analytical solutions for the velocity deficit superposition are proposed based on the mass and momentum conservation principle. These superposition methods assume Gaussian-shaped velocity deficit profiles. In the present article, the approach of Bastankhah et al. (2021) is extended to super-Gaussian wake velocity deficit profiles. Such models, proposed in Shapiro et al. (2019) and later refined in Blondel and Cathelain (2020), allow for the evaluation of the velocity deficit over the full wake. On the contrary, the Gaussian-based approaches are limited to the far wake. Apart from preventing the appearance of unrepresentable numbers, this allows the study of closely packed wind farm layouts. Indeed, some offshore wind farms such as Lillgrund exhibit small wind turbine inter-distances, down to 3.3 wind turbine diameters. Considering such super-Gaussian velocity profiles together with the Bastankhah et al. (2021) superposition method, an integral has no analytical solution, and an approximation is proposed and compared with the numerical solution. It is also shown in Sect. 3 that the method proposed in Bay et al. (2022) leads to similar results in terms of centerline velocity deficit and is suited for wind farm power predictions. The new superposition method has more robust theoretical foundations than the traditionally used local-linear-sum (LLS) superposition technique (method C in Zong and Porté-Agel, 2020), and its applicability is demonstrated based on the large Horns Rev wind farm.

2 Extension of the Bastankhah et al. (2021) model

2.1 Model derivation

In Bastankhah et al. (2021), the conservation of momentum deficit for multiple wakes takes the form

$$\begin{array}{}\text{(1)}& \underset{\stackrel{\mathrm{̃}}{\mathcal{A}}}{\int }\left(u_{\mathrm{0}}{c}_n{f}_n-{\left(c_n{f}_n\right)}^{\mathrm{2}}-\mathrm{2}{c}_n{f}_n\sum _{i=\mathrm{1}}^{n-\mathrm{1}}{c}_i{f}_i\right)d\stackrel{\mathrm{̃}}{A}\approx {\displaystyle \frac{{\stackrel{\mathrm{̃}}{T}}_n}{\mathit{\rho }}},\end{array}$$

with c_n the maximum velocity deficit of turbine n, i the index of the turbines upwind of turbine n, f_n the self-similar function, $\stackrel{\mathrm{̃}}{A}=\mathit{\pi }{\stackrel{\mathrm{̃}}{r}}^{\mathrm{2}}$ the rotor surface with $\stackrel{\mathrm{̃}}{r}=r/d_{\mathrm{0}}$ and d₀ the wind turbine diameter, ${\stackrel{\mathrm{̃}}{T}}_n$ the thrust force of the unit diameter rotor, u₀ the undisturbed wind velocity, and ρ the fluid density. Based on comparisons to numerical results from a large-eddy-simulation (LES) solver, a modified form was proposed in Bastankhah et al. (2021): the factor “2” in the left-hand side of Eq. (1) is dropped.

Let us consider the original form, Eq. (1). Given a super-Gaussian shape function f_n, a solution for c_n is sought. Following Blondel and Cathelain (2020), the shape function reads $f_i=\exp (-{\stackrel{\mathrm{̃}}{r_i}}^k/\mathrm{2}{\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_i}}^{\mathrm{2}})$, with $k=k\left(\stackrel{\mathrm{̃}}{x}\right)$ the super-Gaussian order and i or n the index of a wind turbine. In the following, we assume that the turbines are sorted from the most upwind to the most downwind, and for two turbines i and n, we have i<n.

Here, as indicated by the tilde, the radius and the super-Gaussian characteristic width are both normalized by the wind turbine diameter d₀, such as $\stackrel{\mathrm{̃}}{r_i}=\sqrt{(y-y_i{)}^{\mathrm{2}}+(z-z_i{)}^{\mathrm{2}}}/d_{\mathrm{0}}$ and $\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_i}={\mathit{\sigma }}_i/d_{\mathrm{0}}$. The following integrals are defined in terms of the gamma function Γ:

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& \underset{\stackrel{\mathrm{̃}}{\mathcal{A}}}{\int }{f}_{n}d\stackrel{\mathrm{̃}}{A}={\displaystyle \frac{\mathit{\pi }}{k}}\mathrm{\Gamma }\left({\displaystyle \frac{\mathrm{2}}{k}}\right){\mathrm{2}}^{\mathrm{2}/k+\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_n^{\mathrm{4}/k},\\ & \underset{\stackrel{\mathrm{̃}}{\mathcal{A}}}{\int }{f}_n^{\mathrm{2}}d\stackrel{\mathrm{̃}}{A}={\displaystyle \frac{\mathrm{2}\mathit{\pi }}{k}}\mathrm{\Gamma }\left({\displaystyle \frac{\mathrm{2}}{k}}\right){\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_n^{\mathrm{4}/k},\mathrm{and}\underset{\stackrel{\mathrm{̃}}{\mathcal{A}}}{\int }{f}_n{f}_{i}d\stackrel{\mathrm{̃}}{A}=\mathcal{I}.\end{array}\end{array}$$

No analytical solution could be found for the last integral, denoted ℐ. Inserting Eq. (2) into Eq. (1) leads to

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{u}_{\mathrm{0}}{c}_n{\displaystyle \frac{\mathit{\pi }}{k}}\mathrm{\Gamma }\left({\displaystyle \frac{\mathrm{2}}{k}}\right){\mathrm{2}}^{\mathrm{2}/k+\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_n^{\mathrm{4}/k}& -c_n^{\mathrm{2}}{\displaystyle \frac{\mathrm{2}\mathit{\pi }}{k}}\mathrm{\Gamma }\left({\displaystyle \frac{\mathrm{2}}{k}}\right){\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_n^{\mathrm{4}/k}\\ & -\mathrm{2}{c}_n\sum _i^{n-\mathrm{1}}{c}_i\mathcal{I}\approx {\displaystyle \frac{{\stackrel{\mathrm{̃}}{T}}_n}{\mathit{\rho }}}.\end{array}\end{array}$$

Using the thrust coefficient $C_{{\mathrm{T}}_n}=\mathrm{8}{\stackrel{\mathrm{̃}}{T}}_n/(\mathit{\pi }\mathit{\rho }{\stackrel{\mathrm{̃}}{d_{\mathrm{0}}}}^{\mathrm{2}}<u_{n-\mathrm{1}}{>}_{(n,x_n)}^{\mathrm{2}})$, the operator $\langle {\rangle }_{(n,x_n)}$ denoting the spatial averaging over the frontal projected area of rotor n at x=x_n, and u the streamwise velocity component as in Bastankhah et al. (2021), one obtains

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{c}_n^{\mathrm{2}}& -c_n{\mathrm{2}}^{\mathrm{2}/k}\left(u_{\mathrm{0}}-\mathrm{2}\sum _i^{n-\mathrm{1}}{\displaystyle \frac{c_i}{{\mathrm{2}}^{\mathrm{2}/k}}}{\displaystyle \frac{k\mathcal{I}}{\mathrm{2}\mathit{\pi }\mathrm{\Gamma }\left(\frac{\mathrm{2}}{k}\right){\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_n^{\mathrm{4}/k}}}\right)\\ & +{\displaystyle \frac{k{C}_{{\mathrm{T}}_n}}{\mathrm{16}}}{\displaystyle \frac{{〈u_{n-\mathrm{1}}〉}_{\left(n,x_n\right)}^{\mathrm{2}}}{\mathrm{\Gamma }\left(\frac{\mathrm{2}}{k}\right){\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_n^{\mathrm{4}/k}}}\approx \mathrm{0}.\end{array}\end{array}$$

Let us introduce a modified integral $\mathcal{J}=k\mathcal{I}/\left({\mathrm{2}}^{\mathrm{2}/k}\mathit{\pi }\mathrm{\Gamma }\left(\frac{\mathrm{2}}{k}\right){\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_n^{\mathrm{4}/k}\right)$. After straightforward manipulations, and assuming u₀=u_h, i.e., a constant, shear-free inflow, the solution for c_n reads

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\displaystyle \frac{c_n}{u_h}}& =\left(\mathrm{1}-\sum _{i=\mathrm{1}}^{n-\mathrm{1}}{\displaystyle \frac{c_i}{u_h}}\mathcal{J}\right)\left({\mathrm{2}}^{\mathrm{2}/k-\mathrm{1}}\right.\\ & \left.-\sqrt{{\mathrm{2}}^{\mathrm{4}/k-\mathrm{2}}-{\displaystyle \frac{k{C}_{{\mathrm{T}}_n}{\left(\frac{{〈u_{n-\mathrm{1}}〉}_{\left(n,x_n\right)}}{u_h}\right)}^{\mathrm{2}}}{\mathrm{16}\mathrm{\Gamma }\left(\frac{\mathrm{2}}{k}\right){\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_n^{\mathrm{4}/k}{\left(\mathrm{1}-\sum _{i=\mathrm{1}}^{n-\mathrm{1}}\frac{c_i}{u_h}\mathcal{J}\right)}^{\mathrm{2}}}}}\right).\end{array}\end{array}$$

The modified form is obtained by using a modified 𝒥 together with Eq. (5) and ${\mathcal{I}}^{\mathrm{mod}}=\mathcal{I}/\mathrm{2}$:

$$\begin{array}{}\text{(6)}& {\mathcal{J}}^{\mathrm{mod}}={\displaystyle \frac{k{\mathcal{I}}^{\mathrm{mod}}}{{\mathrm{2}}^{\mathrm{2}/k}\mathit{\pi }\mathrm{\Gamma }\left(\frac{\mathrm{2}}{k}\right){\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_n^{\mathrm{4}/k}}}.\end{array}$$

2.2 Approximate solutions of the integral ℐ

In a first approach, hereafter referred to as the Gauss approach, one may assume a Gaussian behavior of the model to evaluate 𝒥, as done in Bay et al. (2022). One obtains (see Bastankhah et al., 2021)

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\mathcal{J}}_{\mathrm{Gauss}}^{\mathrm{mod}}=& {\displaystyle \frac{\mathit{\pi }{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_i^{\mathrm{2}}{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_n^{\mathrm{2}}}{{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_i^{\mathrm{2}}+{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_n^{\mathrm{2}}}}\exp \left(-{\displaystyle \frac{{\left({\stackrel{\mathrm{̃}}{y}}_n-{\stackrel{\mathrm{̃}}{y}}_i\right)}^{\mathrm{2}}}{\mathrm{2}\left({\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_n^{\mathrm{2}}+{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_i^{\mathrm{2}}\right)}}\right)\\ & \exp \left(-{\displaystyle \frac{{\left({\stackrel{\mathrm{̃}}{z}}_n-{\stackrel{\mathrm{̃}}{z}}_i\right)}^{\mathrm{2}}}{\mathrm{2}\left({\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_n^{\mathrm{2}}+{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_i^{\mathrm{2}}\right)}}\right).\end{array}\end{array}$$

Figure 1 Mean (a) and maximal (b) relative error of the analytical integrals compared with the numerical evaluation of 𝒥 as a function of the maximal considered super-Gaussian order.

Alternatively, in a second approach hereafter referred to as the kEquiv approach, one may first consider aligned turbines (${\stackrel{\mathrm{̃}}{y}}_i-{\stackrel{\mathrm{̃}}{y}}_n=\mathrm{0}$, ${\stackrel{\mathrm{̃}}{z}}_i-{\stackrel{\mathrm{̃}}{z}}_n=\mathrm{0}$) and later correct the integral for the lateral distance between the rotors using a function $\mathit{\delta }(\stackrel{\mathrm{̃}}{y},\stackrel{\mathrm{̃}}{z})$. This function is identified from the Gaussian solution. A second approximation consists in considering an equivalent super-Gaussian order, $k_{\mathrm{eq}}=\mathrm{1}/\mathrm{2}(k_i+k_n)$. Under these hypotheses, the integral ℐ takes the form

$$\begin{array}{}\text{(8)}& \begin{array}{rl}{\mathcal{I}}_{\mathrm{kEquiv}}^{\mathrm{mod}}=& {\displaystyle \frac{\mathit{\pi }\mathrm{\Gamma }\left(\mathrm{2}/k_{\mathrm{eq}}\right){\mathrm{2}}^{\mathrm{2}/k_{\mathrm{eq}}+\mathrm{1}}{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_i^{\mathrm{4}/k_{\mathrm{eq}}}{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_n^{\mathrm{4}/k_{\mathrm{eq}}}}{k_{\mathrm{eq}}{\left({\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_i^{\mathrm{2}}+{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_n^{\mathrm{2}}\right)}^{\mathrm{2}/k_{\mathrm{eq}}}}}\mathit{\delta }\left(\stackrel{\mathrm{̃}}{y},\stackrel{\mathrm{̃}}{z}\right),\\ & \mathrm{with}\mathit{\delta }\left(\stackrel{\mathrm{̃}}{y},\stackrel{\mathrm{̃}}{z}\right)=\exp \left(-{\displaystyle \frac{{\left({\stackrel{\mathrm{̃}}{y}}_n-{\stackrel{\mathrm{̃}}{y}}_i\right)}^{k_{\mathrm{eq}}}}{\mathrm{2}\left({\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_n^{\mathrm{2}}+{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_i^{\mathrm{2}}\right)}}\right)\\ & \exp \left(-{\displaystyle \frac{{\left({\stackrel{\mathrm{̃}}{z}}_n-{\stackrel{\mathrm{̃}}{z}}_i\right)}^{k_{\mathrm{eq}}}}{\mathrm{2}\left({\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_n^{\mathrm{2}}+{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_i^{\mathrm{2}}\right)}}\right),\end{array}\end{array}$$

and Eq. (6) is used to calculate ${\mathcal{J}}_{\mathrm{kEquiv}}^{\mathrm{mod}}$. Another straightforward approach consists in tabulating the integral values (excluding the $\mathit{\delta }(\stackrel{\mathrm{̃}}{y},\stackrel{\mathrm{̃}}{z})$ function) and linearly interpolating between the data, which is the one retained in practice. For a quantitative comparison, the proposed analytical approximations of the integral 𝒥 are compared to the numerical integration. An interval of $\mathrm{0.2}\le {\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_i$, ${\stackrel{\mathrm{̃}}{\mathit{\sigma }}}_n\le \mathrm{2.5}$ is considered for the characteristic width, and several intervals 2≤k_i, k_n<maxk are considered for the super-Gaussian order, with $\mathrm{2}<{\max }_k\le \mathrm{8}$. The bounding values are representative of the very near wake of a wind turbine under laminar flow conditions and the very far wake ($\stackrel{\mathrm{̃}}{x}>\mathrm{15}{d}_{\mathrm{0}}$) under highly turbulent conditions: the typical operating range of a turbine in a wind farm is covered. Among the characteristic width and super-Gaussian order intervals, 15 values are sampled. Regarding the maximum super-Gaussian order, six equally spaced values are sampled. For each set of four inputs, and for a given maximum super-Gaussian order, the analytical approximations are evaluated, and the error is computed ($\mathrm{error}=\left(\left|{\mathcal{J}}_{\mathrm{Analytical}}\right|-\left|{\mathcal{J}}_{\mathrm{Numerical}}\right|\right)/\left|{\mathcal{J}}_{\mathrm{Numerical}}\right|$). The numerical evaluation is based on the SciPy (Virtanen et al., 2020) “integrate.quad” integration routine and extends from 0 to $\mathrm{6}\cdot \max ({\mathit{\sigma }}_i,{\mathit{\sigma }}_n)$. Then, for each max_k, the average and maximal error are computed and reported in Fig. 1. From these results, the so-called kEquiv method seemingly outperforms the Gauss method and should be preferred. However, it will be shown in Sect. 3 that the impact on the velocity deficit is limited.

3 Results

In a recent study, Lanzilao and Meyers (2022) showed that the super-Gaussian model performed poorly compared with other models: for both Horns Rev and the London Array wind farms, the predicted power production is far below the measured power from supervisory control and data acquisition (SCADA) data. Due to these observations, the model is re-calibrated for the present study. The calibration procedure and the notations used hereafter follow the work of Cathelain et al. (2020). The main difference here lies in the use of a Gaussian profile in the far wake, i.e., ${lim}_{\stackrel{\mathrm{̃}}{x}\to \mathrm{\infty }}k\left(\stackrel{\mathrm{̃}}{x}\right)=\mathrm{2}$. The wake characteristic width is assumed to evolve linearly with axial distance:

$$\begin{array}{}\text{(9)}& \stackrel{\mathrm{̃}}{\mathit{\sigma }}=\left(a_{\mathrm{s}}\mathrm{TI}+b_{\mathrm{s}}\right)\stackrel{\mathrm{̃}}{x}+c_{\mathrm{s}}\sqrt{\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\displaystyle \frac{\mathrm{1}+\sqrt{\mathrm{1}-C_{\mathrm{T}}}}{\sqrt{\mathrm{1}-C_{\mathrm{T}}}}}\right)}.\end{array}$$

The three parameters, a_s, b_s, and c_s, are used for both the super-Gaussian and the Gaussian model. The super-Gaussian order follows an exponential decay function:

$$\begin{array}{}\text{(10)}& k=a_{\mathrm{f}}{e}^{b_{\mathrm{f}}\stackrel{\mathrm{̃}}{x}}+c_{\mathrm{f}}.\end{array}$$

A Gaussian profile is assumed in the far wake; thus, c_f=2. The parameter b_f controls the decay of the super-Gaussian order and is taken as a function of the turbulence intensity. a_f is chosen in such a way that the model fulfills the actuator-disk theory (see Cathelain et al., 2020). This can be enforced numerically using a Newton fixed-point algorithm. To facilitate the implementation, this inversion is performed in a pre-processing stage, and a third-order polynomial fit is proposed:

$$\begin{array}{}\text{(11)}& a_{\mathrm{f}}=-\mathrm{8.2635}{C}_{\mathrm{T}}^{\mathrm{3}}+\mathrm{8.5939}{C}_{\mathrm{T}}^{\mathrm{2}}-\mathrm{8.9691}{C}_{\mathrm{T}}+\mathrm{10.7286}.\end{array}$$

The proposed calibration is not meant to be universal but dedicated to the present study. Future work will be dedicated to a calibration that is reliable in both near- and far-wake regions. Table 1 provides the list of the model coefficients used in this study, obtained using a differential evolution algorithm and a set of nine LES simulations.

Figure 2 Schematic view of the simulated wind farm.

Table 1 Coefficients of the super-Gaussian wake model.

3.1 Comparison against large-eddy simulations from Bastankhah et al. (2021)

For the model comparison, the numerical setup based on the aligned wind farm introduced in Bastankhah et al. (2021) is reproduced. A schematic view of the wind farm is given in Fig. 2. It consist of five columns of three wind turbines, which are all included in the simulation. The velocity deficit studied later is extracted over a line passing through the hub of the wind turbines of the central line, as shown in the schema. The first column in our simulations is located at $\stackrel{\mathrm{̃}}{x}=\mathrm{0}$.

The wind farm flow model builds upon the super-Gaussian model as described in Blondel and Cathelain (2020), using the calibration introduced in Sect. 3. The WAT model proposed in Ishihara and Qian (2018) is employed, together with a so-called maximum-value WAT superposition; see Niayifar and Porté-Agel (2016). A correction factor of 1.25 is applied on the maximum of added turbulence to match the results presented in Bastankhah et al. (2021). Following a convergence study, the rotor disks are discretized based on 12×12 polar grids. Velocity deficit and WAT due to upwind rotor wakes are evaluated at every point on the disk. Then, mean velocity and turbulence intensity are computed and used as an input for the wake models and rotor performance evaluation; i.e., the power and thrust coefficients are given as a function of wind speed. Using a polar discretization, the mesh cells are not uniform in size: the ones located near the edge of the disk are significantly larger than the ones near the hub. Thus, when computing the mean quantities, we use a weighted average whose weighting factors are based on the mesh cell surface. In practice, in the case of aligned rotors, this tends to lower the wake effect since the higher velocity deficit is located at the rotor center where the mesh cell's relative areas are the smallest. A blockage correction based on the vortex cylinder flow model is used; see Branlard and Meyer Forsting (2020). The LLS method is compared to the present method, denoted MC (momentum conserving), with the two approximations for 𝒥^mod, as well as a direct numerical evaluation of the integral, denoted ${\mathcal{J}}_{\mathrm{Num}}^{\mathrm{mod}}$.

Figure 3a shows that, compared with the LLS superposition method, the MC model predicts a lower velocity deficit in both near- and far-wake regions, which is more consistent with LES data. Moreover, the proposed analytical approximations of the integral 𝒥^mod are very close to the numerical approximation in the presented test case. At the rotor planes, discontinuities are observed. This can be partially attributed to the use of the modified momentum-conservation method which improves the results in the far wake as shown in Fig. 3b but does not fully respect the conservation laws, as detailed in Bastankhah et al. (2021). Using the unmodified formulation leads to very high near-wake velocity deficits or even unrepresentable numbers in the presented test case. More than three diameters behind the wind turbine, the results based on the ${\mathcal{J}}_{\mathrm{kEquiv}}^{\mathrm{mod}}$ and ${\mathcal{J}}_{\mathrm{Gauss}}^{\mathrm{mod}}$ approximations are superimposed since the super-Gaussian order is close to 2. These observations validate the approach employed in Bay et al. (2022), despite the higher errors noticed in Fig. 1. In practice, using a tabulated version of the integral is a fast and convenient approach. However, it does not circumvent the approximation based on the rotor distance function, $\mathit{\delta }(\stackrel{\mathrm{̃}}{y},\stackrel{\mathrm{̃}}{z})$, since tabulating the complete integral results in large data files that are time-consuming to load. The global agreement against the LES dataset is satisfying. In the first turbine wake, the hub effect prevents a proper analysis of the results. For the second turbine, a good agreement is obtained with the LLS method, while the MC method underpredicts the velocity deficit. This behavior, as noted in Bastankhah et al. (2021), is a consequence of the application of the modified momentum conservation law. For the following three turbines, a good agreement is obtained.

Figure 3 Centerline velocity deficit in the middle column of the farm; LES data scanned from Bastankhah et al. (2021). (a) Comparison between LES, the LLS method together with the super-Gaussian model (LLS-SG), and the momentum-conserving method together with the super-Gaussian formulations (MC-SG). (b) Comparison between LES and the momentum-conserving method and the Gaussian model, using the original (MC-G_Original) and modified (MC-G) formulations, and a super-Gaussian formulation (MC-SG).

Table 2 Root-mean-square error of the analytical models against the LES results from Bastankhah et al. (2021).

For a more quantitative analysis, the root-mean-square error (RMSE) between the different analytical models from $\stackrel{\mathrm{̃}}{x}=\mathrm{0}$ to $\stackrel{\mathrm{̃}}{x}=\mathrm{30}$ and the LES results are given in Table 2. First, the use of Gaussian wake models leads to a rather high error, due to the inaccuracy in the near wake. This behavior is expected, and we are here using the model outside of its definition domain, i.e., the Gaussian model is a far-wake model. Using the super-Gaussian model, the RMSEs fall below 8 %. Whatever the approximation performed on the 𝒥 integral, the momentum-conserving approach outperforms the LLS method: the RMSEs fall again from approximately 8 % to less than 6 %. Using the ${\mathcal{J}}_{\mathrm{Gauss}}^{\mathrm{mod}}$ approximation, the error is slightly higher compared with ${\mathcal{J}}_{\mathrm{kEquiv}}^{\mathrm{mod}}$ and ${\mathcal{J}}_{\mathrm{Num}}^{\mathrm{mod}}$. One should thus prefer one of these two formulations over the so-called ${\mathcal{J}}_{\mathrm{Gauss}}^{\mathrm{mod}}$ approximation.

3.2 Comparison against large-eddy simulations of the Horns Rev wind farm from Portéé-Agel et al. (2013)

The model predictions are also compared with large-eddy simulations of the Horns Rev wind farm, as presented in Porté-Agel et al. (2013). With a minimal inter-turbine distance of seven diameters, this wind farm can not be considered as closely packed. However, the availability of a large set of large-eddy-simulation results makes it a good candidate for validation purposes. The inflow conditions are based on inflow velocity and turbulence intensity profiles scanned from Porté-Agel et al. (2013). Figure 4 compares the wind farm efficiency η (predicted power divided by theoretical power without wake effect) over a wide range of wind directions θ. We use the LES as a reference to avoid the uncertainties of SCADA measurements, mainly due to the wind direction changes during the 10 min averaging in the available data.

Figure 4 Comparison of the normalized Horns Rev wind farm power output based on LES simulations from Porté-Agel et al. (2013), the LLS method together with the super-Gaussian model (LLS-SG), the momentum-conserving method together with the super-Gaussian model (MC-SG), and the momentum-conserving method together with the Gaussian model (MC-G).

The agreement between the analytical model and the LES dataset is overall good. Differences between the momentum-conserving superposition method and the LLS approach are noticed for wind directions where the wake effects are strong, typically at θ≈{222, 270, 312^∘}. Around such directions, the lower velocity deficits predicted by the MC approach lead to lower wake losses and better efficiency of the wind farm, which is more consistent with the LES data. Both the Gaussian and the super-Gaussian models predict the same wind farm efficiency whatever the wind direction: this is due to the large inter-turbine distances in the Horns Rev wind farm. It confirms that the poor results obtained in Lanzilao and Meyers (2022) for the same wind farm are mostly due to inaccuracies in the model calibration introduced in Cathelain et al. (2020).

For a more quantitative comparison, the RMSEs of the different analytical models against the LES results are provided in Table 3. The LLS method together with the super-Gaussian model has the highest level of error, around 3.7 %, while using the momentum-conserving approach, both with a Gaussian or a super-Gaussian model, causes the RMSEs to fall below 2.5 %. Differences between the two aforementioned models appear in the RMSEs only at the fourth decimal. Considering the large inter-turbine spacing in the Horns Rev wind farm, this was expected, since both models use the same characteristic width, and the inter-turbine distances are large enough to have super-Gaussian orders very close to k=2 at the rotor planes.

Table 3 Root-mean-square error of the analytical models against the LES results from Porté-Agel et al. (2013).

4 Conclusions

In this work, the momentum-conserving wake superposition method proposed in Bastankhah et al. (2021) was extended to super-Gaussian-type velocity deficit models. An integral could not be resolved analytically, and an approximation has been proposed. This approximation is closer to numerical evaluations of the integral that the Gaussian assumption used in Bay et al. (2022). Comparisons against large-eddy simulations of wind farms show a satisfactory agreement, allowing the simulation of large wind farms using the super-Gaussian wake model. Further studies will include an extensive validation of the resulting wind farm flow model, including closely packed wind farms.