FLOW Estimation and Rose Superposition (FLOWERS): an integral approach to engineering wake models

Annual energy production (AEP) is often the objective function in wind plant layout optimization studies. The conventional method to compute AEP for a wind farm is to first evaluate power production for each discrete wind direction and speed using either computational fluid dynamics simulations or engineering wake models. The AEP is then calculated by weighted-averaging (based on the wind rose at the wind farm site) the power produced across all wind directions and speeds. We propose a novel formulation for time-averaged wake velocity that incorporates an analytical integral of a wake deficit model across every wind direction. This approach computes the average flow field more efficiently, and layout optimization is an obvious application to exploit this benefit. The clear advantage of this new approach is that the layout optimization produces solutions with comparable AEP performance yet is completed 2 orders of magnitude faster. The analytical integral and the use of a Fourier expansion to express the wind speed and wind direction frequency create a relatively smooth solution space for the gradient-based optimizer to excel in comparison to the existing weighted-averaging power calculation.

1 Introduction

The layout of a wind plant is a primary design element that influences its performance. Optimizing the layout can be thought of as a wake avoidance problem, wherein turbines are placed such that they avoid the wakes from other turbines as much as possible. Power losses from wake interactions can be on the order of 10 %–20 % in wind farms (Barthelmie et al., 2007, 2009). When turbines are placed within about 3 rotor diameters, these power losses can be as high as 40 % (Stanley et al., 2019); even within 15 diameters, wake interactions are non-negligible (Meyers and Meneveau, 2012).

In controls and optimization applications, the wake velocity deficit is approximated with low-fidelity analytical models. The classical top-hat model parameterizes the wake expansion rate and computes the wake deficit as a function of downstream position (Jensen, 1983). Improvements on this approach aim to replace the discrete boundaries of the top-hat model with a continuous profile, such as the Jensen-cosine (Tian et al., 2015, 2017) and Gaussian (Bastankhah and Porté-Agel, 2014) models. More involved engineering models such as the TurbOPark model (Nygaard et al., 2020) account for wake combination and wake-added turbulence more formally in their formulation. The curled wake model is a mid-fidelity numerical model derived from the Reynolds-averaged Navier–Stokes equations (Martínez-Tossas et al., 2019, 2021). The trade-off to explicitly capturing more of the flow physics is the added complexity, both in the calibration of additional parameters and in computational cost. These steady-state wake models are well-suited to estimate wake velocity in simulations with a single wind direction. However, computing average wake velocities or energy production for different wind speeds and directions requires averaging the results of multiple simulations. This process is cumbersome, especially with the more complex models like the curled wake model.

Layout optimization studies leverage these low-fidelity models to approximate the wake velocity within the wind farm. Turbines are placed to minimize wake interactions and thereby maximize annual energy production (AEP) of the plant. Gradient-based optimization algorithms leverage the derivative of the objective function to choose search directions for optimal solutions, while gradient-free optimization only evaluates the objective function (thereby avoiding its derivatives) and is useful for discontinuous and noisy functions. Gradient-free algorithms are common practice in industry for small wind farms (Stanley et al., 2021), but they scale poorly with additional degrees of freedom (Herbert-Acero et al., 2014; Ning and Petch, 2016). Gradient-based optimization, on the other hand, is more robust in systems with a larger number of design variables. The simplest structure for the design variables is to assign the position of each turbine independently (Feng and Shen, 2015; Gebraad et al., 2017). A strategy to reduce the number of design variables is to restrict the layout to a grid (González et al., 2017; Perez-Moreno et al., 2018) or use a combination of placement along the farm boundary and a grid on the interior (Stanley and Ning, 2019). These approaches reduce the cost of the layout optimization study, especially for larger wind farms. However, they restrict the freedom and flexibility of the wind farm developer and produce simplistic layouts that can underperform in practice. Research that addresses the calculation of AEP has focused on statistical methods to improve the efficiency of estimating this quantity (King et al., 2020; Padrón et al., 2019), or on defining the discrete inputs of the wake model (such as the wind rose and power and thrust curves) with analytical functions (Murcia et al., 2015). These mathematical approaches are promising but leave room for a physics-based framework to modify the formulation of AEP.

AEP is an integral quantity. The total power production of a wind plant is calculated based on the wind speed flowing through each turbine. For a single wind speed and direction, this procedure is straightforward. Figure 1 illustrates the flow field around a single turbine, as is often studied in wake modeling problems. The velocity contour plots represent the weighted-averaged flow distribution based on the wind roses shown below each respective contour plot. Figure 1a shows the velocity distribution when the wind rose contains only one predominant wind direction. If there is more than one wind direction, the flow fields from each speed–direction bin must be averaged with weights equal to their normalized frequency as seen in Fig. 1b and c. For example, in Fig. 1b the frequency of the 270^∘ wind direction is greater than that for 225^∘ (and the free-stream wind speed is held constant), and so the velocity contour plot shows the wake directed horizontally with a stronger velocity deficit compared with the angled wake. This procedure is extended across every discrete wind speed–direction bin with the contribution to the sum weighed by the frequency of that bin. Figure 1c illustrates the averaged flow field and how the higher-frequency wind directions manifest as more pronounced wake deficits in the contour plot. The AEP of the wind plant is therefore a numerical integral of the power as a function of wind speed and direction:

$$\begin{array}{}\text{(1)}& \mathrm{AEP}=\sum _{i=\mathrm{1}}^{\mathrm{WD}}\sum _{j=\mathrm{1}}^{\mathrm{WS}}P\left({\mathit{\theta }}_i^{\prime },U_{\mathrm{\infty }}{\left({\mathit{\theta }}_i^{\prime }\right)}_j\right)f_i{f}_j,\end{array}$$

where P is the total power of the wind farm as a function of wind direction θ^′ and free-stream wind speed U_∞, and f_i and f_j  are the frequency of each discrete wind direction and speed, respectively.

Figure 1 (a, d) Velocity contour plot of flow through a wind turbine for a single wind direction. (b, e) Averaging effect of two wind directions. (c, f) The annually averaged velocity flow field. Note that the wind roses (d–f) display the frequency of each wind direction with a constant wind speed of 8 m s⁻¹ for every direction.

The inspiration for the FLOW Estimation and Rose Superposition (FLOWERS) flow field model is to analytically compute the average wake velocity given the frequency and magnitude of the wind speed for every direction. Since the average wake velocity is conceptualized similarly to AEP, extending the FLOWERS approach to calculating AEP would be straightforward. We hypothesize that the analytical integration will considerably reduce the computational cost of average wake velocity and AEP calculations compared to the numerical integration.

