The revised FLORIDyn model: implementation of heterogeneous flow and the Gaussian wake

In this paper, a new version of the FLOw Redirection and Induction Dynamics (FLORIDyn) model is presented. The new model uses the three-dimensional parametric Gaussian FLORIS model and can provide dynamic wind farm simulations at a low computational cost under heterogeneous and changing wind conditions.

Both FLORIS and FLORIDyn are parametric models which can be used to simulate wind farms, evaluate controller performance and can serve as a control-oriented model. One central element in which they differ is in their representation of flow dynamics: FLORIS neglects these and provides a computationally very cheap approximation of the mean wind farm flow. FLORIDyn defines a framework which utilizes this low computational cost of FLORIS to simulate basic wake dynamics. This is achieved by creating so-called observation points (OPs) at each time step at the rotor plane which inherit the turbine state.

In this work, we develop the initial FLORIDyn framework further considering multiple aspects. The underlying FLORIS wake model is replaced by a Gaussian wake model. The distribution and characteristics of the OPs are adapted to account for the new parametric model but also to take complex flow conditions into account. To achieve this, a mathematical approach is developed to combine the parametric model and the changing, heterogeneous world conditions and link them with each OP. We also present a computationally lightweight wind field model to allow for a simulation environment in which heterogeneous flow conditions are possible.

FLORIDyn is compared to Simulator for Offshore Wind Farm Applications (SOWFA) simulations in three- and nine-turbine cases under static and changing environmental conditions. The results show a good agreement with the timing of the impact of upstream state changes on downstream turbines. They also show a good agreement in terms of how wakes are displaced by wind direction changes and when the resulting velocity deficit is experienced by downstream turbines. A good fit of the mean generated power is ensured by the underlying FLORIS model. In the three-turbine case, FLORIDyn simulates 4 s simulation time in 24.49 ms computational time. The resulting new FLORIDyn model proves to be a computationally attractive and capable tool for model-based dynamic wind farm control.

1 Introduction

In recent years, the topic of wind farm control has gained traction as renewable energies become more and more relevant for the current and future energy mix. Maximizing the power generated by a wind farm is not a trivial task as the turbine-to-turbine interaction is characterized by the complex flow, large delay times and an ever-changing environment. In order to describe the wind field, parametric steady-state approximations have been developed. These describe the mean behavior of the flow with parametrized analytical expressions rather than differential equations. A first approach was presented by Jensen (1983), which motivated years later the development of more refined steady-state models, such as the Zone FLORIS model (Gebraad et al., 2014). With these low-computational-cost and easy-to-implement wake descriptions, it is possible to develop a model-based control algorithm. These control strategies have managed to improve the power generated in high-fidelity simulations e.g., Gebraad et al. (2014) and in field experiments (Fleming et al., 2017). The success of parametric steady-state models opens up the question of whether it is possible to overcome one of their great shortcomings: the lack of dynamics. A low-computational-cost dynamic wake description can be used to more accurately describe the wake behavior on smaller timescales, during turbine state changes and during environmental changes. This could lead to more sophisticated control approaches and wind farm analysis methods.

There have been efforts to implement parametric models in a dynamic manner, some of which are described here. For a more in-depth discussion of the current state of the art, the interested reader is referred to the review by Kheirabadi and Nagamune (2019) and more recently, Andersson et al. (2021). In the current literature, we have identified two major trails of publications, which will be briefly discussed below.

The first research trail begins with the Aeolus SimWindFarm toolbox (Grunnet et al., 2010), which is publicly available. The toolbox uses the Jensen model (Jensen, 1983), coupled with a dynamic description of the centerline and a wind field grid. The centerline would imitate the wake meandering effect based on passive tracers, traveling with the synthetically generated turbulent wind speed. A number of limitations have been imposed for this toolbox: the mean wind speed and direction are constant, the flow field is calculated in 2D, and the turbines operate with fixed yaw angles. The toolbox has enabled the work of Poushpas and Leithead (2014), who used the Frandsen multiple wake model (Frandsen et al., 2006) and added a description of turbine dynamics to estimate fatigue loads. The model is then used to perform induction control based on lookup tables of the thrust and power coefficients with the goal to redistribute loads. This work later inspired the dynamic wind farm simulator, introduced in Bossanyi (2018). The model adds wake steering to the fatigue load estimation and induction control capabilities. To model the effect of yawing the turbine, the deflection formulation of Jiménez et al. (2010) is used. Based on data from the in-house code Bladed, the author formulates the effect of yaw misalignment on the power coefficient by a polynomial expression based on the blade pitch. The wind field is represented by low- and a high-frequency wind speed variations. The low-frequency variations are correlated across the wind farm and cause wake meandering and advection. The high-frequency part is uncorrelated between the turbines and is superimposed with the wake deficits. Lastly, the wake model is switched to the Ainslie model (Ainslie, 1988).

A second trail of publications can be found starting with Shapiro et al. (2017), where the authors use the previously mentioned Jensen model and extend it to incorporate the impact of time-varying extraction of kinetic energy of turbines due to induction control. Assuming a constant wind direction and wind speed, the authors derive a linear approximation of the wake advection velocity based on the laws of momentum conservation and mass conservation. The result is a one-dimensional partial differential equation to describe the dynamic wake behavior. The model neglects possible changes of the wake expansion due to a changing thrust coefficient and also does not incorporate yaw angle changes. In Shapiro et al. (2018), the authors extend their model to also take the effects of yawing into account. Most recently, this approach inspired the development of the Floating Offshore Wind Farm Simulator, published in Kheirabadi and Nagamune (2021). The authors extend the momentum conservation equations to incorporate time-varying free-stream wind velocity effects. Additionally, they couple the model to a dynamic description of floating platforms, restricted by mooring lines. The authors closely follow Bastankhah and Porté-Agel (2016) to derive a parametric Gaussian velocity shape for their model.

Alongside the two discussed trails of publications, the dynamic wake meandering (DWM) model was developed. The DWM model, first presented by Larsen et al. (2008) and later calibrated and refined by Madsen et al. (2010), proposes an approach much closer to established CFD methods. The model follows a pseudo-Lagrangian approach and creates turbulence boxes around the wake deficit which is created by the turbine. These boxes are then subject to a synthetic turbulent wind field, which allows the modeling of the wake meandering effect. The DWM model puts a focus on load estimation next to the power generated and simulates the turbine by coupling a CFD actuator disc model with an aeroelastic model. Compared to the other mentioned models, the DWM model presents a synergy of CFD methods with engineering approaches.

Another early attempt to derive a dynamic model from a parametric steady-state model was published by Gebraad and van Wingerden (2014), who utilized the just-published FLORIS model (FLOw Redirection and Induction in Steady state, Gebraad et al., 2014) and created the FLORIDyn model (FLOw Redirection and Induction Dynamics). FLORIDyn creates so-called observation points (OPs) at the rotor plane which travel downstream at hub height with the effective wind speed. Their path follows the zone boundaries described by the FLORIS model. The wake deficit and shape depend on the yaw angle and the induction factor. Changes in these variables travel with the OPs and cause a delayed effect at downstream turbines. The authors derive a state-space representation of the model behavior and validate it in a six-turbine simulation against the high-fidelity large eddy simulation environment SOWFA (National Renewable Energy Laboratory, 2020). The state-space representation is then used to implement a Kalman filter for flow field estimation (Gebraad et al., 2015). The model does have shortcomings: due to the two-dimensional flow, shear and veer effects can not be captured, the simulations only work in one wind direction and they do not capture turbulent effects. Furthermore, due to the way the OPs travel, parts of the wake can overlap and can create a faulty wake representation.

