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 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. We propose a novel formulation 5 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 about 700 times faster. The analytical integral and the use of a Fourier expansion to express the wind speed and wind direction frequency create a more smooth solution space for the gradient-based optimizer to excel 10 compared with the discrete nature of the existing weighted-averaging power calculation.


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(Barthelmie et al., , 2009. 15 When turbines are placed within about 3 rotor diameters, these power losses can be as high as 40% ; 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 tophat 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 tophat model with a 20 continuous profile, such as the Jensen-cosine (Tian et al., 2015(Tian et al., , 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 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 45 the wind rose contains only one non-zero 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 Figure 1b-c. For example, in Figure 1b the frequency of the 180 • wind direction is greater than that for 225 • (and the freestream 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 50 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.
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 coming from every direction. Since the 55 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 Section 2, including its application in the wind plant layout optimization problem. In Section 3, the AEP calculations from FLOWERS are 60 compared to the conventional numerical integral method, which is the standard method in the NREL FLORIS model (NREL, 2021). Finally, in Section 4, the performance of the FLOWERS AEP in the optimization of wind farm layouts is compared against the traditional optimizer. To derive a mathematical formulation for the time-averaged flow distribution, we use the classical Jensen (tophat) wake deficit model (Jensen, 1983): where x and y are the streamwise and spanwise position, respectively, normalized by the rotor radius. 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 70 assumed to be constant for simplicity, k is the wake expansion coefficient, U is the wake speed and U ∞ is the freestream wind speed. W (x, y) represents the Jensen wake region: W (x, y) = 1 if |y| ≤ kx + 1 and is zero elsewhere.
We transform Eq. 1 from Cartesian to polar coordinates denoted by r and θ, where r = x 2 + y 2 and tan(θ) = y/x. We allow the wind direction θ to be variable in Eq. 2:

75
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 U (r, θ) is written as: We define two new variables to simplify this expression. Let u = θ − θ be the angular position relative to the wind direction.
Also, the wake region W (r, θ, θ ) is zero for θ outside of the wake region defined by θ c : θ = θ ± θ c . In our coordinate system, the wake geometry is defined by the line sin(θ c ) = k cos(θ c ) + 1/r. We solve this equation for r > 1 since we are interested in the wake speed at positions outside of the rotor area: (4) 85 We must include the interaction between multiple wakes for this formulation to be useful in systems of several turbines. The wake velocity is defined as the difference between the freestream velocity and the velocity deficit. Eq. 3 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:

90
where r i and θ i are the relative radius and polar angle with respect to each turbine position (r i = (x − x i ) 2 + (y − y i ) 2 and tan(θ i ) = 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. The product U ∞ (θ )f (θ ) is a vector with length equal to the number of wind direction bins. We define g(θ ) = 1 2π U ∞ (θ )f (θ ). If we expand g(θ ) with a Fourier series: the wind rose is defined continuously rather than discretely, where a 0 , a n and b n are Fourier coefficients and can be easily found for a given g(θ ). For a wind rose with B wind direction bins, the maximum number of terms in this discrete Fourier 100 transform is N = ceiling(B/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. 5 using a Taylor series: The first integral in the right-hand side of Eq. 7 represents the weighted-average of the freestream velocity, denoted by U ∞ hereafter. Using Eq. 6, The second integral on the right-hand side of Eq. 7 represents the average wake velocity deficit ∆U i (r i , θ i ). Using Eq. 6, a n cos(n(θ i − u)) + b n sin(n(θ i − u)) 1 (kr i + 1) 2 + kr i u 2 (kr i + 1) 3 du.
Solving the above integral yields: The time-averaged wake speed U (r, θ) at a given location (r, θ) is the difference between these two terms: where U ∞ and ∆U i (r i , θ i ) can be found from Eq. (8) and Eq. (9), respectively.

Annual Energy Production
The power P produced by a turbine is a function of the incoming wind speed: 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 velocity at the location of the center of the turbine when it is not present. The average turbine power would be 120 P (U ) = 1 2 ρAC p (U )U 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 3 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 125 average wake speed U , and we substitute U 3 into Eq. 11 instead of evaluating U 3 . Therefore, the AEP for a given turbine is

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 2M design variables to independently specify the position of each turbine in 2D space.

130
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 two 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 : 135 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 140 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, 145 just as it will in the conventional 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 standard technique. A more accurate AEP estimate can be added as a final post-processing step once the layout optimization is complete. The FLOWERS AEP is 1.8% 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 U 3 and U 3 contributes to this small discrepancy. Also, the aligned 160 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 freestream velocity and assume a constant C p across all turbines. 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.
One important feature of the FLOWERS solution is its smoothness. Despite using the discrete tophat model, the flow field in

Generalized Case
We now examine the differences in AEP between FLOWERS and the numerical integration approach more broadly. Forty 170 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 25D, and a minimum separation of 4D between the turbines is enforced. Each wind rose is randomly selected from the WIND Toolkit (Draxl et al., 2015). Figure 3 displays the number of turbines and annual average wind speed for each case.
We compare the computation time and percent difference in AEP between the two methods in    The discrepancy in AEP between the two methods is more pronounced in these randomized cases. The freestream wind speed is not held constant here as it was in the first example. More variations in U across different wind directions results in more 180 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 freestream 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 8%.
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 185 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 conventional method. More precise AEP estimates can always be generated as a final post-processing step after the layout optimization is complete.

190
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 conventional 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 5a 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 195 presented here is for a 7x3 grid of turbines with 5D spacing in all directions and five wind roses randomly sampled from the WIND Toolkit.
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 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 5b illustrates that the computational time for FLOWERS is proportional to the number of Fourier modes included in the solution.

205
The lowest-frequency modes are used. In contrast to the 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 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 is aggregated in uniform bins 210 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 almost 150 times faster to compute than the numerically-integrated Jensen method.
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 225 series if it only increases the computational cost of the FLOWERS solution.

Optimization Comparison
Consider 9 turbines placed within a square boundary of side length 12D. 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 conventional AEP via numerical integration. All other inputs and parameters in the optimization are identical: initial layout 230 of the wind plant, wind direction, speed, and frequency distribution, wake expansion coefficient, and convergence threshold. 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 Section 3. In this section, 235 we use the numerical integration method to compute the AEP produced by the initial and final layouts of both optimizers for a

FLOWERS and Jensen
The first comparison is against the integrated Jensen model. 5   To investigate this result more generally, we consider nine additional multistart cases with randomized initial conditions (ten in total). Figure 8 displays the computation time and AEP gain for these ten cases. The previous example from Figure   7 corresponds to Case 3 here. On average, FLOWERS is about 59 times faster than the Jensen optimizer. Also, FLOWERS achieves an AEP gain that is on average 1% higher than Jensen. 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 260 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 model, as seen in Figure 6a); we maintain N = 5 Fourier terms for FLOWERS. Figure 9 displays the computation time and AEP gain for these ten new cases. FLOWERS is now about 200 times faster than the Jensen optimization on average; the relative improvement in computation time is due to the Jensen AEP calculation covering five times more wind direction bins.

265
The average AEP gain in FLOWERS is about 1.8% 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.
We should note that the 10 randomized initial layouts tested here represent a limited sample size for a multistart study.

270
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.

FLOWERS and Gauss
Comparing FLOWERS to a smoother wake model, in this case the Gaussian wake model (Bastankhah and Porté-Agel, 2014): further highlights the characteristics of the FLOWERS method relative to the conventional numerical integration approach.
Instead of discrete boundaries of the wake in the tophat 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, The results of one of the optimization studies is shown in Figure 10. 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.9% and Gauss is about 7.3%. The fact that the Gaussian gain is 0.6% lower than that for FLOWERS is not necessarily a sign that the Gauss 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. However, the trade-off for improved performance by using the Gauss 290 wake model with numerical integration is in computational cost. The Gaussian optimization took 16,700 s (4.6 hr), while the FLOWERS optimization only required 46.4 s. This is an improvement by a factor of 360 for FLOWERS. Figure 11 displays the results for the 10 multistart cases. Figure 10 is Case 1 in these plots. On average, the Gaussian optimization takes 690 times longer than FLOWERS. The average AEP gain in FLOWERS is 0.1% higher than Gauss. The takeaway from this small difference in AEP gain is that FLOWERS and Gauss produce comparable solutions, not that FLOW-

295
ERS 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. a similarly small change in the wake deficit because of the smooth profile. This continuity produces a more smooth solution space, which enables the optimizer to move along more subtle gradients and achieve more optimal solutions than the Jensen optimizer.

305
While the quality of solutions between FLOWERS and the Gauss optimizer are comparable, there is no contest in terms of cost. The FLOWERS optimization is almost three 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 Figure 3a, the computational cost scales with the number of turbines more sharply with the numerical integration approach. It is possible that this factor of 690 could grow to a factor of 1,000 or more for a larger wind farm. 310 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 315 Gaussian optimization, as seen in Figure 6c. The results in this case are shown in Figure 12. The reduction in wind rose resolution brings the cost of the Gaussian optimization down significantly such that the FLOWERS optimizer is only about 6 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% higher than Gauss on average, which is an improvement of a factor of 2.5. This experiment proves that the Gaussian optimizer cannot achieve greater computational efficiency by manipulating the resolution 320 of the wind rose without substantially impacting the quality of its solutions.

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 integration of the Jensen and Gaussian models did not inhibit its application to the optimization problem. However, we could enhance the accuracy of 325 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 U 3 rather than U 3 . We would reexamine the integration of Eq. 11 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 330 incorporate C p as a function of the wake velocity for each wind direction such that average power is computed exactly: -Local flow conditions: We can define the wake velocity deficit relative to the local flow velocity rather than the freestream, which will better capture the influence of upstream turbines and development of the flow as it moves through the wind farm. This improvement is particularly expected to improve results in aligned cases such as the one discussed 335 in Section 3.1.
-Gaussian model: We currently integrate a classical tophat wake deficit model. We expect that the AEP estimates would be more accurate by integrating a Gaussian wake model instead.

Conclusions
The objective of this paper was to develop a novel analytical formulation of annually-averaged wake velocity to use in a layout 340 optimization problem and demonstrate its effect on reducing the computational cost of these studies. We derived the equations for the analytical integration of the tophat 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 345 wakes or the wind speed varies significantly across different wake directions. Also, the local wind speed's effect on power production was not accounted for, and a constant power coefficient was assumed. These simplifications introduced error that led to the AEP computed in FLORIS differing from the numerical integration approach by about 8%.
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 350 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 as 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 over two 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.

355
This achievement could translate to the difference between running an optimization study in 10 minutes versus 5 days, 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.

360
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. We also plan to validate the performance of the FLOWERS optimization against high-fidelity simulations. Competing interests. The authors declare that they have no conflict of interest.