In this paper, we first derive the equations for the time-averaged wake velocity and a new formulation for AEP in Sect. 2, including its application in the wind plant layout optimization problem. In Sect. 3, the AEP calculations from FLOWERS are compared to the numerical integration method, which is the standard approach in the NREL FLORIS model (NREL, 2021). Finally, in Sect. 4, the performance of the FLOWERS AEP in the optimization of wind farm layouts is compared against the numerical integration-based optimizer.

2 Mathematical formulation

2.1 Time-averaged wake speed

To derive a mathematical formulation for the time-averaged flow distribution, we use the classical Jensen (top-hat) wake deficit model (Jensen, 1983):

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{U(x,y)}{U_{\mathrm{\infty }}}}=\mathrm{1}-{\displaystyle \frac{\mathrm{1}-\sqrt{\mathrm{1}-C_{\mathrm{T}}}}{(kx+\mathrm{1}{)}^{\mathrm{2}}}}W(x,y),\end{array}$$

where x and y are the streamwise and spanwise positions, respectively, normalized by the rotor radius with the origin at the turbine location. In this coordinate frame, the wind is coming from the negative x direction. We only consider the 2D plane at hub height such that wake deficit is not a function of a vertical position z. C_T is the thrust coefficient, which is realistically a function of the local inflow speed; we simplify it to be a function of the free-stream wind speed U_∞ such that it is uniform across all turbines in the wind farm but still varies around the wind rose. U is the wake speed, k is the wake expansion coefficient, and W(x,y) represents the Jensen wake region: $W(x,y)=\mathrm{1}$ if $\left|y\right|\le kx+\mathrm{1}$ and x≥0 and is zero elsewhere. This geometry is illustrated in Fig. 2a.

Figure 2 (a) A schematic of the top-hat wake model in Cartesian coordinates with the wind aligned with the x direction. (b) The wake region in polar coordinates, rotated by a wind direction θ^′, with the wake boundary defined by the angle θ_c.

We transform Eq. (2) from Cartesian to polar coordinates denoted by r and θ, where x=rcos (θ) and y=rsin (θ). We allow the wind direction θ^′ to be variable in Eq. (3):

$$\begin{array}{}\text{(3)}& \begin{array}{rl}U(r,\mathit{\theta },{\mathit{\theta }}^{\prime })=& U_{\mathrm{\infty }}\left({\mathit{\theta }}^{\prime }\right)\left[\mathrm{1}-{\displaystyle \frac{\mathrm{1}-\sqrt{\mathrm{1}-C_{\mathrm{T}}\left(U_{\mathrm{\infty }}\left({\mathit{\theta }}^{\prime }\right)\right)}}{\left(kr\cos \right(\mathit{\theta }-{\mathit{\theta }}^{\prime })+\mathrm{1}{)}^{\mathrm{2}}}}\right.\\ & \left.W\left(r\cos \right(\mathit{\theta }-{\mathit{\theta }}^{\prime }),r\sin (\mathit{\theta }-{\mathit{\theta }}^{\prime }\left)\right)\right].\end{array}\end{array}$$

To clarify, θ is the angular position in polar coordinates where we wish to compute the average wake velocity, and θ^′ is the wind direction defined within that coordinate frame. We integrate across all wind directions θ^′ to compute the average wake speed at a given location (r,θ). In doing so, we weight each wind direction with its frequency f(θ^′), so the weighted-averaged wind speed denoted by $\stackrel{\mathrm{‾}}{U(r,\mathit{\theta })}$ is written as

$$\begin{array}{}\text{(4)}& \begin{array}{rl}\stackrel{\mathrm{‾}}{U(r,\mathit{\theta })}=& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\underset{-\mathit{\pi }}{\overset{\mathit{\pi }}{\int }}{U}_{\mathrm{\infty }}\left({\mathit{\theta }}^{\prime }\right)f\left({\mathit{\theta }}^{\prime }\right)\\ & \left[\mathrm{1}-{\displaystyle \frac{\mathrm{1}-\sqrt{\mathrm{1}-C_{\mathrm{T}}\left(U_{\mathrm{\infty }}\left({\mathit{\theta }}^{\prime }\right)\right)}}{\left(kr\cos \right(\mathit{\theta }-{\mathit{\theta }}^{\prime })+\mathrm{1}{)}^{\mathrm{2}}}}W(r,\mathit{\theta },{\mathit{\theta }}^{\prime })\right]\mathrm{d}{\mathit{\theta }}^{\prime }.\end{array}\end{array}$$

We define two new variables to simplify this expression. Let $u=\mathit{\theta }-{\mathit{\theta }}^{\prime }$ be the angular position relative to the wind direction. Also, the wake region $W(r,\mathit{\theta },{\mathit{\theta }}^{\prime })$ is zero for θ outside of the wake region defined by θ_c: $\mathit{\theta }={\mathit{\theta }}^{\prime }\pm {\mathit{\theta }}_c$. In our coordinate system, the wake geometry is defined by the line $\sin \left({\mathit{\theta }}_c\right)=k\cos \left({\mathit{\theta }}_c\right)+\mathrm{1}/r$, as shown in Fig. 2b. We solve this equation for r>1 since we are interested in the wake speed at positions outside of the rotor area:

$$\begin{array}{}\text{(5)}& \tan \left({\mathit{\theta }}_c\right)={\displaystyle \frac{\frac{\mathrm{1}}{r}+k\sqrt{\mathrm{1}+k^{\mathrm{2}}-{\left(\frac{\mathrm{1}}{r}\right)}^{\mathrm{2}}}}{-\frac{k}{r}+\sqrt{\mathrm{1}+k^{\mathrm{2}}-{\left(\frac{\mathrm{1}}{r}\right)}^{\mathrm{2}}}}}.\end{array}$$