In this paper, we aim to overcome these issues and bring the FLORIDyn approach into a form where it can incorporate heterogeneous and changing flow conditions, wind shear, and added turbulence levels. To achieve these changes, we rework the framework to use a Gaussian FLORIS model (Bastankhah and Porté-Agel, 2016). This requires a new formulation of the OP behavior. Due to these changes, the wakes can also incorporate locally different and changing flow conditions, such as wind speed, direction and ambient turbulence intensity. To drive the model, a concept of a wind field model is presented as well. The framework is then compared to the simulation environment SOWFA in three- and nine-turbine cases. Furthermore, in order to allow for collaboration and extension, the code is published in its entirety (Becker, 2022a). The resulting Gaussian FLORIDyn model is a capable, open-source alternative to the few other existing in-house parametric dynamic models, developed for wind farm control purposes.

The remainder of this paper is organized as follows: Sect. 2 discusses the relevant characteristics of the former FLORIDyn framework and how it is adapted. The simulation results are presented in Sect. 3, which also discusses the computational performance. Section 4 concludes the paper and gives recommendations for future work.

2 A new parametric dynamic wind farm model

In this section, the new Gaussian FLORIDyn model is introduced. To prevent confusion, we will refer to the models of Gebraad et al. as the Zone FLORIS model (Gebraad et al., 2014) and the Zone FLORIDyn model (Gebraad and van Wingerden, 2014). The Gaussian model by Bastankhah and Porté-Agel (2016) will be referred to as Gaussian FLORIS model.

Figure 1 Sketched shape of the wake with the different sections, the deflection and areas of equal relative reduction by the Gaussian shape.

As the new Gaussian FLORIDyn model is building upon previous work, Sect. 2.1 and 2.2 briefly introduce the terminology and properties of the underlying Gaussian FLORIS model and the Zone FLORIDyn framework. The novel Gaussian FLORIDyn model makes changes to the Zone FLORIDyn framework. These are discussed in Sect. 2.3. Section 2.4 describes how heterogeneous environmental conditions are taken into account. To get the power coefficient (C_P) and the thrust coefficient (C_T) values closer to the validation platform SOWFA, a lookup table was generated (Sect. 2.5). Lastly, a basic wind field model is given in Sect. 2.6. It is built to provide the heterogeneous field conditions to evaluate the FLORIDyn model.

In the wake coordinate system, 𝒦₁, x₁ describes the downwind direction, y₁ the horizontal crosswind direction and z₁ the vertical crosswind direction (Fig. 1). In this coordinate frame, the rotor center is always located at $[\mathrm{0},\mathrm{0},\mathrm{0}{]}^{\top }$. This coordinate system is not to be confused with the longitudinal (x₀), latitudinal (y₀) and vertical (z₀) world coordinate system 𝒦₀.

2.1 The Gaussian FLORIS model

The core of the used Gaussian FLORIS model is based on the work of Bastankhah and Porté-Agel (2016). This work describes a parametric, three-dimensional wake with a Gaussian-shaped wind speed recovery. As it has been applied and described in previous publications (e.g., Farrell et al., 2021), only the basic terminology is introduced here as well as the wake shape. In the present work of this paper, the model has been extended with the calculation of added turbulence as proposed by Crespo and Hernández (1996). The power calculation has been extended by the $\cos (\mathit{\gamma }{)}^{p_{\mathrm{p}}}$ adaptation to the yaw angle (Medici, 2005) and an efficiency term η for tuning (Gebraad et al., 2014). Figure 1 depicts an illustration of the wake with its three areas: the potential core, the near-wake area and the far-wake area. For all areas a reduction factor $r=\mathrm{\Delta }u/u_{\mathrm{free}}$ can be calculated, where u_free is the free wind speed and Δu is the wind speed deficit. The potential core is a region from jets in a coflow (Lee and Chu, 2003). Here, it is used to approximate the immediate region behind the rotor plane. Within the potential core, r is constant. In the near- and far-field area r reduces to 0, following a Gaussian shape with the extremum at the centerline or border of the potential core. The recovery rate is based on σ_y and σ_z in the respective crosswind directions. The potential core width is described by w_y,pc and w_z,pc, which continuously decrease for the length of the potential core x_c. Lastly, the deflection δ returns the position of the centerline.

The mentioned variables are dependent on turbine states, such as the thrust coefficient C_T and the yaw angle γ, the ambient turbulence intensity I₀, and a set of 10 parameters. The parameters adjust wake properties such as the recovery rate, the expansion rate, the sensitivity to added turbulence levels and the influence of the yaw angle. The values of the parameters are listed in Table 1 in Sect. 3.

2.2 The Zone FLORIDyn model

An initial FLORIDyn model was published in Gebraad and van Wingerden (2014). The model is based on the previously published Zone FLORIS model, which approximates the wake shape with three zones: near field, far field and mixing zone (Gebraad et al., 2014). Every zone has a formulation of the velocity recovery in downstream direction. To introduce dynamics, observation points (OPs) are created at the rotor plane at each time step.

The OPs serve the purpose to describe the local FLORIS wake characteristics at their location. To do that, they inherit the turbine states at the time of their creation which are necessary to calculate the FLORIS wake. With time, each OP travels downstream, representing a mass of air traveling in the wind. Their travel path is determined by the borders of the FLORIS wake zones. The speed they travel with is equal to the effective wind speed they represent. Figure 2 shows the basic concept. Initial OPs are colored black to stress that they inherited the same state. The OPs created after the yaw step are colored white, showing that their inherited state differs.

Figure 2 Creation and propagation of the OPs: in panel (a) a set of OPs is created, inherits the turbine state and travels downstream, following the FLORIS wake shape, shown in panel (b). In panel (c) the turbine state changed and the new OPs inherit a different state (now colored white) and follow the new, dark-indicated wake shape in panel (d).

With this framework, the steady-state wake represents the known FLORIS wake, but other than in FLORIS, changes propagate through the wake instead of instantly affecting turbines downstream. If, for instance, the yaw angle of the turbine changes, the new generation of OPs will inherit the new angle while old OPs still travel according to the previous angle.

In the case of overlapping wakes, an OP travels into the wake of another turbine. The OP locates the closest up- and downstream OPs from the foreign wake and interpolates their reduction factor at its location. In this model, the resulting reduction of the free wind speed is calculated as follows:

$$\begin{array}{}\text{(1)}& u_{\mathrm{eff},\mathrm{OP}}\left(u_{\mathrm{free},\mathrm{OP}},r_{\mathrm{own}},r_{\mathrm{f},\mathrm{OP}}\right)=u_{\mathrm{free},\mathrm{OP}}\left(\mathrm{1}-r_{\mathrm{own}}\right)\underset{r_{\mathrm{f},\mathrm{OP}}}{\underbrace{\prod _{i=\mathrm{1}}^{n_{\mathrm{T}}}\left(\mathrm{1}-r_{\mathrm{i}}\right)}},\end{array}$$

