The revised FLORIDyn model: implementation of heterogeneous ﬂow 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 ﬂow dynamics: FLORIS neglects these and provides a computationally very cheap approximation of the mean wind farm ﬂow. FLORIDyn deﬁnes 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.


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 lowcomputational-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 SimWind-Farm 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 for-mulation 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 highfrequency 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 onedimensional 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 sixturbine 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 FLORI-Dyn 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.

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.
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 FLORI-Dyn 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, K 1 , x 1 describes the downwind direction, y 1 the horizontal crosswind direction and z 1 the vertical crosswind direction (Fig. 1).In this coordinate frame, the rotor center is always located at [0, 0, 0] .This coordinate system is not to be confused with the longitudinal (x 0 ), latitudinal (y 0 ) and vertical (z 0 ) world coordinate system K 0 .

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(γ ) p 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 = u/u 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 https://doi.org/10.5194/wes-7-2163-2022 Wind Energ.Sci., 7, 2163-2179, 2022 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 0 , 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.

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

Changes to the FLORIDyn approach
Due to the changed underlying FLORIS model, the FLORI-Dyn 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 K 1 , 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.

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 (ν y , ν z ) ∈ [−0.5, 0.5] for a given number of OPs.We used 50 OPs per time step.To Wind Energ.Sci. , 7, 2163-2179, 2022 https://doi.org/10.5194/wes-7-2163-2022cover 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: y 1,OP ν y,OP , σ y , w y,pc , δ = ν y,OP 6σ y + w y,pc + δ, (2) z 1,OP ν z,OP , σ z , w z,pc = ν z,OP 6σ z + w z,pc . (3) 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 0 < x 1 ≤ x c and σ y,fw for x 1 > x c .Respectively, σ z is defined the same way.
The variable δ describes the deflection of the centerline.If OPs travel below z 1 = 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.

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: (5) 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
nOP ] the effective wind speed at the rotor plane is calculated as follows: 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 1 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.

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

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 K 0 and the wake coordinate system K 1 .The free-flow conditions are described in K 0 , whereas the wake properties are described in K 1 .An OP links these two coordinate systems.
The underlying FLORIS model is described in K 1 , where the origin x 1 = y 1 = z 1 = 0 is located in the center of the rotor plane.The downwind distance is denoted as x 1 , y 1 describes the horizontal crosswind distance and z 1 describes https://doi.org/10.5194/wes-7-2163-2022 Wind Energ.Sci., 7, 2163-2179, 2022  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 vertical one.K 0 does not have a special orientation apart from z 0 = 0 being the ground level and the z 0 axis pointing upwards.In this work, x 0 describes the west-east axis, and y 0 describes the south-north axis.To transform a location vector r 1 , described in K 1 of a turbine with the rotor-center location t 0 , into r 0 the rotational matrix R 01 is used: 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 K 1 .K 0 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 K 1 : 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 0 and y 0 direction.With x 1,OP (k + 1) the new crosswind locations y 1,OP (k + 1) and z 1,OP (k + 1) can be calculated with the Eqs.( 2) and (3), respectively.This completes the transition r 1,OP (k) → r 1,OP (k + 1).Note that only x 1,OP (k) is needed to determine the OP's location in K 1 .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 K 1 has to be translated into K 0 : 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 0 axis.The OP calculates its step in the wake coordinate system (dotted arrow) and updates its location vectors.These are here simplified to r 0 and r 1 .In Fig. 3c the wind direction changes, and the former FLORIS wake description is invalid and greyed out.
With the new wind direction R 01 (ϕ 0,OP ) is calculated differently.The OP can calculate its step in the wake coordinate system as before, but its translation K 1 → K 0 changed.Note that neither r 0 nor r 1 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.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 highfidelity 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.
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 √ 1 − C 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 T | a=1/3 = 0.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: 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 −1 , the power coefficient is constant at C P = 0.4929.The effect of γ is approximated by multiplying C P with cos(γ ) p 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.

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 nearestneighbor 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 : 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.However, the presented method is only meant for values changing in the x 0 and y 0 direction.The wind speed is the only field measurement which is also changed in the z 0 direction; wind direction and ambient turbulence intensity are assumed to be constant in the vertical direction.Following Farrell et al. (2021) 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.

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

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 −1 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 3 × 1 × 1 km and was simulated with a time step of t = 0.04 s.The base cells of the flow field are 10 × 10 × 10 m.The refinement areas are centered in the domain and have no offset from the ground.The first refinement is 2.4 × 0.8 × 0.5 km with 5 × 5 × 5 m cells, the second one is 2.2 × 0.6 × 0.35 km with 2.5 × 2.5 × 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).
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 −1 .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.
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 gen-Wind Energ.Sci., 7, 2163-2179, 2022 https://doi.org/10.5194/wes-7-2163-2022erated.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 noncausal 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.

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 FLORI-Dyn framework is inherently dependent on the capabilities of the employed FLORIS model.

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 −1 .Figure 7 shows the unfiltered SOWFA data in comparison to the unfiltered FLORI-Dyn 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.https://doi.org/10.5194/wes-7-2163-2022 Wind Energ.Sci., 7, 2163-2179, 2022  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 ≈ 320±8 s, T2 at t ≈ 434±8 s.Given the distance to T0, this translates to a travel speed of the first influence between [6.98, 7.97] m s −1 to T1 and [7.38, 7.90] m s −1 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.

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 −1 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 10 × 10 × 10 m cells, and the SOWFA time step length is set to 0.5 s.The flow field spans 3 × 3 × 1 km.The average wind speed during the simulation is 8.2 m s −1 , 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 0.5z h < z OP < 1.5z 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 FLORI-Dyn 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.
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   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.
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 gener- Total number of OPs 2 × 10 4 3 × 10 4 4 × 10 4 9 × 10 4 t comp.time per step (s) 2.44 × 10 −2 5.87 × 10 −2 1.09 × 10 −1 6.13 × 10 −1 ated power.This can be seen in both simulations.Table 2 shows the largest timing error between SOWFA and FLORI-Dyn 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.

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 2 to 5.4 × 10 3 s per core per time step, depending on the setup (van den Broek and van Wingerden, 2020).The FLORI-Dyn 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 Ma-cOS Catalina (10.15.7).The simulation environment is MAT-LAB 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 2 T − n 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.

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 FLORI-Dyn 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 https://doi.org/10.5194/wes-7-2163-2022 Wind Energ.Sci., 7, 2163-2179, 2022 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 FLORI-Dyn 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 pro-vides 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.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 −1 .However, between the turbine rows, the wind speed increases to a higher level, up to 9.44 m s −1 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, FLORI-Dyn (and FLORIS) will not be able to accurately describe the wind field.

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

Figure 2 .
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).

Figure 3 .
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.

Figure 4 .
Figure 4. Greedy control settings of the un-yawed 10 MW DTU reference turbine based on the effective rotor wind speed.
the FLORIS and FLORIDyn parameters used in the simulations.The FLORIS parameters k a and k b are fromNiayifar 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 ofBastankhah 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.1and 3.1.2respectively).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.

Figure 5 .
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.

Figure 6 .
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.

Figure 7 .
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 .
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 • .

Figure 9 .
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.

Figure 10 .
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 .
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.

Figure A1 .
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.

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

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

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