We must include the interaction between multiple wakes for this formulation to be useful in systems of several turbines. The wake velocity deficit is defined as the difference between the free-stream velocity and the wake velocity. Equation (4) can be split into these two components. Then, we linearly superimpose the deficits to form a relation for the total velocity deficit caused by all turbines:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}\stackrel{\mathrm{‾}}{U(r,\mathit{\theta })}& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\underset{-\mathit{\pi }}{\overset{\mathit{\pi }}{\int }}{U}_{\mathrm{\infty }}\left({\mathit{\theta }}^{\prime }\right)f\left({\mathit{\theta }}^{\prime }\right)\mathrm{d}{\mathit{\theta }}^{\prime }\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}}\sum _i\underset{-{\mathit{\theta }}_c}{\overset{{\mathit{\theta }}_c}{\int }}{U}_{\mathrm{\infty }}\left({\mathit{\theta }}_i-u\right)f\left({\mathit{\theta }}_i-u\right)\\ & \left[{\displaystyle \frac{\mathrm{1}-\sqrt{\mathrm{1}-C_{\mathrm{T}}\left(U_{\mathrm{\infty }}\left({\mathit{\theta }}_i-u\right)\right)}}{{\left(k{r}_i\cos \left(u\right)+\mathrm{1}\right)}^{\mathrm{2}}}}\right]\mathrm{d}u,\end{array}\end{array}$$

where r_i and θ_i are the relative radius and polar angle with respect to each turbine position ($r_i=\sqrt{{\left(x-x_i\right)}^{\mathrm{2}}+{\left(y-y_i\right)}^{\mathrm{2}}}$ and $\tan \left({\mathit{\theta }}_i\right)=y_i/x_i$, where x_i and y_i represent the position of the center of the ith turbine). The relevant information about the wind conditions is the wind direction θ^′, the wind speed U_∞(θ^′), and the wind direction frequency f(θ^′). Note that for simplicity a single average wind speed is used for each wind direction. These quantities are specified for a particular location by the wind rose. In practice, the wind rose is a discrete data set in which wind directions (and their associated average speeds and frequencies) are binned. We define $g\left({\mathit{\theta }}^{\prime }\right)=\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}{U}_{\mathrm{\infty }}\left({\mathit{\theta }}^{\prime }\right)f\left({\mathit{\theta }}^{\prime }\right)$ and $h\left({\mathit{\theta }}^{\prime }\right)=\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}\left[\mathrm{1}-\sqrt{\mathrm{1}-C_{\mathrm{T}}\left(U_{\mathrm{\infty }}\left({\mathit{\theta }}^{\prime }\right)\right)}\right]U_{\mathrm{\infty }}\left({\mathit{\theta }}^{\prime }\right)f\left({\mathit{\theta }}^{\prime }\right)$. If we expand g(θ^′) and h(θ^′) with Fourier series,

$$\begin{array}{}\text{(7)}& {\displaystyle }& {\displaystyle }g\left({\mathit{\theta }}^{\prime }\right)={\displaystyle \frac{c_{\mathrm{0}}}{\mathrm{2}}}+\sum _{n=\mathrm{1}}^N{c}_n\cos \left(n{\mathit{\theta }}^{\prime }\right)+d_n\sin \left(n{\mathit{\theta }}^{\prime }\right),\text{(8)}& {\displaystyle }& {\displaystyle }h\left({\mathit{\theta }}^{\prime }\right)={\displaystyle \frac{a_{\mathrm{0}}}{\mathrm{2}}}+\sum _{n=\mathrm{1}}^N{a}_n\cos \left(n{\mathit{\theta }}^{\prime }\right)+b_n\sin \left(n{\mathit{\theta }}^{\prime }\right),\end{array}$$

then the wind rose is defined continuously rather than discretely. The Fourier coefficients a₀, a_n, b_n, c₀, c_n, and d_n can be easily found for a given g(θ^′), h(θ^′), and N. For a wind rose with B wind direction bins, the maximum number of terms in this discrete Fourier transform is $N=\mathrm{ceiling}(B/\mathrm{2})$, where “ceiling” indicates that we round up to the nearest integer. Also, we approximate the fraction in the second term in the right-hand side of Eq. (6) using a second-order Taylor expansion:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}\stackrel{\mathrm{‾}}{U(r,\mathit{\theta })}=& \underset{-\mathit{\pi }}{\overset{\mathit{\pi }}{\int }}g\left({\mathit{\theta }}^{\prime }\right)\mathrm{d}{\mathit{\theta }}^{\prime }-\sum _i\underset{-{\mathit{\theta }}_c}{\overset{{\mathit{\theta }}_c}{\int }}h\left({\mathit{\theta }}_i-u\right)\\ & \left[{\displaystyle \frac{\mathrm{1}}{(kr+\mathrm{1}{)}^{\mathrm{2}}}}+{\displaystyle \frac{kr{u}^{\mathrm{2}}}{(kr+\mathrm{1}{)}^{\mathrm{3}}}}\right]\mathrm{d}u.\end{array}\end{array}$$

The first integral in the right-hand side of Eq. (9) represents the weighted average of the free-stream velocity, denoted by $\stackrel{\mathrm{‾}}{U_{\mathrm{\infty }}}$ hereafter. Using Eq. (7),

$$\begin{array}{}\text{(10)}& \begin{array}{rl}\stackrel{\mathrm{‾}}{U_{\mathrm{\infty }}}& =\underset{-\mathit{\pi }}{\overset{\mathit{\pi }}{\int }}g\left({\mathit{\theta }}^{\prime }\right)\mathrm{d}{\mathit{\theta }}^{\prime }=\underset{-\mathit{\pi }}{\overset{\mathit{\pi }}{\int }}\left[{\displaystyle \frac{c_{\mathrm{0}}}{\mathrm{2}}}+\sum _{n=\mathrm{1}}^N{c}_n\cos \left(n{\mathit{\theta }}^{\prime }\right)\right.\\ & \left.+d_n\sin \left(n{\mathit{\theta }}^{\prime }\right)\right]\mathrm{d}{\mathit{\theta }}^{\prime }=c_{\mathrm{0}}\mathit{\pi }.\end{array}\end{array}$$

The second integral on the right-hand side of Eq. (9) represents the average wake velocity deficit $\stackrel{\mathrm{‾}}{\mathrm{\Delta }{U}_i(r_i,{\mathit{\theta }}_i)}$. Using Eq. (8),

$$\begin{array}{ll}{\displaystyle }\stackrel{\mathrm{‾}}{\mathrm{\Delta }{U}_i\left(r_i,{\mathit{\theta }}_i\right)}=& {\displaystyle }\underset{-{\mathit{\theta }}_c}{\overset{{\mathit{\theta }}_c}{\int }}\left[{\displaystyle \frac{a_{\mathrm{0}}}{\mathrm{2}}}+\sum _{n=\mathrm{1}}^N{a}_n\cos \left(n\left({\mathit{\theta }}_i-u\right)\right)\right.\\ {\displaystyle }& {\displaystyle }\left.+b_n\sin \left(n\left({\mathit{\theta }}_i-u\right)\right)\right]\left[{\displaystyle \frac{\mathrm{1}}{{\left(k{r}_i+\mathrm{1}\right)}^{\mathrm{2}}}}+{\displaystyle \frac{k{r}_i{u}^{\mathrm{2}}}{{\left(k{r}_i+\mathrm{1}\right)}^{\mathrm{3}}}}\right]\mathrm{d}u.\end{array}$$