where u_free,OP is the free wind speed at the OP's location. This wake interaction model could also be exchanged for another formulation. The wind speed reduction r_own is based on the OP's own wake, and r_i is the interpolated reduction of one of the n_T upwind turbines.

To calculate the effective wind speed at the rotor plane, the model calculates an effective velocity reduction factor r_T for every turbine at every time step. The algorithm combines the reduction of each upstream turbine by a root sum square. Within one wake, the reduction factors of the zones are summed, weighted by the overlapping area with the rotor plane.

2.3 Changes to the FLORIDyn approach

Due to the changed underlying FLORIS model, the FLORIDyn approach needs to be adapted. Specifically, the move to a three-dimensional flow field requires a fitting distribution of the OPs, which is discussed in Sect. 2.3.1. This opens up the possibility to reformulate the calculation of the effective wind speed at the rotor plane, which is presented in Sect. 2.3.2. The travel speed of the OPs is addressed in Sect. 2.3.3. In this section, we use the wake coordinate system 𝒦₁, indicated by the lower index 1 (e.g., y_1,OP). The relation between world and wake coordinate system will be explained in Sect. 2.4.

2.3.1 Distribution of the observation points

By changing the underlying FLORIS model, the travel path of the OPs and their distribution has to be rethought. The Gaussian FLORIS model does not have defined borders, and it is three-dimensional. To cover the crosswind wake area regularly for any number of OPs, an algorithm based on the sunflower distribution was used (Vogel, 1979). The algorithm returns a relative crosswind coordinate $({\mathit{\nu }}_y,{\mathit{\nu }}_z)\in [-\mathrm{0.5},\mathrm{0.5}]$ for a given number of OPs. We used 50 OPs per time step. To cover the majority of the Gaussian wake influence, the wake width was chosen to be ±3σ_y and ±3σ_z from the centerline and the potential core. The following equation is used to calculate the position of an OP in the wake coordinate system:

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }{y}_{\mathrm{1},\mathrm{OP}}\left({\mathit{\nu }}_{y,\mathrm{OP}},{\mathit{\sigma }}_y,w_{y,\mathrm{pc}},\mathit{\delta }\right)={\mathit{\nu }}_{y,\mathrm{OP}}\left(\mathrm{6}{\mathit{\sigma }}_y+w_{y,\mathrm{pc}}\right)+\mathit{\delta },\text{(3)}& {\displaystyle }& {\displaystyle }{z}_{\mathrm{1},\mathrm{OP}}\left({\mathit{\nu }}_{z,\mathrm{OP}},{\mathit{\sigma }}_z,w_{z,\mathrm{pc}}\right)={\mathit{\nu }}_{z,\mathrm{OP}}\left(\mathrm{6}{\mathit{\sigma }}_z+w_{z,\mathrm{pc}}\right).\end{array}$$

Note that this model only assumes a horizontal deflection. To add a vertical deflection, due to rotor tilt for instance, Eq. (3) needs to be adapted accordingly. For simplicity's sake σ_y is used, which represents σ_y,nw for $\mathrm{0}<x_{\mathrm{1}}\le x_{\mathrm{c}}$ and σ_y,fw for x₁>x_c. Respectively, σ_z is defined the same way. The variable δ describes the deflection of the centerline. If OPs travel below z₁=0 they are ignored. Since ν_y and ν_z are not changed during the simulation, they can be calculated a priori. They are then used in every time step for the new generation of OPs. OPs with the same relative coordinate follow each other and form what is called a chain. The number of chains is equal to the number of OPs created at each time step.

2.3.2 Wind speed at the rotor plane

Since OPs are created at the rotor plane and they interact with foreign wakes, they can be used to estimate the effective wind speed for the power generation. To do that, they have to be distributed across the rotor plane rather than the wake area:

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle }{y}_{\mathrm{1},\mathrm{OP}}\left({\mathit{\nu }}_{y,\mathrm{OP}},\mathit{\gamma }\right){\mid }_{x_{\mathrm{1}}=\mathrm{0}}={\mathit{\nu }}_{y,\mathrm{OP}}D\cos \mathit{\gamma },\text{(5)}& {\displaystyle }& {\displaystyle }{z}_{\mathrm{1},\mathrm{OP}}\left({\mathit{\nu }}_{y,\mathrm{OP}}\right){\mid }_{x_{\mathrm{1}}=\mathrm{0}}={\mathit{\nu }}_{z,\mathrm{OP}}D.\end{array}$$

The next step is to determine the area represented by every OP. This is done offline by generating a Voronoi pattern (Voronoi, 1908a, b) with the OPs' relative location as seeds and a circular boundary with radius 0.5. The area of the resulting polygons is normalized by the rotor area and used as weight. All weights are stored in the vector w.

During the simulation, the OPs calculate the reduction of foreign wakes r_f,OP on themselves as shown in Eq. (1). Stored in a vector ${\mathit{r}}_{\mathrm{f}}=[r_{\mathrm{f},\mathrm{1}},\mathrm{\cdots },r_{\mathrm{f},\mathrm{nOP}}{]}^{\top }$ the effective wind speed at the rotor plane is calculated as follows:

$$\begin{array}{}\text{(6)}& u_{\mathrm{eff}}={\mathit{w}}^{\top }\left({\mathit{r}}_{\mathrm{f}}\circ \mathit{u}\right),\end{array}$$

where ∘ stands for the element-wise multiplication and u represents a vector of the free wind speeds at the locations of the OPs. An OP considers itself influenced by a foreign wake if the closest foreign OP is less than $\frac{\mathrm{1}}{\mathrm{4}}D$ away. This is an arbitrary chosen threshold to reduce the number of OPs for the interaction interpolation. As the outer wake OPs represent the most recovered sections of the wake, this still results in a smooth influence transition.

2.3.3 Travel speed

In the former version of the FLORIDyn model, the OPs travel with the effective wind speed they represent. Regions in the center of the wake with lower effective wind speeds therefore propagate the changes slower than the outer areas. While this seems an intuitive choice, it leads to problems. Initial simulation results showed that, in comparison to the SOWFA simulation, the effects of a state change arrive noticeably slower in FLORIDyn at downstream turbines. Also, due to the difference in OP travel speed, the outer regions adapt their shape earlier in a downstream location, which leads to overlapping areas with the slow regions, which have not adapted yet. This makes the wake representation not injective anymore: multiple OPs occupy and describe the same space at the same time with varying properties.

In this article, the OPs are assumed to propagate with the speed of the free-stream wind rather than the effective wind speed in accordance with Taylor's frozen turbulence hypothesis (Taylor, 1938). The decision is supported by experimental results from Schlipf et al. (2010) and has also been used by other similar codes, e.g., Grunnet et al. (2010). This also solves the issue of the overlapping wake areas since neighboring OPs travel at the same speed and follow the same state changes. Another implication of this adaptation is that OPs no longer need to calculate the influence of foreign wakes at every time step. This would be used to determine their effective wind speed and thus how far they travel downstream in one time step. The only OPs which need to calculate the foreign influence are the ones at the rotor plane in order to determine the effective wind speed according to Eq. (6). These model assumptions also significantly decrease the computational load during the simulation. The downside of the change is that the effects of state changes now arrive too fast and abrupt at downstream turbines, which will be seen and discussed with the simulation results in Sect. 3. In future work, the wake propagation speed could be a tuning parameter which is set depending on atmospheric conditions such as the turbulence intensity for instance (Andersen et al., 2017).