Solving the above integral yields

$$\begin{array}{}\text{(11)}& \begin{array}{rl}\stackrel{\mathrm{‾}}{\mathrm{\Delta }{U}_i\left(r_i,{\mathit{\theta }}_i\right)}& ={\displaystyle \frac{a_{\mathrm{0}}{\mathit{\theta }}_c\left[k{r}_i\left({\mathit{\theta }}_c^{\mathrm{2}}+\mathrm{3}\right)+\mathrm{3}\right]}{\mathrm{3}{\left(k{r}_i+\mathrm{1}\right)}^{\mathrm{3}}}}+\sum _{n=\mathrm{1}}^N\\ & {\displaystyle \frac{\mathrm{2}\left[a_n\cos \left(n{\mathit{\theta }}_i\right)+b_n\sin \left(n{\mathit{\theta }}_i\right)\right]}{{\left[n\left(k{r}_i+\mathrm{1}\right)\right]}^{\mathrm{3}}}}\\ & \left[\sin \left(n{\mathit{\theta }}_c\right)\left\{n^{\mathrm{2}}\left[k{r}_i\left({\mathit{\theta }}_c^{\mathrm{2}}+\mathrm{1}\right)+\mathrm{1}\right]-\mathrm{2}k{r}_i\right\}\right.\\ & \left.+\mathrm{2}n{\mathit{\theta }}_{c}k{r}_i\cos \left(n{\mathit{\theta }}_c\right)\right].\end{array}\end{array}$$

The time-averaged wake speed $\stackrel{\mathrm{‾}}{U(r,\mathit{\theta })}$ at a given location (r,θ) is the difference between the free-stream and the wake-averaged velocity:

$$\begin{array}{}\text{(12)}& \stackrel{\mathrm{‾}}{U(r,\mathit{\theta })}=\stackrel{\mathrm{‾}}{U_{\mathrm{\infty }}}-\sum _i\stackrel{\mathrm{‾}}{\mathrm{\Delta }{U}_i\left(r_i,{\mathit{\theta }}_i\right)},\end{array}$$

where $\stackrel{\mathrm{‾}}{U_{\mathrm{\infty }}}$ and $\stackrel{\mathrm{‾}}{\mathrm{\Delta }{U}_i(r_i,{\mathit{\theta }}_i)}$ can be found from Eqs. (10) and (11), respectively.

2.2 Annual energy production

The power P produced by a turbine is a function of the incoming wind speed:

$$\begin{array}{}\text{(13)}& P\left(U\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{C}_{\mathrm{P}}\left(U\right)\mathit{\rho }A{U}^{\mathrm{3}},\end{array}$$

where ρ is the density of air (assumed to be constant), C_P is the power coefficient, A is the swept area of the rotor, and U is the speed at the location of the center of the turbine when it is not present. The average turbine power would be $\stackrel{\mathrm{‾}}{P\left(U\right)}=\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }A\stackrel{\mathrm{‾}}{C_{\mathrm{P}}\left(U\right)U^{\mathrm{3}}}$. However, the integral to compute this average is intractable in our formulation. The power coefficient is dependent on the local wind speed, but U is not known as a function of wind direction prior to the integration; to obtain it, we would need to calculate the wake speed independently for each wind direction, which negates the purpose of the FLOWERS approach. Also, including the nonlinear U³ term in the integral introduces complications with the wake superposition and the definition of the independent wake regions. We make two simplifications to compensate: C_P is calibrated as a function of the average wake speed $\stackrel{\mathrm{‾}}{U}$ (which for each turbine is a constant across all wind directions), and we substitute  ${\stackrel{\mathrm{‾}}{U}}^{\mathrm{3}}$ into Eq. (13) instead of evaluating $\stackrel{\mathrm{‾}}{U^{\mathrm{3}}}$. Therefore, the AEP for a given turbine is

$$\begin{array}{}\text{(14)}& \mathrm{AEP}\left(\stackrel{\mathrm{‾}}{U}\right)\propto C_{\mathrm{P}}\left(\stackrel{\mathrm{‾}}{U}\right){\stackrel{\mathrm{‾}}{U}}^{\mathrm{3}}.\end{array}$$

2.3 Layout optimization problem

We apply this novel formulation of time-averaged wake velocity and AEP to the wind plant layout optimization problem. For M turbines in a wind plant, there are 2 M design variables to independently specify the position of each turbine in 2D space. The objective function is the total AEP of the plant, which we aim to maximize. The turbines are constrained within a specified boundary with a minimum separation of 2 rotor diameters. Gradient-based optimization is performed with the sequential least squares programming (SLSQP) algorithm. This optimizer solves a minimization problem, so we technically minimize −AEP:

$$\begin{array}{}\text{(15)}& \begin{array}{ll}{\min }_{x_i,y_i}& -\mathrm{AEP}\left(x_i,y_i,U_{\mathrm{\infty }}\left({\mathit{\theta }}^{\prime }\right),f\left({\mathit{\theta }}^{\prime }\right)\right)\\ \mathrm{subject}\mathrm{to}& \mathrm{boundary}\mathrm{constraints}\\ & S_{ij}\ge \mathrm{2}\mathrm{D},\end{array}\end{array}$$

where x_i and y_i are the center of the turbine i in Cartesian coordinates, and S_ij represents the separation between the centers of turbines i and j.

The primary benefits of FLOWERS lie in its suitability to drive layout optimization as a wake avoidance problem, despite the fact that simplifications made to develop FLOWERS might induce some errors in the predicted magnitude of AEP. The optimizer relies on the objective function to provide a quantitative metric to compare possible solutions; in this case with gradient-based optimization, the ratio of the objective function evaluated for two different solutions is of importance. In other words, the objective function's output itself is not necessarily critical as long as the mapping between inputs and outputs in the function remains consistent. If we think of the layout optimization problem as a wake avoidance problem, then the objective function must be able to approximate wake magnitudes and downstream influence to minimize their interactions. Turbines aligned with predominant wind directions and turbines with close spacing will reduce AEP in the FLOWERS optimization, just as it will in the numerical integration-based optimization. Wake avoidance can be achieved despite a less accurate estimate of AEP because the factors that cause positive or negative changes to AEP are still present. The gradient throughout the optimization space will be different because the objective functions are not identical. However, with a sufficiently strict convergence criterion, we predict that the FLOWERS optimization will still find a similar quality solution to the numerical integration technique. A more accurate AEP estimate can be added as a final post-processing step once the layout optimization is complete.

3 AEP comparison

3.1 Example cases

We start by comparing the AEP estimates for an illustrative test case of three turbines aligned with a predominant wind direction. The AEP for the numerical integration approach is computed using the Jensen wake deficit model with the same nominal (i.e., based on ambient turbulence intensity and excluding wake-added turbulence) wake expansion coefficient as in FLOWERS (k=0.05). The rotor diameter is D=126 m throughout this paper based on the NREL 5 MW turbine (Jonkman et al., 2009). Figure 3 illustrates the wind rose and farm layout for this case and the flow fields generated from the FLOWERS and Jensen models.

Figure 3 Annually averaged flow field comparison for three turbines aligned with the predominant wind direction with 6 D spacing. The free-stream wind speed is a constant 8 m s⁻¹ (a). The number of wind direction bins B used in (c) is B=72, and the number of Fourier coefficients used to plot results in (b) is N=37.

Figure 4 Flow field comparison for a 60-turbine wind farm with 6 D spacing with a realistic wind rose (corresponding to an offshore site in the northwest US) (a). Again, the number of wind direction bins B in the numerical method (c) is B=72 , and the number of Fourier coefficients in FLOWERS (b) is N=37.

The FLOWERS AEP is 2.9 % lower than the result from the numerical integration approach. Substantial wakes only exist for three discrete wind directions in this example, so the profile of U versus θ^′ is mostly uniform except for a sharp decrease in these aligned orientations. As a result, the difference between ${\stackrel{\mathrm{‾}}{U}}^{\mathrm{3}}$ and $\stackrel{\mathrm{‾}}{U^{\mathrm{3}}}$ contributes to this small discrepancy. Also, the aligned placement of the turbines leads to strong wake deficits. In FLOWERS, we neglect local wind conditions when computing the wake velocity and instead normalize the velocity deficit by the global free-stream velocity and assume a uniform C_P and C_T across all turbines; the Jensen numerical integration implementation, on the other hand, normalizes the velocity deficit with and computes C_P and C_T as a function of the local wind velocity. When the magnitude of the wake velocity deficit is pronounced, these assumptions in the FLOWERS formulation become more noticeable. Despite the differences in modeling assumptions, the AEP estimates match closely.

For a more realistic wind rose, the discrepancies in AEP between FLOWERS and the Jensen numerical integration approach are more substantial. We consider a larger wind farm of 60 turbines with 6 D spacing in the streamwise and spanwise directions and an offshore wind rose sampled from the WIND Toolkit (Draxl et al., 2015) with higher average free-stream wind speeds than those considered in the previous example. Figure 4 displays the resulting flow fields from the FLOWERS and Jensen integration models. The AEP computed with FLOWERS is about 17 % lower than that from the Jensen model. This greater difference between the two approaches in this case can be attributed to the assumption that C_T and C_P do not depend on the local flow velocity at each turbine, which breaks down for larger wind farms with more heterogeneity in the flow velocity. Including the local flow velocity in the FLOWERS formulation is part of future work.

One important feature of the FLOWERS solution is its smoothness. Despite using the discrete top-hat model, the flow field in Fig. 3b is continuous. The Fourier transform of the wind rose information and the analytical integration results in a smooth solution. On the other hand, the numerical integration in the conventional approach relies on discrete wind direction bins. The contour plot in Fig. 3c clearly illustrates the discrete boundaries of the wakes for the three dominant wind directions.

3.2 Generalized case

We now examine the differences in AEP between FLOWERS and the Jensen numerical integration more broadly. A total of 40 randomized test cases were generated. A random number of turbines between 4 and 50 was chosen for each. The layout of the turbines is randomized within a square boundary of side length 25 D, and a minimum separation of 4 D between the turbines is enforced. Each wind rose is randomly selected from the WIND Toolkit. Figure 5 displays the number of turbines and annual average wind speed for each case.

Figure 5 A total of 40 test cases with random layouts and random wind roses to compare AEP more generally.

We compare the computation time and percent difference in AEP between the two methods in Fig. 6. The FLOWERS AEP computation is on average about 18 times faster than the Jensen numerical integration. This difference scales with the size of the wind farm, in part because FLOWERS is only computing the velocity through each turbine at a single point instead of an array of points on the rotor area (25 points per rotor in the conventional wake model).

Figure 6 Comparison of computational cost and relative difference in AEP between FLOWERS and Jensen for the randomized cases.

We should note that the implementation of the Jensen numerical integration in FLORIS Version 2 is non-vectorized (i.e., calculations are mostly performed through for-loops instead of vector operations). Comparisons in computation time with a vectorized code such as PyWake would likely show a smaller discrepancy (Pedersen et al., 2019). However, the code implementation of FLOWERS is not fully vectorized either, so it is a fair comparison of FLOWERS to the Jensen numerical integration to showcase the performance of our newly developed model.

The discrepancy in AEP between the two methods is more pronounced in these randomized cases. The free-stream wind speed is not held constant here as it was in the first example. More variations in U across different wind directions and across the wind farm result in more error between the two methods due to the approximations built into the FLOWERS formulation. The common characteristics of the cases with percent difference greater than 20 % are a free-stream wind speed consistently less than 5 m s⁻¹ and a small number of predominant wind directions. The average difference between FLOWERS and the Jensen integration is 13 %.

This difference in AEP between the two methods is not necessarily a fatal flaw. FLOWERS is likely not a reliable prediction of AEP for a wind farm, but it is difficult to expect a highly accurate and precise estimate of AEP from a low-fidelity wake model anyway. However, as we will illustrate shortly, it is still possible to use the FLOWERS AEP in the layout optimization problem. In fact, the FLOWERS AEP calculation is better suited for layout optimization problems than the numerical integration method. More precise AEP estimates can always be generated as a final post-processing step after the layout optimization is complete.

3.3 Improving computational efficiency