2.4 Including directional dependency and observation point propagation

In this section, we address how the OPs, and therefore the wakes, react to a wind direction change. We assume that a wind direction change only affects the wake orientation and that the wake structure and downstream evolution (as defined by the underlying FLORIS model) can be seen independent from the free-stream behavior. It is therefore possible to split the two aspects into two coordinate systems: the world coordinate system 𝒦₀ and the wake coordinate system 𝒦₁. The free-flow conditions are described in 𝒦₀, whereas the wake properties are described in 𝒦₁. An OP links these two coordinate systems.

Figure 3 This figure visualizes the working of Eq. (9), which is applied for each OP individually. In panels (a) and (b), the position update of an OP in a time step with a constant wind direction is depicted. Panels (c) and (d) show the position update when the wind direction changes. In this case, the wake coordinate system is rotated around the OP's location to match the new downstream direction. This causes the apparent origin of the wake in the world coordinate system to change, which is visualized by the grey turbine.

The underlying FLORIS model is described in 𝒦₁, where the origin $x_{\mathrm{1}}=y_{\mathrm{1}}=z_{\mathrm{1}}=\mathrm{0}$ is located in the center of the rotor plane. The downwind distance is denoted as x₁, y₁ describes the horizontal crosswind distance and z₁ describes the vertical one. 𝒦₀ does not have a special orientation apart from z₀=0 being the ground level and the z₀ axis pointing upwards. In this work, x₀ describes the west–east axis, and y₀ describes the south–north axis. To transform a location vector r₁, described in 𝒦₁ of a turbine with the rotor-center location t₀, into r₀ the rotational matrix R₀₁ is used:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\mathit{r}}_{\mathrm{0}}& =\left[\begin{array}{c}{x}_{\mathrm{0}}\\ y_{\mathrm{0}}\\ z_{\mathrm{0}}\end{array}\right]={\mathit{t}}_{\mathrm{0}}+{\mathbf{R}}_{\mathrm{01}}\left(\mathit{\phi }\right){\mathit{r}}_{\mathrm{1}}=\left[\begin{array}{c}{x}_{\mathrm{0},\mathrm{T}}\\ y_{\mathrm{0},\mathrm{T}}\\ z_{\mathrm{0},\mathrm{T}}\end{array}\right]\\ & +\left[\begin{array}{llc}\cos \mathit{\phi }& -\sin \mathit{\phi }& \mathrm{0}\\ \sin \mathit{\phi }& \cos \mathit{\phi }& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}\end{array}\right]\left[\begin{array}{c}{x}_{\mathrm{1}}\\ y_{\mathrm{1}}\\ z_{\mathrm{1}}\end{array}\right].\end{array}\end{array}$$

This equation assumes a uniform wind direction φ at every location. This will not be the case for the formulation used for the OP propagation later on in Eq. (9). Each OP has two location vectors, r_0,OP and r_1,OP, one for each coordinate system. The OP's position update and its reduction factor is calculated in 𝒦₁. 𝒦₀ is used to calculate the wake interaction and to determine the wind speed, the wind direction and the ambient turbulence intensity. At the OP's creation, r_1,OP is determined by the Eqs. (4) and (5) for the crosswind coordinates; the downwind coordinate is set to 0. Its world location, r_0,OP, is then determined by Eq. (7) with the wind direction φ_0,T at the turbine location. To iterate the location of an OP from time step k to time step k+1, the downwind step is calculated first in 𝒦₁:

$$\begin{array}{}\text{(8)}& x_{\mathrm{1},\mathrm{OP}}(k+\mathrm{1})=x_{\mathrm{1},\mathrm{OP}}\left(k\right)+u_{\mathrm{OP}}\mathrm{\Delta }t,\end{array}$$

where Δt is the time step duration and u_OP is the magnitude of the wind vector u_0,OP at the OP's location r_0,OP. The direction will be applied in Eq. (9). For the scope of this work, u_0,OP can only have non-zero components in the x₀ and y₀ direction. With $x_{\mathrm{1},\mathrm{OP}}(k+\mathrm{1})$ the new crosswind locations $y_{\mathrm{1},\mathrm{OP}}(k+\mathrm{1})$ and $z_{\mathrm{1},\mathrm{OP}}(k+\mathrm{1})$ can be calculated with the Eqs. (2) and (3), respectively. This completes the transition ${\mathit{r}}_{\mathrm{1},\mathrm{OP}}\left(k\right)\to {\mathit{r}}_{\mathrm{1},\mathrm{OP}}(k+\mathrm{1})$. Note that only x_1,OP(k) is needed to determine the OP's location in 𝒦₁. At the cost of calculating y_1,OP(k) and z_1,OP(k) again at each time step, they do not have to be stored as states. To update r_0,OP(k), the step which the OP took in 𝒦₁ has to be translated into 𝒦₀:

$$\begin{array}{}\text{(9)}& {\mathit{r}}_{\mathrm{0},\mathrm{OP}}(k+\mathrm{1})={\mathit{r}}_{\mathrm{0},\mathrm{OP}}\left(k\right)+{\mathbf{R}}_{\mathrm{01}}\left({\mathit{\phi }}_{\mathrm{0},\mathrm{OP}}\right)\left[{\mathit{r}}_{\mathrm{1},\mathrm{OP}}(k+\mathrm{1})-{\mathit{r}}_{\mathrm{1},\mathrm{OP}}\left(k\right)\right],\end{array}$$

where φ_0,OP is the wind direction at r_0,OP(k). Note that φ_0,OP refers to one OP's individual wind direction; other OPs may have different values. This means that each OP propagates on its own and non-uniform wind directions can be simulated. Figure 3 shows the OP step in the wake and world coordinate system. In Fig. 3a and b the wind direction is constant, indicated by the arrow left to the y₀ axis. The OP calculates its step in the wake coordinate system (dotted arrow) and updates its location vectors. These are here simplified to r₀ and r₁. In Fig. 3c the wind direction changes, and the former FLORIS wake description is invalid and greyed out. With the new wind direction R₀₁(φ_0,OP) is calculated differently. The OP can calculate its step in the wake coordinate system as before, but its translation 𝒦₁→𝒦₀ changed. Note that neither r₀ nor r₁ is influenced by the changed wind direction. Their magnitude and orientation remain the same in their respective coordinate systems; however, their orientation towards each other changes.

2.5 Calculation of C_T and C_P

The thrust coefficient C_T is often approximated following the actuator disc theory: $C_{\mathrm{T}}\left(a\right)=\mathrm{4}a(\mathrm{1}-a)$, where a is the axial induction factor. To circumvent this approximation, simulations or experiments can be used to create lookup tables. Since most equations of the Gaussian FLORIS model are dependent directly on C_T rather than a, we used lookup tables generated in SOWFA to align FLORIDyn's thrust coefficient with what the turbines in the validation environment experience. For completeness, we also use lookup tables for the power coefficient C_P. The tables are generated for the DTU 10 MW reference turbine (Bak et al., 2013). It has to be added that these tables are generated from a grid of high-fidelity simulations, where the coefficients were read after the simulation converged to a steady state. The tables can therefore only approximate the effect a changing turbine state and changing wind field conditions onto C_T and C_P. Control approaches for axial-induction-based controllers, such as the one presented by Annoni et al. (2016), successfully use similar lookup tables, which is why we assume these to be sufficient. Nevertheless, an extension for dynamic circumstances would be a valuable addition for future work but is also connected to a significant computational effort.