Before addressing the utility of the FLOWERS formulation in the layout optimization problem, we can explore how to further improve the computational time. For the Jensen numerical integration method, a common approach to reduce the cost of calculating AEP is to reduce the number of wind speed–direction bins, thereby reducing the number of simulations that must be run. Figure 7a illustrates how the computational time of the AEP calculation is roughly proportional to the number of bins. Each bin adds an identical set of function calls and operations to the summing process, so the cost scales linearly. The data presented here are for a seven-by-three grid of turbines with 5 D spacing in all directions and five wind roses randomly sampled from the WIND Toolkit.

Figure 7 The effect of the resolution of wind direction bins in the Jensen model (a) and number of Fourier modes in FLOWERS (b) on cost (top panels) and accuracy (bottom panels) of the AEP calculations. Each colored line represents one of five sampled wind roses. AEP is computed for a seven-by-three grid of turbines with a spacing of 5 D. Values are normalized in each plot by the highest-resolution data point. For context, computation times are on the order of seconds.

The trade-off of sparse sampling of the wind rose is that the AEP from numerical integration is highly sensitive to the number of bins chosen. The AEP varies by as much as 40 % as we reduce the number of wind direction bins from 72 to 9. To reduce the computational cost by a factor of 2, AEP fluctuates by about 2 %; to reduce the cost by a factor of 5, AEP changes up to 10 %. The sensitivity manifests as both overestimates and underestimates of AEP, so it is not possible to assume a conservative underestimate of AEP, for example.

The equivalent idea in the FLOWERS formulation is to reduce the number of Fourier series modes. Each term in the discrete Fourier series is a single arithmetic expression, so the cost should also scale linearly with the number of terms. Figure 7b illustrates that the computational time for FLOWERS is proportional to the number of Fourier modes included in the solution. The lowest-frequency modes are used. In contrast to the Jensen numerical integration method, the AEP for FLOWERS is less sensitive to the number of Fourier modes. For all five wind roses, AEP remains virtually unchanged when using half of the maximum number of terms, and within 1 % when using only five terms. The Fourier transform is well-suited to approximate the wind speed and frequency data for each discrete wind direction and is a common tool to represent complex signals in a compact format. Reducing the number of wind direction bins is less successful because the decrease in resolution is indiscriminate: data are aggregated in uniform bins and averaged without considering the important features of the signal. Another useful feature of the Fourier expansion is that reducing the number of modes does not change the smoothness of the superimposed signal. A single Fourier mode is still a continuous sinusoidal function. On the other hand, reducing the number of wind direction bins directly causes a more discrete numerical integral. The more coarse the wind direction, speed, and frequency data are, the more sensitive the AEP calculation becomes. This has significant implications for the quality and robustness of layout optimization solutions.

With better understanding of the low sensitivity of AEP accuracy to the number of Fourier terms in FLOWERS, we have the opportunity to further reduce the computational cost. To hone in on the appropriate number of Fourier terms to use, we return to the 40 randomized test cases from Sect. 3.2. Table 1 compares the percent difference in AEP and the ratio of computation time between the Jensen numerical integration and FLOWERS with the maximum number of Fourier terms (N=37) and a truncated Fourier series (N=5). The mean and standard deviation of the AEP difference between FLOWERS and Jensen remain virtually unchanged. Only the FLOWERS AEP is computed differently between these two cases, which indicates that the FLOWERS AEP is insensitive to the number of Fourier modes included in the solution. The advantage is that the FLOWERS solution with five Fourier terms is about 110 times faster to compute than the Jensen numerical integration.

Table 1 Influence of number of Fourier modes in the FLOWERS solution for the 40 randomized AEP test cases.

We therefore recommend using this truncated FLOWERS solution. By using only one-eighth of the Fourier terms, there is a reduction in cost of roughly a factor of 6 with virtually no trade-off in accuracy. There is no reason to use the extended Fourier series if it increases the computational cost of the FLOWERS solution with no benefit to accuracy.

Figure 8 The wind rose used for the layout optimization studies performed in Sect. 4. The resolution of the wind direction bins varies from 1^∘ (a) to 5^∘ (b) to 40^∘ (c). The wind speed is a constant 8 m s⁻¹.

4 Optimization comparison

Consider nine turbines placed within a square boundary of side length 12 D. The wind is coming from the left, with a fixed speed of 8 m s⁻¹. We compare two optimizations with different objective functions: one with the FLOWERS AEP and the other with the Jensen numerical integration AEP. All other inputs and parameters in the optimization are identical: initial layout of the wind plant; wind direction, speed, and frequency distribution; wake expansion coefficient; and convergence threshold. Figure 8 shows this wind rose with three different resolutions: 1 , 5, and 40^∘ wind direction bins. These different resolutions are used in the following studies.

The AEP that drives the gradient-based optimization is different between both optimizers. However, we wish to compare the quality of the optimal solutions for both without confounding the differences in AEP discussed in Sect. 3. In this section, we use the numerical integration method to compute the AEP produced by the initial and final layouts of both optimizers for a straightforward comparison. The objective function of the FLOWERS optimizer still uses the FLOWERS AEP formulated in Sect. 2.2, but the AEP that we are reporting here is computed in the same fashion as the Jensen numerical integration AEP. The metric to compare the quality of optimal solutions is AEP gain:

$$\begin{array}{}\text{(16)}& G_{\mathrm{AEP}}={\displaystyle \frac{{\mathrm{AEP}}_{\mathrm{opt}}-{\mathrm{AEP}}_{\mathrm{init}}}{{\mathrm{AEP}}_{\mathrm{init}}}}*\mathrm{100}\mathit{\%},\end{array}$$

where AEP_init is the AEP of the initial layout and AEP_opt is the AEP of the final solution.

4.1 FLOWERS and Jensen

The first comparison is against the Jensen numerical integration model. The 5^∘ wind direction bins (i.e., B=72) are used for the Jensen model (see Fig. 8b), and N=5 Fourier terms are used for FLOWERS.

Figure 9 shows the results of the optimization for a randomized initial layout of the wind plant. The pattern of turbine placement is qualitatively different between the FLOWERS and Jensen optimizers. The Jensen optimization focuses on maximizing the streamwise spacing of the turbines by placing all turbines on either the leading or lagging edge of the wind farm's boundary. One reason why the Jensen optimizer favors this type of solution is that wake-added turbulence is included in the modeling framework, so maximizing the streamwise spacing of the turbines improves the wake recovery for downstream turbines. We have neglected wake-added turbulence in the FLOWERS formulation to present the simplest form of this new model, but we compare against the Jensen numerical integration that includes wake-added turbulence because we expect this effect to be modeled in most wake modeling codes. On the other hand, the FLOWERS solution focuses on the spanwise spacing of the turbines, placing them such that eight turbines are unaligned with respect to the predominant wind direction from the left. As expected, the FLOWERS optimization is performed faster than Jensen by a time factor of about 26 (851 s versus 31.9 s). The FLOWERS solution actually achieves a more optimal AEP gain of 14.3 % versus 12.9 % for Jensen.

Figure 9 The optimal layouts for the FLOWERS and Jensen optimizers (B=72). The FLOWERS solution required 31.9 s and achieved an AEP gain of 14.3 %. The Jensen solution required 851 s and achieved 12.9 % AEP gain. The black dots denote the initial layout.

Figure 10 A comparison of cost and performance for 10 cases with randomized initial conditions. Five Fourier modes are used for the FLOWERS solution, and 72 wind direction bins are used for the Jensen numerical integration. FLOWERS is on average about 48 times faster than Jensen and achieves an AEP gain that is 1.5 % higher.

To investigate this result more generally, we consider nine additional multistart cases with randomized initial conditions (10 in total). Figure 10 displays the computation time and AEP gain for these 10 cases. The previous example from Fig. 9 corresponds to Case 4 here. On average, FLOWERS is about 48 times faster than the Jensen optimizer. Also, FLOWERS achieves an AEP gain that is on average 1.5 % higher than Jensen.

The superior performance of FLOWERS compared to the Jensen optimizer connects back to the smooth nature of the formulation. The FLOWERS optimization space is smooth and continuous because of the Fourier transform and analytical integration. On the other hand, the Jensen optimization space is coarse because of the discrete model and numerical integration. The gradient-based optimizer thrives in the smoother optimization space of FLOWERS. More refined adjustments of the turbine positions are possible, and the optimizer is less likely to become stuck in local optimal solutions in the smooth landscape. In the discrete space of the Jensen optimization, it is more difficult for the optimizer to explore the optimization space with equivalent precision and efficiency.

To test the effect of wind rose resolution on the optimization performance, we use B=360 wind direction bins for the Jensen numerical integration, as seen in Fig. 8a; we maintain N=5 Fourier terms for FLOWERS. Figure 11 displays the computation time and AEP gain for these 10 new cases. FLOWERS is now about 680 times faster than the Jensen optimization on average; the relative improvement in computation time is due to the Jensen numerical integration AEP calculation covering 5 times more wind direction bins. The average AEP gain in FLOWERS is about 1.7 % higher than that for Jensen. The particular AEP gains that are achieved by each optimizer in this limited sample size of 10 cases are sensitive to the initial layout of the wind farm. Regardless, the overall result is that there is a negligible change in the quality of the solution that the Jensen optimizer achieves with a higher-resolution wind rose.

Figure 11 A multistart study, now for the Jensen model with B=360 wind direction bins. The Jensen optimization now takes about 680 times longer than FLOWERS and achieves an AEP gain that is 1.7 % lower.

Figure 12 The optimal layouts for the FLOWERS and Gauss optimizers (with B=360 wind direction bins). The FLOWERS solution required 50.2 s and achieved an AEP gain of 7.3 %. The Gauss solution required 10 800 s and achieved 7.4 % AEP gain. The black dots denote the initial layout.

We should note that the 10 randomized initial layouts tested here represent a limited sample size for a multistart study. The results still allow for a comprehensive discussion of the differences between the FLOWERS and Jensen models. Future work could expand on the scope of the multistart experiments to investigate whether FLOWERS converges to more consistent solutions than Jensen and also whether FLOWERS converges to a solution in fewer iterations than Jensen.

4.2 FLOWERS and Gauss

Comparing FLOWERS to a smoother wake model, in this case the Gaussian wake model (Bastankhah and Porté-Agel, 2014),

$$\begin{array}{}\text{(17)}& {\displaystyle \frac{U(x,y)}{U_{\mathrm{\infty }}}}=\left(\mathrm{1}-\sqrt{\mathrm{1}-{\displaystyle \frac{C_{\mathrm{T}}}{\mathrm{8}(kx+\mathit{ϵ}{)}^{\mathrm{2}}}}}\right)\exp \left(-{\displaystyle \frac{y^{\mathrm{2}}}{\mathrm{2}(kx+\mathit{ϵ}{)}^{\mathrm{2}}}}\right),\end{array}$$

further highlights the characteristics of the FLOWERS method relative to the numerical integration approach. Instead of discrete boundaries of the wake in the top-hat model, the Gaussian model is smooth in space. However, as a trade-off for the improved detail of the Gaussian profile, the cost of a function evaluation is higher. We use a sum-of-squares wake deficit superposition for the Gaussian wakes. Every other parameter of the optimization study is the same as the previous experiment, including 1^∘ wind direction bins (i.e., B=360).

Figure 13 Multistart experiment for FLOWERS against the Gauss model with B=360 wind direction bins. FLOWERS is about 470 times faster than Gauss on average, and the AEP gain between the two optimizers is within 0.5 % on average.

The results of one of the optimization studies are shown in Fig. 12. The Gauss optimal solution is smooth – without the discrete boundaries of the Jensen model – and the layouts are qualitatively more similar to FLOWERS with most turbines placed in angled rows. The similarity of these solutions suggests that the optimization spaces are similar despite using different wake models and AEP calculations. The similarity is not merely qualitative: the AEP gain for FLOWERS is 7.3 % and Gauss is 7.4 %. The fact that the FLOWERS gain is 0.1 % lower than that for Gauss is not necessarily a sign that the FLOWERS optimization is subpar but rather a byproduct of the many local optima. There are many possible layouts that satisfy the spacing constraints and achieve similar AEP performance, so the solutions given here likely represent local optima, not a global optimum. Since the two optimizations evaluate different objective functions and operate in different solution landscapes, they cannot be expected to arrive at identical solutions. Finding a global optimum in this unrealistic test case would require a stricter convergence criterion. We note that a more realistic wind rose would likely not enable so many local optima in the solution space, so this artifact should not be as pronounced in practice.

The trade-off for improved performance by using the Gauss numerical integration is in computational cost. The Gaussian optimization took 10 800 s (3 h), while the FLOWERS optimization only required 50.2 s. This is an improvement by a factor of 216 for FLOWERS.

Figure 13 displays the results for the 10 multistart cases. Figure 12 is Case 1 in these plots. On average, the Gauss optimization takes 470 times longer than FLOWERS. The average AEP gain in FLOWERS is 0.5 % higher than Gauss. The takeaway from this small difference in AEP gain is that FLOWERS and Gauss produce comparable solutions, not that FLOWERS is better at finding a more optimal solution. A stricter convergence criterion would force the optimizers to search longer for the global solution, which would likely cause these AEP gains to grow even more similar.