In the tables, the coefficients are described dependent on the blade pitch angle β and the tip speed ratio λ(ω,u_eff), where ω is the angular velocity of the rotor. However, neither FLORIS nor FLORIDyn can provide λ and β. What they can provide is u_eff. Combined with the assumption that each turbine follows a greedy control strategy and maximizes C_P(λ,β) for the given wind, we can formulate the coefficients dependent only on u_eff: first maximize C_P within the physical limitations of the wind turbine for all wind speeds, and then use the λ_P,max and β_P,max to calculate the respective C_T. The resulting curves can be seen in Fig. 4.

Figure 4 Greedy control settings of the un-yawed 10 MW DTU reference turbine based on the effective rotor wind speed.

Unfortunately, the resulting C_T(u_eff) values can get very high, especially for low wind speeds. This conflicts with some FLORIS equations which comprise the term $\sqrt{\mathrm{1}-C_{\mathrm{T}}}$ and become complex for C_T values above 1. To avoid these issues, C_T(u_eff) is limited to its value at the Betz limit: $C_{\mathrm{T}}{|}_{a=\mathrm{1}/\mathrm{3}}=\mathrm{0}.\stackrel{\mathrm{‾}}{\mathrm{8}}$ (Bianchi et al., 2007). Another complication is the calculation of the added turbulence levels as it is the only equation which requires the axial induction factor. In this case, the calculation of C_T(a) was inverted to determine a(C_T), based on the actuator disc theory, as follows:

$$\begin{array}{}\text{(10)}& a={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{1}-\sqrt{\mathrm{1}-C_{\mathrm{T}}}\right).\end{array}$$

Yaw effects on C_T and a are neglected here. In future work this expression could be substituted, for instance by the polynomial approximation of Madsen et al. (2020). It extends a(C_T) to C_T values above 1. However, as C_T is limited in this work, this extension is not necessary. The power coefficient is the remaining aspect which was used unaltered from the lookup tables. For the tested wind speeds below 11 m s⁻¹, the power coefficient is constant at C_P=0.4929. The effect of γ is approximated by multiplying C_P with $\cos (\mathit{\gamma }{)}^{p_{\mathrm{p}}}$. For simplicity's sake we assume p_p to be a constant value. This could be extended by the work presented by Liew et al. (2020) which takes the presence of other wakes into account. Similarly, Howland et al. (2020) presents an adaptation for locally varying wind profiles.

2.6 Wind field model

In order to drive the FLORIDyn model, the wind field needs to be able to simulate heterogeneous, changing environmental conditions. The implemented solution is inspired by the work of Farrell et al. (2021). The basic assumption is that measurements of the wind field variables are available at certain locations. This could be due to satellite data, lidar measurements, met masts or other sensors. The value of a measurement for the location of an OP is then interpolated between the measurements available. To reduce the computational effort of an interpolation at every time step, a nearest-neighbor interpolation (NNI) is desirable. To get a sufficient resolution of the measurements to justify a NNI, the sparse measurements m have to be mapped to dense measurement grid points m_g:

$$\begin{array}{}\text{(11)}& {\mathit{m}}_{\mathrm{g}}=\mathbf{M}\mathit{m},\end{array}$$

where the matrix M describes the mapping and can be calculated offline. The ith row in M describes the percental composition of m_g,i from m. As a result, the sum of every row in M is equal to 1. This way, a more complex interpolation can be reduced to a matrix multiplication and a NNI at runtime. In this work, a linear interpolation is used to map the measurements to the grid points, which are spaced in a 20×20 m grid. OPs outside of the grid defined by m_g use the closest grid point. This method is also independent from the quantity measured. In this work, the wind speed, the wind direction and the ambient turbulence intensity were interpolated with the presented method.

Table 1 Parameters used in the simulation with the values they influence.

However, the presented method is only meant for values changing in the x₀ and y₀ direction. The wind speed is the only field measurement which is also changed in the z₀ direction; wind direction and ambient turbulence intensity are assumed to be constant in the vertical direction. Following Farrell et al. (2021) the power law is applied:

$$\begin{array}{}\text{(12)}& u\left(z_{\mathrm{0}}\right)={\left({\displaystyle \frac{z_{\mathrm{0}}}{z_{\mathrm{0},\mathrm{m}}}}\right)}^{{\mathit{\alpha }}_{\mathrm{s}}}u\left(z_{\mathrm{0},\mathrm{m}}\right),\end{array}$$

where z_0,m is the height of the measurement and α_s is the shear coefficient. The shear coefficient approximates the combined effect of atmospheric stability and surface roughness. A small value describes unstable flow conditions. Examples for characteristic α_s values due to surface roughness are 0.11 over water, 0.16 over grass, 0.20 over shrubs, 0.28 over forests and 0.40 over cities (Emeis, 2018). In this work z_0,m is equal to the hub height z_h of the turbine.

3 Simulation results

In this section, the Gaussian FLORIDyn model is compared to SOWFA with the focus on turbine interaction. Two wind farm layouts are considered for comparison: three consecutive turbines and a nine-turbine cluster arranged in a 3×3 configuration. The DTU 10 MW reference turbine is used for all simulations. Table 1 summarizes the FLORIS and FLORIDyn parameters used in the simulations. The FLORIS parameters k_a and k_b are from Niayifar and Porté-Agel (2015), k_f,a to k_f,d are set based on FLORISSE_M (Doekemeijer et al., 2021), and α^* and β^* follow the findings of Bastankhah and Porté-Agel (2016). The efficiency η was tuned based on turbine T0 in the three-turbine baseline case; p_p was tuned based on the three-turbine yaw case (Sect. 3.1.1 and 3.1.2 respectively). For FLORIDyn, n_c relates to the number of OP chains per turbine and n_OP to the number of OPs per chain. The value of n_OP was set to cover the entire relevant downstream domain of a turbine; n_c was set to maintain a sufficient density of OPs at the location of other turbines. In FLORIDyn, one time step is 4.0 s long.

Table 1 also includes the wind shear coefficient, α_s, which was approximated based on the free flow in SOWFA. The inflow boundary conditions for SOWFA are provided by a precursor simulation which simulates a horizontally homogenous, conventionally neutral atmospheric boundary layer including Coriolis effects. The SOWFA settings differ for the three-turbine case and the nine-turbine case and will be explained in the respective sections.

3.1 Three-turbine case

The three turbines are placed 5D=892 m apart from each other in downwind direction. Turbine T0 is located at (608, 500 m), and T1 and T2 are at (1500, 500 m) and (2392, 500 m) respectively. The mean wind speed at hub height is approximately 8.2 m s⁻¹ with an ambient turbulence intensity of roughly 6 %. The mean wind direction is constant along the x axis. The full SOWFA flow field domain spans $\mathrm{3}\times \mathrm{1}\times \mathrm{1}$ km and was simulated with a time step of Δt=0.04 s. The base cells of the flow field are $\mathrm{10}\times \mathrm{10}\times \mathrm{10}$ m. The refinement areas are centered in the domain and have no offset from the ground. The first refinement is $\mathrm{2.4}\times \mathrm{0.8}\times \mathrm{0.5}$ km with $\mathrm{5}\times \mathrm{5}\times \mathrm{5}$ m cells, the second one is $\mathrm{2.2}\times \mathrm{0.6}\times \mathrm{0.35}$ km with $\mathrm{2.5}\times \mathrm{2.5}\times \mathrm{2.5}$ m cells. Figure 5 shows the to-scale layout including the areas of cell refinement. In SOWFA, the turbines are modeled with the built-in actuator line method (ALM) implementation (Sorensen and Shen, 2002).

Figure 5 Scaled layout of the three-turbine case with the wind direction indicated by an arrow on the left. The 0, 10 and 20^∘ yaw orientations from T0 are indicated as turbine symbols with the according orientation. The colored background areas indicate the zones of cell refinement.

To give a better idea of the low-frequency, less-turbulent dynamics, the power generated in SOWFA is also presented filtered by a zero-phase second-order low-pass filter. This non-causal filter is added to aid the visual interpretation of the simulation results. The filter has a damping ratio of d=0.7 and a natural frequency of ω=0.03 s⁻¹. This allows for a more equal comparison as FLORIDyn is sampled at a lower frequency and turbulence is only included as a flow field metric.

Figure 6 Wind farm start-up and steady-state comparison of the power generated in SOWFA and in FLORIDyn. The unfiltered data are plotted in panel (a), filtered data are plotted in panel (b). The SOWFA data are filtered with a zero-phase (noncausal) low-pass filter and FLORIDyn with a causal low-pass filter.

A regular second-order low-pass filter with the same d and ω is used for the FLORIDyn data. This causal filter visualizes how low-pass filtering would affect the predicted power generated. This could have advantages due to the changes made to the OP travel speed in Sect. 2.3.3, which can lead to a very abrupt wake interaction, as will be discussed in Sect. 3.1.1. However, the filter also naturally adds a phase shift to the data, an effect which might not be desired.

Note that the two filters have different purposes: the non-causal SOWFA filter aims to help to interpret the simulation results, while the causal FLORIDyn filter explores if and when the use of a low-pass filter would be advantageous or if it would decrease the quality of the results.

3.1.1 Comparison of the wind farm start-up and steady state

In Fig. 6 the power generated by the turbines in FLORIDyn is compared to the SOWFA simulation. The dynamics in this simulation are the turbulent wind field and the settling of the wake. In the unfiltered data, the interaction in FLORIDyn sets in earlier and more abrupt than in SOWFA. This is due to the OPs traveling at the free wind speed, as explained in Sect. 2.3.3. The slight curvature of the drop at t≈100 s can be explained by the wind shear: OPs at a lower altitude travel at a slower free wind speed than OPs at a higher altitude and therefore arrive later at the downstream turbine and therefore affect the turbine later. There are two major aspects to address in order to close the gap between the SOWFA and the FLORIDyn start-up: on the one hand, the way state changes propagate through the wake; on the other hand, how a downstream turbine reacts to the new wind field. To give an idea of how a change of the latter aspect would influence the plot, Fig. 6 also shows low-pass-filtered FLORIDyn data in comparison to zero-phase-filtered SOWFA data. The FLORIDyn data align much more with the SOWFA data but still show discrepancies in terms of dynamic response and steady-state quality of the solution. It should be emphasized that the in-FLORIDyn-applied filter does not affect the wake, it only adds an artificial dynamic response to the power calculation. This is important when heterogeneous and changing wind directions are taken into account.

The power generated after the wind farm start-up remains steady in FLORIDyn. This is because there are no turbulent wind speed changes in FLORIDyn. In this state, FLORIDyn is equal to the underlying FLORIS model; therefore, errors in this state need to be solved by adapting the FLORIS model. This could be done by parameter tuning, which has only been done partially in this work (η and p_p; see introduction Sect. 3).

To incorporate the turbulent wind speed changes, at least to a certain degree, an estimation of the wind speed at the turbine location would be necessary. This could be done by including wind speed sensor data or estimating the wind based on the power generated (Gebraad et al., 2015) or based on a torque balance equation (Ortega et al., 2013).

A notable aspect of this simulation is the influence of the added turbulence. Because T0 adds turbulence to the wind field, the wake of turbine T1 recovers faster. Turbine T2 thus experiences higher wind speeds and generates more power than it would without the additional turbulence. This effect can also be observed in the SOWFA data. The old Zone FLORIDyn model is not able to capture this effect due to the underlying Zone FLORIS model. It shows how the FLORIDyn framework is inherently dependent on the capabilities of the employed FLORIS model.

3.1.2 Comparison during a yaw angle change

In this simulation, the yaw angle γ of turbine T0 is changing from 0 to 20^∘ in steps of 10^∘, starting at t=200 and 800 s. The change rate of γ is set to 0.3^∘ s⁻¹. Figure 7 shows the unfiltered SOWFA data in comparison to the unfiltered FLORIDyn data on the left, as well as the filtered data on the right. Filtering was performed as described in the introduction of Sect. 3.1.

Figure 7 Comparison of the power generated in SOWFA and in FLORIDyn with changing yaw angles. The transparent bars indicate the time window in which turbine T0 increases its yaw angle by 10^∘. Panel (a) shows the unfiltered data, and panel (b) shows the filtered data. The SOWFA data are filtered with a zero-phase (noncausal) low-pass filter and FLORIDyn with a causal low-pass filter.

Figure 8 Comparison between the zero-phase-filtered (zp.f.) SOWFA data in the baseline case (bl.) and in the yaw case. (a) In absolute values; (b) with the difference between the yaw case and the baseline case. The transparent bars indicate the time window in which turbine T0 increases its yaw angle, first from 0 to 10^∘ and then to 20^∘.

In FLORIDyn, turbine T1 shows a slight reaction to the yaw changes of turbine T0 at roughly t≈320 s and more significantly at t≈920 s. The influence of the state change then travels further and impacts T2 at t≈430 s and, as well more significantly, at t≈1030 s. In SOWFA, the reaction is obscured by turbulent influences. However, an increase in average power can be seen for T1 and T2 throughout the entire simulation. Figure 8 shows the baseline simulation in comparison to the SOWFA simulation, in absolute values and the difference. The data of both simulations can be compared since they use the same wind field. It allows for a more accurate determination of the reaction time to the upstream change. Turbine T1 starts to react to the first yaw angle change at $t\approx \mathrm{320}\pm \mathrm{8}$ s, T2 at $t\approx \mathrm{434}\pm \mathrm{8}$ s. Given the distance to T0, this translates to a travel speed of the first influence between [6.98, 7.97] m s⁻¹ to T1 and [7.38, 7.90] m s⁻¹ to T2. This indicates that first effects of the yaw angle change do travel at almost free-stream velocity and the times align with the FLORIDyn prediction. However, FLORIDyn does lack the dynamic nature of the interaction, which means the response of the wake to the state change of the upstream turbine and the response of the downstream turbine to changes in the wake. Given that all OPs travel at their free-stream velocity, turbine state changes are directly picked up by the OPs and transported, and the FLORIDyn turbine reacts immediately when the OP arrives. The low-pass-filtered results provide an idea of how a dynamic response could change the results. The unfiltered difference between the SOWFA simulations is given in the Appendix A1. Figure 7 also shows that FLORIDyn underestimates the gain in generated energy in the steady-state region. The error likely lies with the underlying FLORIS model as it is a steady-state error. Better parameter tuning would likely decrease or even eliminate the error.