This experiment suggests that the smoothness of the Gaussian model compared with the Jensen model is the most likely explanation for the improved performance. While the number of wind direction bins is unchanged, the flow field for each simulation is smoother with the Gaussian model. When a turbine's position is adjusted, there is no binary switch between being within the wake or outside of it; this discrete change in wake velocity would cause the AEP to be sensitive to slight perturbations in the turbine positions. On the other hand, in the Gaussian model, a small change in turbine position results in a similarly small change in the wake deficit because of the smooth profile. This continuity produces a smoother solution space, which enables the optimizer to move along more subtle gradients and achieve more optimal solutions than the Jensen optimizer.

While the quality of solutions between FLOWERS and the Gauss optimizer is comparable, there is no contest in terms of cost. The FLOWERS optimization is 2 orders of magnitude faster than the Gauss optimizer and can produce optimal solutions with equivalent performance. Moreover, we are only comparing the results for a wind plant with nine turbines. As illustrated in Fig. 6a, the computational cost scales with the number of turbines more sharply with the numerical integration approach. It is possible that this factor of 470 could grow to a factor of 1000 or more for a larger wind farm.

We have demonstrated that the FLOWERS AEP is insensitive to the number of Fourier series terms and have used the truncated series to achieve similar performance to the Gauss optimization. We also previously showed that the AEP calculated from numerical integration is extremely sensitive to the resolution of the wind direction bins. For a fair and comprehensive comparison, the Gauss optimization should be performed with a limited number of wind direction bins to mimic the reduction in cost that was implemented in FLOWERS. To match N=5 Fourier terms, B=9 wind direction bins are now used in the Gaussian optimization, as seen in Fig. 8c. The results in this case are shown in Fig. 14. The reduction in wind rose resolution brings the cost of the Gaussian optimization down significantly such that the FLOWERS optimizer is only about 5 times faster on average. However, the AEP gain of the optimal solutions for this Gaussian optimizer with a coarse wind rose is poor. The AEP gain for FLOWERS is 7.5 % higher than Gauss on average, which is an improvement of a factor of 2.7. This experiment proves that the Gaussian optimizer cannot achieve greater computational efficiency by manipulating the resolution of the wind rose without substantially impacting the quality of its solutions.

Figure 14 The 10 randomized cases for FLOWERS against the Gauss model with B=9 wind direction bins. FLOWERS is only about 5 times faster than Gauss now, on average. The average AEP gain for FLOWERS is about 12 %, but only about 4.5 % for Gauss.

5 Potential model improvements

As we have demonstrated, FLOWERS is able to match the performance of conventional layout optimization methods despite simplifications in its formulation. The discrepancies in AEP estimates between FLOWERS and the Jensen and Gauss numerical integrations did not inhibit its application to the optimization problem. However, we could enhance the accuracy of FLOWERS by improving the following.

• Power integral. We introduced a simplification to make the integration of turbine power tractable by computing AEP as a function of ${\stackrel{\mathrm{‾}}{U}}^{\mathrm{3}}$ rather than $\stackrel{\mathrm{‾}}{U^{\mathrm{3}}}$. We would reexamine the integration of Eq. (13) to avoid this simplification, which introduces errors when wind speed varies across different wind directions.

• Coefficient of power. C_P is currently defined as a function of the average wake velocity, making it a constant. We aim to incorporate C_P as a function of the wake velocity for each wind direction such that average power is computed exactly: $\stackrel{\mathrm{‾}}{P\left(U\right)}=\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }A\stackrel{\mathrm{‾}}{C_{\mathrm{P}}\left(U\right)U^{\mathrm{3}}}$.

• Local flow conditions. We can define the wake velocity deficit relative to the local flow velocity rather than the free stream, which will better capture the influence of upstream turbines and development of the flow as it moves through the wind farm. This improvement will require an iterative approach and is particularly expected to improve results in aligned cases such as the ones discussed in Sect. 3.1.

• Gaussian model. We currently integrate a classical top-hat wake deficit model. We expect that the AEP estimates would be more accurate by integrating a Gaussian wake model instead.

6 Conclusions

The objective of this paper was to develop a novel analytical formulation of annually averaged wake velocity to use in a layout optimization problem and demonstrate its effect on reducing the computational cost of these studies. We derived the equations for the analytical integration of the top-hat wake deficit model. The wind speed and wind direction frequency distributions were expressed as a Fourier series to facilitate the integration.

The annually averaged wake velocity was used to compute AEP. We approximated the average power by using the average wake speed cubed rather than an average of the cube of the wake speed, which introduces error when there are pronounced wakes or the wind speed varies significantly across different wind directions. Also, the local wind speed's effect on turbine power production and thrust was not accounted for. These simplifications introduced error that led to the AEP computed in FLOWERS differing from the Jensen numerical integration approach by about 13 %.

Fortunately, these limitations in the accuracy of the FLOWERS AEP do not preclude its use in the optimization problem. The FLOWERS optimizer built around the Jensen wake model finds optimal wind plant layouts with AEP comparable to an optimizer that numerically integrates a Gaussian wake model. This finding is unexpected but promising because it implies that the mathematical formulation behind FLOWERS compensates for a more simplistic wake model to achieve similar results to a more sophisticated wake model. The clear advantage of FLOWERS, then, is the robust layout optimization performance while achieving a reduction in computational cost of 2 orders of magnitude. We believe that this improvement in computation time will scale better with wind farms containing more than the nine turbines studied here.

This achievement could translate to the difference between running an optimization study in 10 min versus 5 d, or between running the study on a personal laptop versus a high-performance computer cluster. This technique could open the door for other areas of research in layout optimization, including optimization under uncertainty, by making these studies more accessible and less costly. Moreover, the new conceptualization of the wake velocity deficit could inspire brand new areas of research in wake modeling and wind plant control and optimization.

This paper serves as a foundation for future work on the FLOWERS formulation. Since the motivation of this approach was to improve computational cost, one avenue to explore is further optimization of the FLOWERS code. Wind plant layout and yaw steering co-design is a popular area of research, and another potential application for FLOWERS if yaw deflection models could be included in the formulation. Future work will also focus on studying the effects of superposition methods and model uncertainty in the FLOWERS formulation. We also plan to validate the performance of the FLOWERS optimal solutions with high-fidelity simulations.