Figure 9 Complete nine-turbine FLORIDyn flow field in comparison to SOWFA at t=600 s. The wind direction change is indicated in the lower-left corner of the SOWFA plot.

3.2 Nine-turbine case

In order to test the model in a changing environment, a simulation with nine turbines was performed. The turbines are arranged in a three-by-three grid with 900 m distance to the closest turbines and 600 m to the edge. The setup is presented in Fig. 9 as well as the numbering of the turbines. The wind field performs a 60^∘ uniform wind direction change from 15 to 75^∘, as indicated in Fig. 9. The change starts at t=600 s with 0.2^∘ s⁻¹ and ends at t=900 s. The change in wind direction is achieved by using SOWFA's built-in utility to specify the wind speed and wind direction at a certain height and time. For the remainder of the simulation, the wind field conditions remain steady. To keep the computational load of the SOWFA simulation low, the DTU 10 MW turbines were simulated with the actuator disc method (ADM). This also allows a coarser grid and time resolution: the domain is discretized in $\mathrm{10}\times \mathrm{10}\times \mathrm{10}$ m cells, and the SOWFA time step length is set to 0.5 s. The flow field spans $\mathrm{3}\times \mathrm{3}\times \mathrm{1}$ km. The average wind speed during the simulation is 8.2 m s⁻¹, and the ambient turbulence intensity is approximately 6 %. During the simulation, the turbines maintain a yaw angle of 0^∘ and turn with the wind. For simplicity we assume ideal wind direction tracking capabilities and apply a prescribed motion. For more information see the dataset which contains the SOWFA files for the case and the precursor simulation (Becker, 2022b).

Figure 10 shows the flow field during the wind direction change, starting at the time instances t=700, t=800 and 900 s. The SOWFA slices are taken at hub height. To show the center of the FLORIDyn flow field, only OPs between $\mathrm{0.5}{z}_{\mathrm{h}}<z_{\mathrm{OP}}<\mathrm{1.5}{z}_{\mathrm{h}}$ are plotted. As the wake expands, more chains leave these bounds, which leads to a sparser description of the wake in the plot. Due to the changes to FLORIDyn described in Sect. 2.3.3, the OPs do not influence each other in the field and OPs with a higher velocity can appear among OPs with a lower velocity. The net effect of multiple OPs is only calculated at the rotor plane. The grey arrows in the FLORIDyn plots indicate the current wind direction. The plots of both simulations visualize how the wakes slowly transition to the new wind direction, forming a bow shape in the process. FLORIDyn seems to describe the general path of the SOWFA wakes quite well. It also capture some effects like shorter, wider wakes of T4 and T5 at t=900 s.

Figure 10 Nine-turbine flow field at hub height during the wind direction change at t=700 s (a, b), t=800 s (c, d) and t=900 s (e, f). The FLORIDyn flow field is on the left and includes grey arrows as an indicator of the current wind direction. On the right is the corresponding snapshot from the SOWFA simulation.

Figure 11 Power generated in the nine-turbine case. The grey area marks the time window in which the wind direction linearly changes by 60^∘. The plots are arranged to fit the layout of the wind farm in Fig. 9. The data show the zero-phase-filtered (zp. f.) and the unfiltered (unf.) SOWFA data, as well as the filtered (f.) and unfiltered FLORIDyn data.

Table 2 Points in time at which the power generated in the nine-turbine case is minimal due to wake interaction.

To more accurately judge the timing of FLORIDyn, Fig. 11 shows the generated power of all nine turbines. The plots are arranged in the same way as the turbines in the flow field plots. All plots show the filtered and unfiltered data of FLORIDyn and SOWFA. The filtering is identical to the filtering in Sect. 3.1. The grey area marks the time window of the wind direction change. Looking at the magnitude of the generated power, FLORIDyn predicts the average of the free-stream turbines T0, T3, T6, T7 and T8 quite well, but the remaining turbines show a noticeable offset in generated power. This could be due to speed-up effects and is briefly discussed in Appendix A2. The interesting aspect is the timing of the wake interaction from upstream turbines with downstream turbines. The generated power by T4 shows the passing of the wake of T6 during the wind direction change. Noticeable is the accuracy with which the unfiltered FLORIDyn data align the unfiltered SOWFA data. Table 2 lists the points in time at which the power generated is minimal in SOWFA and in FLORIDyn, as well as the difference. This shows that FLORIDyn predicts the maximal wake influence 5.5 s later than in SOWFA.

Table 3 Computational performance.

The filtering of the FLORIDyn data significantly worsens the quality of the result, in contrast to the simulations of the three-turbine case. The filtering was applied on the data from the rotor plane. Thus, only modifying the way a turbine perceives the incoming, foreign wake in FLORIDyn will not improve the simulation for changing environments. As a result, future research has to improve the way a turbine dynamically influences its own wake. Turbine T1 shows similar behavior to T4: the generated power shows the wake influence of T3 first, but it also shows, after the wind direction stopped changing, the influence of the outskirts of the wake of T6. While the timing of this interaction shows good agreement, the magnitude of the interaction is considerably lower in FLORIDyn than in SOWFA. This could be due to a too fast recovering FLORIS wake, an inadequate wake superposition method or due to local turbulence levels, which FLORIDyn can not capture. T5 shows the overlapping influence of the wakes of T6 and T7. The two overlapping Gaussian influences do form a longer period of reduced generated power. This can be seen in both simulations. Table 2 shows the largest timing error between SOWFA and FLORIDyn for this wake influence. This could stem from an inaccurate wake interaction model and the way added turbulence is treated. T2 shows the most overlapping influences by the wakes of T3, T4, T6 and T7 in this order. While the first two overlapping interactions show good agreement, FLORIDyn shows poorer agreement with SOWFA for the influence of T6 and T7. Again, the SOWFA simulation suggests a larger decrease in generated power. A reason could be that the way wakes combine is not described accurately enough. The timing of the wake of T7 seems to be a bit too late as well. However, the SOWFA simulation recovers to its steady state at about the same time as FLORIDyn. Additionally, all turbines with the exception of T6 experience a small influence of the upstream turbines in the steady-state configuration. This effect is not noticeable in the SOWFA simulations. Concluding, FLORIDyn describes the timing of passing wake influences quite well. However, there are discrepancies in terms of magnitude and possibly in the way wakes combine their effects on downstream turbines.

3.3 Computational performance

Table 3 contains the average computational time per time step, which is equivalent to 4 s simulation time. This can be compared to SOWFA, which can take around 5.8×10² to 5.4×10³ s per core per time step, depending on the setup (van den Broek and van Wingerden, 2020). The FLORIDyn measurements were performed for two and three consecutive turbines and a 2×2 and 3×3 turbines wind farm. The times exclude plotting and the simulation setup time. A setup can take up to 3 s, depending on how much data need to be imported. The measurements were taken on a MacBook Pro (2019) with a 2.3 GHz eight-core Intel i9 CPU, 32 GB of 2667 MHz DDR4 RAM, an SSD, and MacOS Catalina (10.15.7). The simulation environment is MATLAB 2020a without the use of toolboxes, such as the parallel computing toolbox, and without precompiled code, besides what is built into the simulation environment. These results naturally vary with the layout, atmospheric behavior, simulation settings, etc. and are only meant to give an estimation of the performance.

A first takeaway is that FLORIDyn simulates all cases faster than real time: within 4 s, the simulation can perform between 164 and 6.5 simulation steps, depending on the number of turbines simulated. This results in 656 s_Sim in-simulation time for two turbines and goes down to 26 s_Sim for nine turbines in 4 s of real time. This opens up the needed computational headroom for a model-based real-time control strategy and the necessary optimization. On the other hand, the times also do not offer a large time window for optimization. For instance, in the three-turbine yaw case, it takes roughly 300 s_Sim in simulated time until the yaw changes have propagated from the first turbine to the last turbine. To optimize the control actions with this model for the near future, parallel computing would be needed. With an increasing number of turbines, this time window decreases.

Generally, the computational time increases quadratically with $n_{\mathrm{T}}^{\mathrm{2}}-n_{\mathrm{T}}$, as n_T turbines need to determine if they are in the wake of another turbine and calculate the influence. This growth in computational effort is assumed to decrease with larger wind farms: as the spacial dimension grows, not every turbine needs to consider all other turbines for interaction. Nevertheless, the simulation times will exceed what is practical for a wind farm with ≫3 turbines.

There are multiple opportunities to improve performance which have not been utilized so far. The main aspect which increases computational effort is the interaction among the turbines. A first step can be to calculate the turbine interactions in parallel. A second step is to find a way to efficiently determine if a turbine is influenced by a wake. Furthermore, the number of OPs per turbine can be decreased and tuned: not all chains need to be equally long, and OPs which wander out of the domain can be disregarded. Then, there is the fundamental question of whether the proposed structure of FLORIDyn can be improved, for instance by using less but more efficient OPs. Eventually, the programming platform can be switched to a choice which allows more specific optimization, e.g., C, C$++$ or Julia.

4 Conclusions and recommendations

In this paper, a new FLORIDyn model is presented and compared to SOWFA simulations. This model utilizes a Gaussian FLORIS model and concepts from the previously published FLORIDyn framework by Gebraad and van Wingerden (2014) to create a three-dimensional, dynamic and computationally lightweight wind farm model. The new FLORIDyn model is further capable of simulating its wakes under heterogeneous and time-varying flow conditions. To achieve this, we presented a mathematical approach to decouple wake and flow characteristics into two coordinate system which are connected by observation points. To simulate changing environmental conditions, a method to map sparse flow field measurements to a finer grid was presented which avoids the interpolation cost at runtime. The new FLORIDyn model shows good performance compared to SOWFA in terms of timing and is able to predict accurately when a downstream turbine is experiencing influences from upstream turbines.

Despite the considerable advancements over the old FLORIDyn implementation, there are still several aspect of the model which can be improved. The central aspect is how turbines influence wakes and how wakes are perceived by turbines. In this work we have decoupled the OP propagation speed from the effective wind speed, which effectively leads to a simpler, lightweight model while the wake behavior is still dynamically described. However, this way state changes reach downstream turbines too soon and in a sudden manner. Ideally this can be overcome by finding better, computationally lightweight methods to model the influence of changing turbine states on the wake needs and also how a turbine reacts to dynamic changes in the flow. Another aspect that can be improved is related to the interface between FLORIDyn and FLORIS. FLORIS has been subject to many developments and improvements, and FLORIDyn can utilize these improvements if it improves the interface: with a generic interface, newer developments can be included and existing code can be used in a sustainable manner. The simulations also show that parameter tuning has to be more accessible and possibly needs to be performed online in some cases. The next aspect which could be improved is the coupling of FLORIDyn with the turbulent environment of the real wind farm (or its surrogate). Combined with the changes to the OP propagation speed from this work, this can lead to a more uneven OP distribution with dense areas where high wind speeds decrease and sparse areas where low wind speeds increase. An extension to the model could feature a method to combine and generate OPs, depending on the density of OPs. Although, this could also lead to undesired information loss, depending on the implementation. To achieve better results, the wind field model has to be replaced or enhanced by estimators. The latter would provide a more accurate estimate of wind speed, direction and ambient turbulence intensity for the FLORIDyn simulation. In the long term, a dynamic description of the environment could become part of the FLORIDyn model. This could also include effects like induction zones and speed-up between the turbines. The last aspect to consider for improvement is the performance. In its current implementation, FLORIDyn delivers its results at a low computational cost. This has to be maintained, if not improved, to allow its use for dynamic real-time closed-loop control algorithms in the future. The simulation also needs to be structurally improved to keep its low computational cost for wind farms with large numbers of turbines.

Concluding, the new FLORIDyn model is a promising concept with unique strengths. With FLORIS in its core, it utilizes an existing, successfully employed model and provides a new dimension in a challenging environment at a low computational cost. The model can already be adapted to work in a closed-loop control design and shows more potential if the mentioned aspects are improved.

Appendix A: Additional plots and aspects of the simulation results

A1 Unfiltered difference between yawed and baseline case

Figure A1 shows the difference between the power generated in SOWFA in the yaw case (Sect. 3.1.2) and the steady-state baseline case (Sect. 3.1.1). Both simulations are performed in the same turbulent environment, something which would be impossible to achieve in realistic conditions. This way, the difference allows for a clearer interpretation of the influence of the yaw step, at least the timing. In comparison to Fig. 8, Fig. A1 shows the unfiltered data as well as the filtered data for all three turbines. T1 shows between t=312 and 329 s a first reaction due to the changed wake of T0. T2 shows a first reaction between t=426 and 442 s. The filtered data show a slightly earlier influence due to the nature of a zero-phase filter.

Figure A1 Difference between the yaw case in SOWFA and the baseline case with unfiltered and zero-phase-filtered data. Filtering was performed before calculating the difference. The dotted lines mark the start of the yaw angle changes of T0.

A2 Averaged velocity in the nine-turbine case

Figure A2 shows the averaged wind speed in the nine-turbine case from t=500 to 600 s at hub height in SOWFA. A total of 34 slices were used to average. The wind speed is binned into 11 wind speed sections which are plotted as contours. During the time of averaging the wind direction is constant. The average wind speed of the incoming air is approximately 8 m s⁻¹. However, between the turbine rows, the wind speed increases to a higher level, up to 9.44 m s⁻¹ in some places. This could be explained by speed-up effects: the turbines act as resistances in the flow field, and the wind speed in the place of least resistance, between the turbines, increases. The effect has been observed and described in Bastankhah et al. (2021) as well for instance. Due to the speed-up, the turbines further downstream experience higher wind speeds than the ones in free stream and generate more energy. Figure 11 quantifies the effect, where T2, T4 and T5 generate significantly more energy than T6 for instance. After the wind direction change, the effect leads again T2 and T4 to generate more power. T1 is now in the situation T5 was in initially and also generates more energy. T5, however, drops to a lower level. Without an added model for effects like these, FLORIDyn (and FLORIS) will not be able to accurately describe the wind field.

Figure A2 Averaged wind speed from t=500 to 600 s, divided into 11 speed sections.