A control-oriented dynamic wind farm model: WFSim

. Wind turbines are often sited together in wind farms as it is economically advantageous. Controlling the ﬂow within wind farms to reduce the fatigue loads, maximize energy production and provide ancillary services is a challenging control problem due to the underlying time-varying non-linear wake dynamics. In this paper, we present a control-oriented dynamical wind farm model called the WindFarmSimulator (WFSim) that can be used in closed-loop wind farm control algorithms. The three-dimensional Navier–Stokes equations were the starting point for deriving the control-oriented dynamic wind farm model. Then, in order to reduce computational complexity, terms involving the vertical dimension were either neglected or estimated in order to partially compensate for neglecting the vertical dimension. Sparsity of and structure in the system matrices make this model relatively computationally inexpensive. We showed that by taking the vertical dimension partially into account, the estimation of ﬂow data generated with a high-ﬁdelity wind farm model is improved relative to when the vertical dimension is completely neglected in WFSim. Moreover, we showed that, for the study cases considered in this work, WFSim is potentially fast enough to be used in an online closed-loop control framework including model parameter updates. Finally we showed that the proposed wind farm model is able to estimate ﬂow and power signals generated by two different 3-D high-ﬁdelity wind farm models.


Introduction
Optimizing the control of wind turbines in a farm is challenging due to the aerodynamic interactions among turbines.These interactions come from the fact that downwind turbines are often operating in the wakes of upwind ones (Barthelmie et al., 2009).Two important wake characteristics are (1) a flow velocity deficit and (2) an increase in turbulence intensity.The former reduces power production of the farm while the latter leads to a higher dynamic loading on downstream turbines but also induces wake recovery.Individual turbine control variables can influence the wake's flow velocity, turbulence intensity and also location.Hence, by changing the control variables of individual turbines, power production of and loading on these controlled turbines and the downwind turbines can be manipulated.Wind farm control aims to find control variables under changing atmospheric conditions such that demanded power production and/or a minimization of the loading can be guaranteed, improving the cost and quality of wind energy.State-of-theart closed-loop dynamic wind farm controllers are based on computationally expensive wind farm models, which make these methods suitable for analysis though unsuitable for online control.The latter is important because it allows for model adaptation to the time-varying atmospheric conditions using supervisory control and data acquisition (SCADA) measurements.As a consequence, more reliable control settings can be evaluated.A survey on wind farm control can be found in Knudsen et al. (2015) and Boersma et al. (2017), for example.In the latter, a clear distinction is made between model-based and model-free control algorithms.This paper is focussed on the former in which it is assumed that controllers are based on a mathematical model of the system.Consequently, the controller performance depends highly on the model quality.Modelling is therefore a crucial step towards successful implementation of model-based wind farm control.
Overviews on wind farm models can be found in Crespo et al. (1999), Vermeer et al. (2003), Sanderse (2009), Sanderse et al. (2011), Annoni et al. (2014), Göçmen at al. (2016) and Boersma et al. (2017).The spectrum of these models ranges from low fidelity to high fidelity.The latter tries to capture relatively precise wind farm flow and turbine dynamics, while the former tries to capture only the dominant characteristics (dynamic or static) in a wind farm.Examples of high-fidelity wind farm models are Simulator fOr Wind Farm Applications (SOWFA) (Churchfield et al., 2012), UTD Wind Farm (UTDWF) (Martinez-Tossas et al., 2014), SP-Wind (Meyers and Meneveau, 2010) and PArallelized LES Model (PALM) (Maronga et al., 2015).These three dimensional (3-D) high-fidelity wind farm models can easily have 10 6 or more states.The resulting computation time can be of the order of days or weeks using distributed computation for simulation times less than the computation time.In other words, the computation time needed for largeeddy simulations (LESs) is in general more than the total time that is simulated.Clearly, these types of models are not applicable for online model-based control.Rather, these models serve as analysis or validation tools.
One way to reduce the high complexity of wake modelling is by using two-dimensional (2-D) heuristic models that only capture specific wake and turbine characteristics in a wind farm in the horizontal plane at hub height.These types of models are found on the low-fidelity side of the spectrum.Most of these wake models exclusively estimate a steady state situation for given atmospheric conditions.Examples of static models are the Frandsen model (Frandsen et al., 2006) and the Jensen-Park model (Jensen, 1983;Katic et al., 1986).One extension of the Jensen model resulted in the parametric model called FLOw Redirection and Induction in Steady-state (FLORIS) (Gebraad et al., 2014b).Two examples of low-fidelity dynamic models are SimWind-Farm (Grunnet et al., 2010) and the model used in Shapiro et al. (2017a), in which relatively simple approximations of the flow deficit are computed using heuristic expressions.
Medium-fidelity models can be found in the middle of the spectrum as they trade off the accuracy of high-fidelity models with the computational complexity of low-fidelity models.These are in general based on simplified versions of the Navier-Stokes equations.For example, in the 2-D dynamic wake meandering (DWM) model (Larsen et al., 2007), assumptions are made regarding the thin shear layer such that the Navier-Stokes equations can be approximated using less computational effort.The authors in Trabucchi et al. (2016) present a model, which is also based on the thin shear layer approximation, but according to the authors it is applicable for non-axisymmetric wind turbine wakes.Wake-Farm (also referred to as FarmFlow) simulates the wind turbine wakes by solving the steady parabolized Navier-Stokes equations in three dimensions (Crespo et al., 1988;Özdemir et al., 2013).Other wind farm models based on the 3-D Reynolds-averaged Navier-Stokes (RANS) equations are Avila et al. (2013) and van der Laan et al. (2015).In Annoni and Seiler (2015), time averaging is applied to the Navier-Stokes equations, resulting in the 2-D RANS equations.The number of states is then reduced by employing a state reduction technique.
Also considered as medium-fidelity models are the ones presented in Boersma et al. (2016b) and Soleimanzadeh et al. (2014).These wind farm models are based on the discretized 2-D Navier-Stokes equations.However, these models do not contain a turbulence model that allows for wake recovery.In addition, these 2-D models do not take any neglected 3-D effects into account and no yaw actuation of the individual turbines is included.
In this paper, a model will be presented that can be considered as a building block for the closed-loop control framework as illustrated in Fig. 1.
In current practice, signals such as power can be measured from a wind farm, but current research is also focussing on estimating wake characteristics using a lidar device (Raach et al., 2017).These and other wind farm measurements are called SCADA data and can be used by an estimator that is able to adapt the model parameters to current atmospheric conditions and/or estimate the full state space, e.g.all the flow velocities at hub height in the farm.The work presented in Doekemeijer et al. (2016) illustrates the latter and employs the dynamic wind farm model presented in this paper.Subsequently, the estimation can then be used to compute optimal control variables using a model predictive controller.The work presented in Vali et al. (2016) illustrates the application of such a model predictive wind farm controller using the dynamic model presented in this work.
The online closed-loop control paradigm as depicted in Fig. 1 demands a control-oriented dynamic wind farm model that will be presented in this paper.Characteristics of such control-oriented models are 1.low computational cost such that online model update, state estimation and control signal evaluation is possible; 2. their dynamical nature such that they can deal with varying atmospheric conditions within relatively small timescales.
The dynamic control-oriented wind farm model presented in this paper, referred to as WindFarmSimulator (WFSim), is applicable in the framework discussed above and satisfies the two points above.It is based on corrected 2-D Navier-Stokes equations and contains a heuristic turbulence model.The Navier-Stokes equations are modified in order to partially correct for the neglected vertical dimension.Each turbine is modelled using the actuator disk model (ADM) and features yaw and axial induction actuation.An important model fea- ture is the exploitation of the sparse system matrices, leading to computational efficiency.WFSim will be compared to high-fidelity flow data and used in a practical control application.
The remainder of this paper is organized as follows.In Sect.2, the mathematical background of the medium-fidelity wind farm model will be explained including a discussion on the rotor and turbulence model.This section ends with an analysis regarding the wind farm model's computation time.In Sect.3, WFSim will be validated in two cases using flow velocities in the longitudinal and lateral directions at hub height and turbine power signals computed with two different LES-based wind farm models.The first case considered is a two-turbine wind farm with turbines modelled using the ADM.The second case is a nine-turbine wind farm with turbines modelled employing Fatigue, Aerodynamics, Structures and Turbulence (FAST) (Jonkman and Buhl, 2005).This paper is concluded in Sect. 4.

Formulation of a dynamic control-oriented wind farm model
In the current section, a simplified wind farm model is formulated that is sufficiently fast for online control but retains some of the elemental features of three-dimensional turbulent flows.In order for the model to be fast, we envisage a 2-D-like model, but adapted to account for three-dimensional flow relaxation.We will dub the resulting model WFSim (WindFarmSimulator).
As a starting point we use the standard incompressible three-dimensional filtered Navier-Stokes equations, as used in LES, i.e. ( The velocity field v = ( v 1 , v 2 , v 3 ) T and pressure field p represent filtered variables; ∇ = (∂/∂x, ∂/∂y, ∂/∂z) T ; the air density ρ, which is assumed to be constant; and τ M represents the subgrid-scale model, which will be defined in Sect.2.1.As is common in LESs of high-Reynolds-number atmospheric simulations with grid resolutions in the metre range, direct effects of viscous stresses on the filtered fields are negligible, so that these terms are left out.Finally, the term f represents the effect of turbines on the flow, as further detailed in Sect.2.2.Although LES filters are usually implicitly tied to the LES grid and filter length scale in the subgrid-scale model, we presume here that v corresponds to a top-hat filtered velocity field, with filter width D, where D is the turbine diameter.Thus, Therefore, we focus on formulating a 2-D-like set of equations for v(x, y, z h ).To this end, we define and w = v 3 (x, y, z h ) and p = p(x, y, z h )/ρ.Moreover, we assume that w ≈ 0, so that the LES equations given in Eq. ( 1) can be reformulated in terms of u as with ∇ H = (∂/∂x, ∂/∂y) T , τ H a 2-D tensor containing the horizontal components of the subgrid stresses τ M , and e 1 and e 2 the unit vectors in the x and y directions, respectively.Finally, we further simplify the equations above using two additional assumptions.First of all, we presume ∂w/∂z ≈ ∂v/∂y.When centred at the turbine axis, this is one of the conditions required for axial symmetry, though axial symmetry also requires further conditions on ∂v/∂z and ∂w/∂y, which are not imposed.In general, ∂w/∂z ≈ ∂v/∂y implies equal divergence-convergence of streamlines in y and z directions.Although this is a very simple condition, we presume it to be good enough to resolve the lack of relaxation of purely 2-D models.If necessary, a more general form (w ≈ 0 and ∂w/∂z ≈ c • ∂v/∂y), with c a tuning parameter (e.g.obtained through state estimation) could be considered, but results in the current work indicate that this may not be necessary.Secondly, we simply neglect the right-hand side of Eq. ( 5).Though this is a rather crude assumption, the rationale is that the modelling term τ H will suffice in the context of a control model, where model coefficients can be updated online based on feedback (see also the discussion in Sect.2.1).Hence our final 2-D-like model corresponds to We emphasize here that the model above is not a classical 2-D model due to the difference in formulation of the continuity equation.In contrast to a standard 2-D model, this allows for flow relaxation in the third direction when, for example, encountering slow down by a wind turbine.This can be seen in Fig. 2 be discussed in more detail in Sect.3.2.1.Here we depict the normalized flow deficit in the wake at 5D downstream of the turbine along the cross-stream axis.The figure illustrates that the standard 2-D Navier-Stokes equations lead to a significant speed up at the wake edges.This is a result from conservation of mass in two dimensions and the flow deceleration in front of the turbine, pushing part of the air around the turbine.In the WFSim model, this speed up is smaller, as mass can also flow around the turbine in the third dimension.In Fig. 2 it can be seen that LES data are better estimated when imposing flow relaxation in the third dimension.Finally, note that partially modelling the missing vertical dimension as proposed above is novel with respect to the work presented in Boersma et al. (2016a).This section is further organized as follows.First, in Sect.2.1, the subgrid-scale model will be introduced.Then, in Sect.2.2, the turbine model will be explained.The discretization of the equations is presented in Sect.2.4, and boundary and initial conditions are discussed in Sect.2.5.

Turbulence model
In the literature, many subgrid-scale models are documented, and to date, model accuracy remains a challenge in LES research (see e.g.Sagaut, 2006).However, in the current paper, an important factor in the selection of a model is simplicity and computational efficiency, rather than accuracy.In fact, in contrast to conventional modelling, in a controloriented model completeness of the turbulence model is not a major issue since unknown model coefficients can be calibrated online using measurements and feedback (Shapiro et al., 2017b), thus also controlling the overall error.Therefore, in this work we fall back to one of the simplest and first known turbulence models, Prandtl's mixing length model.We formulate the stress tensor τ H using an eddy-viscosity assumption, i.e.
with S = 1 2 (∇ H u + (∇ H u) T ) the 2-D rate-of-strain tensor, and ν t the eddy viscosity.The latter is further modelled as in Prandtl (1925): where l u (x, y) is the mixing length.It could be interesting to define the mixing length for each position in the wind farm separately, but this will lead to too many tuning variables.Moreover, in Iungo et al. (2015), the authors illustrate that in a turbine's near wake the mixing length is roughly invariant for different downstream locations, but in the far wake, the mixing length increases linearly with downstream distance.
We use this to formulate a simple heuristic parametrization for the mixing length model so that the number of decision variables will be reduced drastically.From now on we assume that the wind is coming from the east, but can have a direction defined by ϕ.Then, the wind farm will be divided in segments as illustrated in Fig. 3.Each segments has its own (x n , y n ) coordinate system located in the global (x, y) coordinate system.Now we propose the following mixing length parametrization: with G(x, y) a (smoothing) pillbox filter with radius 3, * the 2-D spatial convolution operator, X = {x : x n } and ϕ defined as the mean wind direction (see Fig. 3), which we bound by |ϕ| ≤ 45 • in this work.In addition we constraint d by cos(ϕ)d ≤ |x q −x n |, with x n a turbine x coordinate and x q its downwind turbines x coordinate.We can see l n u (x n , y n ) as the local mixing length that belongs to turbine n and denote it as with and tuning parameter l s that defines the slope of the (linearly increasing) local mixing length parameter.In fact, this parameter could be related to turbulence intensity, i.e., the amount of wake recovery.In this work we will not investigate this relation further.With the formulation above, the number of tuning variables that belong to the turbulence model (l s , d, d ) is reduced to 3ℵ.Additionally, we assume that l s , d and d are equal for each turbine in the farm, which reduces the number of tuning variables that belong to the turbulence model to three, a quantity that could be dealt with by an online estimator.However, in order to have only three tuning variables, the included turbulence model is defined as a simplified mixing length model found heuristically using and adapting information from Iungo et al. (2015).

Turbine model
Turbines are modelled using a classical non-rotating actuator disk model (ADM).In this method, each wind turbine is represented by a uniformly distributed force acting on the grid points where the rotor disk is located.Figure 4 depicts a schematic top-view representation of a turbine with yaw angle γ .
Using such an approach, the force exerted by the turbines can be expressed as www.wind-energ-sci.net/3/75/2018/Wind Energ.Sci., 3, 75-95, 2018  with s = (x, y) T , H[•] the Heaviside function and δ[•] the Dirac delta function, e ⊥,n the unit vector perpendicular to the nth rotor disk with position t n .Furthermore, we have C T n the disk-based thrust coefficient following Meyers and Meneveau (2010), which can be expressed in terms of the classical thrust coefficient C T n using the following relation: with a n the axial induction factor of the nth turbine.Interestingly, the coefficient C T n can be directly related to the turbine set point in terms of blade pitch angle and rotational speed (see, e.g.Appendix A in Goit and Meyers, 2015).In the WFSim model, C T n and yaw angle γ n are considered as the control variables and can thus be used to regulate the wakes and hence wind farm performance.Furthermore, the scalar c f in Eq. ( 13) can be regarded as a tuning variable and will in this work be set equal for all turbines in the farm.

Power
From the resolved flow velocity components, the power generated by the farm is computed as It is stated in Goit and Meyers (2015) (Appendix A) that when there is no drag and swirl is added to the wake, C T n = C P n .Since this is an idealized situation, a loss factor will be introduced such that C P n = c p C T n .The scalar c p can be seen as a tuning variable and will be set equal for all turbines in the farm.In the power expression above, we have the factor cos(γ n ) 3 with exponent 3.In literature such as Gebraad et al. (2014a) and Medici (2005, p. 37), for example, numerical values for the exponent were given according to LES and wind tunnel data, respectively.However, to date, the exact value for it is still under research and since this is outside the scope of this study, the value of the exponent will be 3.This concludes the formulation of the WFSim model.In order to resolve for flow velocity components and wind farm power, the governing equations given in Eqs. ( 7) and ( 8) need to be discretized, a topic that will be discussed in the following section.

Discretization
The set of equations are spatially discretized over a staggered grid following Versteeg and Malalasekera (2007).This is carried out by employing the finite volume method and the hybrid differencing scheme.Temporal discretization is performed using the implicit method that is unconditionally stable (Versteeg and Malalasekera, 2007).This comes down to deriving the integrals: with V the volume of one cell (see Fig. 5) and t the sample period.One obtains, for each cell, the following fully discretized Navier-Stokes equations (for detailed derivation we refer to Appendix A): x momentum equation for the (i, J )th cell (black in Fig. 5): y momentum equation for the (I, j )th cell (yellow in Fig. 5): continuity equation for the (I, J )th cell (pink in Fig. 5): • and the forcing terms f • •,• depend on the state at time k.Detailed definitions of these coefficients are given in Appendix A, Table A1.Note the appearance of the previously explained factor 2 (see Eq. 8) in Eq. ( 18).
Next, the state vectors u k , v k , and p k and control variable vectors ν k and γ k at time step k will be defined: with N x and N y the number of cells in the x and y directions, respectively, and ℵ the number of turbines in the wind farm.
Each component in u k , v k and p k represents a flow velocity and pressure at a point in the field defined by the subscript.
For clarity reasons, an example of a staggered grid is depicted in Fig. 6.

Boundary and initial conditions
All the components that are not contained in the vector, u k , v k and p k , but do appear in the staggered grid need to be defined.For the flow velocity components, first-order conditions on the west side of the grid are prescribed assuming the wind is coming from the east.These Dirichlet inflow boundary conditions are related to the ambient inflow defined as u b and v b and can vary over time.Zero stress (also referred to as Neumann) boundary conditions are prescribed on the other boundaries.Therefore, for the flow velocity components on the boundaries we define for j = 3, 4, . .., N y − 1.
For the initial conditions, we define all longitudinal and latitudinal flow velocity components in the field as u b and v b , respectively, which are the boundary velocity values.The initial pressure field is set to zero.Note that by defining the boundary conditions as given above, the assumption is that the wind is coming from the east in Fig. 6, which coincides with the definition of the mixing length (see Sect. 2.1).Finally, the equations given in Eqs. ( 7) and ( 8) can be transformed to the difference algebraic equation1 : with containing all flow velocities in the longitudinal and lateral S. Boersma et al.: WFSim Δy 2,3 δy 4,5 Figure 6.Example of a staggered grid with cells each having volume V .In WFSim, the grid is quadrilateral.
Table 1.Mean computation time per simulation time step t cpu vs. number of states n q for the WFSim representation as given in Eq. ( 21).
Computations are performed on a regular notebook with one core.
By defining N x , N y , x I,I +1 , y J,J +1 , and the turbine positions, a wind farm topology is determined.Next, ambient flow conditions u b and v b , tuning parameters c f , c p , d, d , l s and the control variable w k need to be specified.The system given in Eq. ( 20) is then fully defined and can be solved.

Computation time
When discretizing partial differential equations, a trade-off has to be made between the computation time and grid resolution.Typically, a higher resolution results in more precise computation of the variables but also increasing computation time.In WFSim, computational cost is reduced by exploiting sparsity and by applying the reverse Cuthill-McKee algorithm (George and Liu, 1981). 2 The latter is applicable due to the fact that the matrix structure is fixed.The inter-  vs. number of states n q .The red dashed line is WFSim as presented in Eq. ( 20) and blue is WFSim as presented in Eq. ( 21).Note that the number of cells is approximately n q /3 with n q the number of states.
ested reader is referred to Doekemeijer et al. (2016) for more information on the Cuthill-McKee algorithm in WFSim.
In this section, the mean computation time needed for one time step t cpu will be analysed.The presented results are obtained on a regular notebook with an Intel Core i7-4600U 2.7 GHz processor employing one core and MATLAB.Since the objective is to do online control, i.e., it is desired to reduce computational complexity, this section introduces a second WFSim representation.The first representation was given in Eq. ( 20) while the second is defined as The difference can be found in the descriptor matrix.In the representation above, the elements A xy (u k ), A yx (u k ) that occur in Eq. ( 20) are set to zero.This can be justified by the fact that their contribution is negligible since these matrices contain elements that, for our case studies, are of the order of O(1) while the elements in A x (u k , v k ) and A y (u k , v k ) are of the order of O(3).Consequently, no significant change in the computed flow field has been observed, but a decrease in t cpu has (see Fig. 7).Therefore, the remainder of this paper will continue with the WFSim representation given in Eq. ( 21).Table 1 depicts more numerical values of t cpu for this WFSim representation.
From Table 1 we can conclude that t cpu increases between quadratic and linear with respect to the number of states n q for n q < 221 × 10 3 .It depends on the computer properties how much you can increase the number of states until the CPU is out of memory.

Simulation results
In this section, WFSim flow and power data will be compared with LES data and it is organized as follows.In Sect.3.1, quality measures are introduced.In Sect.3.2.1,WFSim data are compared with PALM data, and in Sect.3.2.2WFSim is validated against SOWFA data.In both simulation cases, the thrust coefficients C T are varied while the yaw angles are set to zero.

Quality measures
Suppose we have at time k a measurement of one quantity z k ∈ R N and its estimation ẑk ∈ R N .Define the prediction error e k = ẑk − z k .The quality measure RMSE is, for time step k, defined as This measure is used to compare the flow centreline velocity U c k (x) and power signals from LES and WFSim data for different model parameters.The flow centreline is, for one time step, defined as the laterally averaged longitudinal flow velocity throughout the simulation domain across the rotor diameter.Mathematically this can, for LES data at time step k at longitudinal position x i , be defined as with y s the y coordinate of one cell across the line y ⊂ y, which contains N y number of cells and has a length equal to the rotor diameter.From WFSim data, the flow velocity component at the rotor centre will be taken across the position x.In this work we compare lateral and longitudinal flow velocity components at hub height and power signals calculated with LES with lateral and longitudinal flow velocity components and power signals calculated with WFSim.3

Axial induction actuation
Studies such as Shapiro et al. (2017a), Munters and Meyers (2017), Vali et al. (2017) and van Wingerden et al. (2017) illustrate that axial induction actuation can be used in active power control where the objective is to provide grid facilities.In order to utilize the WFSim model in active power control, it is important to first validate it when exciting the thrust coefficient.
In the following, WFSim is compared with simulation data from PALM (Maronga et al., 2015) and SOWFA (Churchfield et al., 2012), both high-fidelity wind farm models that were briefly discussed in Sect. 1.The latter includes the actuator line model (ALM) while the former employs the ADM.4

PArallelized LES Model (PALM) and WFSim
PALM predicts the 3-D flow velocity vectors and turbine power signals in a wind farm using LES and is based on the 3-D incompressible Navier-Stokes equations. 5Table 2 gives a summary of the two-turbine wind farm simulated in WFSim.A summary of the PALM simulation set-up can be found in Appendix B. The applied control signals are depicted in Fig. 8 and are chosen such that different system dynamics are excited.The final values for the tuning parameters are obtained using a grid search.Figures 9 and 10 show a comparison of the mean flow centreline and the wind farm power, respectively.A flow field evaluated with both the WF-Sim model and PALM can be found in Appendix B.
In Fig. 9, the mean flow centrelines through the farms of WFSim and PALM are relatively similar.The PALM data exhibit more turbulent fluctuations due to the presence of a more sophisticated turbulence model, which allows for better capturing small-scale dynamics such as turbine-induced turbulence.However, the WFSim model is capable of estimating wake recovery similar to the PALM model.The recovery in the WFSim model is due to the turbulence model as presented in Sect.2.1.It is in fact the slope of the local mixing length parameters that can determine the amount of wake recovery, or more precisely, the larger this slope, the higher the wake recovery.It is therefore interesting to link this tuning variable to the turbulence intensity in the farm.Furthermore, it can be seen in Fig. 10  is capable of estimating the wind farm power.Since both the WFSim model and PALM employ the ADM, fast fluctuations in the power signal can be observed.This is due to the lack of rotor inertia in both simulation cases.The simulation case presented in this section illustrates that the WFSim model, in which the third dimension is partially neglected, is able to estimate wind farm flow and power signals computed with a 3-D LES wind farm model.In Sect.3.2.2, a SOWFA case study will be presented, a LES model that includes turbine dynamics.

Simulator fOr Wind Farm Applications (SOWFA) and WFSim
SOWFA predicts the 3-D flow velocity vectors in a wind farm using LES and is based on the 3-D incompressible Navier-Stokes equations.For turbine modelling it employs the ALM, which is a more sophisticated model than the ADM (Sanderse et al., 2011).In addition, the FAST model from the National Renewable Energy Laboratory (NREL) is implemented (Jonkman and Buhl, 2005).This model calculates turbine power production, blade forces on the flow and structural loading on the turbine.In the SOWFA sim-ulation presented, the NREL 5 MW wind turbine is simulated (Jonkman et al., 2009).
The SOWFA data set used in this work for validation is equivalent to the set used in van Wingerden et al. (2017).The thrust coefficient C T is not a control variable in SOWFA due to the employment of the ALM and therefore has to be estimated.This will be discussed in the following paragraph.

Turbine operating settings
Since F sowfa k , U sowfa k and ρ can be obtained from SOWFA data and the yaw angles are given, all the variables in Eq. ( 24) www.wind-energ-sci.net/3/75/2018/Wind Energ.Sci., 3, 75-95, 2018   are known; hence the control variable C T can be estimated from SOWFA data for each turbine.6It will be used, together with the yaw angle, as an input to the WFSim model.
In the following, flow data at hub height from a nineturbine SOWFA simulation case will be compared with WF-Sim data.See Fig. 12a for the simulated wind farm topology.The turbines are excited with thrust coefficients as depicted in Fig. 11.These excitation signals are estimated from SOWFA data using the relation defined in Eq. ( 24).Table 3 presents the WFSim parameters used during simulations.The tuning variables of the WFSim model are found using a grid search and the inflow conditions u b , v b are estimated from SOWFA data.
Figures 12b and 13 depict a mean flow centreline (see Eq. 23) comparison for each row at four time instances.It can be concluded that the mean flow centreline derived from WFSim data approximates the mean flow centreline derived from SOWFA data.In Fig. 14, time series of the power signals from SOWFA and WFSim are depicted.The signals from the latter are more oscillating than the power signals from SOWFA.This is due to the fact that the power expression in WFSim is a non-linear static map depending on the C T .Thus, no turbine dynamics are taken into account, which is contrary to SOWFA in which the FAST turbine model is simulated.However, important characteristics can be cap-tured with WFSim.A flow field evaluated with both the WF-Sim model and SOWFA can be found in Appendix C.
WFSim is capable of estimating dominant wake dynamics, the objective of the control-oriented model WFSim.Smallerscale and stochastic effects can be measured by sensors and incorporated using an estimator based on WFSim, as has been shown in Doekemeijer et al. (2016Doekemeijer et al. ( , 2017)).

Conclusions
Current literature on wind farm control can be categorized into model-free and model-based methods.This paper focused on the latter category.Here, a distinction can be made between type of model employed, a steady-state or dynamic wind farm model.In order to use the closed-loop control paradigm, and account for model uncertainties, we think it is important to utilize a dynamic wind farm model for controller design and possible online wind farm control.In this paper, such a control-oriented dynamic wind farm model, referred to as WFSim, has been presented. 7It is a wind farm model that can predict flow fields and power production and includes turbines that are modelled using actuator disk theory and is based on modified two-dimensional Navier-Stokes equations.Completely neglecting the third (vertical) dimension is a too crude assumption to describe the flow in a wind farm accurately enough for control purposes.In this paper, we included a correction term in the continuity equation.It has been illustrated that the inclusion of this factor reduces the effect of neglecting the third (vertical) dimension.More precisely, it has been shown that the speed-up effect of the flow on the right and left downwind of a turbine will be reduced when solving for the corrected Navier-Stokes equations compared to the standard two-dimensional Navier-Stokes equations.It has been shown that this resulted in a better approximation of LES data.
In addition, a turbulence model was included, taking into account the desired wake recovery.The heuristically found turbulence model is based on Prandtl's mixing length hypotheses, where the mixing length parameter is made de-pendent on the downstream distance from the turbine rotors and also dependent on the mean wind direction.After theoretically formulating the WFSim model, this paper followed by illustrating that the computed flow velocities and power signals from the 2-D-like WFSim model can estimate flow velocity data and power signals from the 3-D high-fidelity wind farm models PALM and SOWFA.The necessary computation time of the WFSim model is a fraction of what is needed to perform LES, making it suitable for online control.This work focussed on axial induction actuation, but future work will also include the validation of yaw actuation and wind direction changes.For the simulation cases presented, no grid convergence studies have been performed, but future work should entail this.In addition, future work will entail the online update of the tuning variables c f , c p , d, d , l s by an observer and the employment of the presented dynamic wind farm model in an online closed-loop control scheme.This section will present the necessary derivations to go from Eqs. ( 15) to (20), i.e. it will elaborate on the discretization of the NS equations.In the following subsections, all terms in the NS equations will subsequently be dealt with.
A1 Discretizing the convection (non-linear) terms The non-linear term that occurs in the momentum equations can be spatially discretized by deriving where (u 2 δy) e , (u 2 δy) w are the quantities u 2 at the east and west sides of the cell with surface δy e , δy w , respectively.Similarly, (uv x) n , (uv x) s are the quantities uv at the north and south sides of the cell with surface x n , x s , respectively.Assuming δy = δy e = δy w and x = x n = x s , the above can be written as F ex = ρu e δy, F wx = ρu w δy, F nx = ρv n x, and F sx = ρv s x.In Versteeg and Malalasekera (2007) this is referred to as a convective mass flux approximation.The above can then be written as In Fig. 5 we observe that u e , u w , u n , u s , v n , v s are not defined for the black cell.Applying central differencing approximates the terms as follows: We can now write In Eq. (A1), central differencing is applied.A disadvantage of this method is that it does not use prior knowledge on the flow direction.The upwind differencing scheme, however, employs this prior knowledge as explained in Versteeg and Malalasekera (2007).A combination of the central and upwind differencing scheme is the hybrid differencing scheme.When applying this, the above can be written as with , 0] and c px i,J = c ex i,J + c wx i,J + c nx i,J + c sx i,J + F ex i,J − F wx i,J + F nx i,J − F sx i,J .In WFSim, the coefficients c • i,J and F • i,J are evaluated for time k while the other flow velocity components are computed for time k + 1.

y momentum equation
Deriving the non-linear term in the y momentum equation yields The intermediate steps are omitted here since they are similar to the steps presented when handling the non-linear term Wind Energ.Sci., 3, 75-95, 2018 www.wind-energ-sci.net/3/75/2018/ in the x momentum equation.Note, however, that the discretization is evaluated using the yellow cell (see Fig. 5).
When applying the hybrid differencing scheme, the above can be written as The pressure components are evaluated for time k + 1.

Evaluate
and δy = δy j,j +1 .Substituting these expressions yields The second term evaluates as x.
Here we have which can be rearranged to with T px i,J = T ex i,J +T wx i,J +T nx i,J +T sx i,J .The coefficients T • i,J will be computed for time k while the flow components will be evaluated for time k + 1.Here we have and x = δx i,i+1 .Substituting these expressions yields V ∂ ∂y l u (x, y) 2 ∂u ∂y ∂v ∂y dV = l u (x I , y J ) 2 (u i+1,J − u i+1,J −1 )δx i,i+1 y J −1,J δy j,j +1 T ny I,j ×(v I,j +1 − v I,j ). . .−l u (x I , y J −1 ) 2 (u i,J − u i,J −1 )δx i,i+1 y J −1,J δy j −1,j T sy I,j (v I,j − v I,j −1 ). (A7) The second term evaluates as V .
Table A1.Fully discretized Navier-Stokes equations and all their coefficients.
All the coefficients derived above are given in Table A1.

Figure 1 .
Figure 1.General dynamical closed-loop control framework with measurements z k and its estimation ẑk and state estimation qk .The signals r k and w k are reference and control signals, respectively.In this paper we present a dynamic model that is compatible with this framework.

Figure 3 .
Figure 3. Schematic illustration of the mixing length.

Figure 5 .
Figure5.One cell for the x momentum equation (grey, u i,J ), one for the y momentum equation (yellow, v I,j ) and one for the continuity equation (pink, p I,J ).All three cells have equal dimensions and overlap.

Figure 7 .
Figure 7. Mean computation time per simulation time step t cpu vs. number of states n q .The red dashed line is WFSim as presented in Eq. (20) and blue is WFSim as presented in Eq. (21).Note that the number of cells is approximately n q /3 with n q the number of states.

Figure 8 .Figure 9 .
Figure 8. Excitation signals for the two-turbine simulation case.The yaw angles are set to zero.

For
estimating the control signals C T n , the turbine's fore-aft bending moment M sowfa k calculated with FAST is exploited.Using the relation M sowfa k = F sowfa k z h with z h the hub height, the (indirect) measured thrust force F sowfa k can be derived.An estimation from SOWFA data of the rotor flow velocity U sowfa k is obtained by averaging the flow velocity components across the rotor.Using the standard ADM yields for each turbine = 0.1 (s) Atmospheric conditions u b = 12, v b = 0 (m s −1 ), ρ = 1.2 (kg m −3 ) Force and power factor c f = 5 2 , c p = 1.1 Turbulence model d = 635, d = 76.2 (m) l s = 0.17

Figure 12 .
Figure 12.Topology simulated wind farm (a) and mean flow centreline at four time instances through the first row (b).

+
Deriving the term in the x momentum equation (first element in the vector above) yields (uv x) n − (uv x) s , uv) n x − (uv) s x .
u e − F wx u w + F nx u n − F sx u s .
. Similar as before, the coefficients c • I,j and F • I,j are evaluated for time k while the other flow velocity components are computed for time k + 1.A2 Discretizing the pressure gradientFor the pressure gradient we evaluate V ∂p ∂x ∂p ∂y dV = p I,J − p I −1,J δy p I,J − p I,J −1 δx .
j − v I −1,j x I −1,I , and y = y J −1,J .Substituting these expressions yields x i+1 , y j) 2 u i,J − u i,J −1 x (u i,J − u i,J−1 ), with T py I,j = T ey I,j + T wy I,j + T ny I,j + T sy I,j .The coefficients T • I,jwill be computed for time k while the flow components will be evaluated for time k + 1.

Table 2 .
Summary of the WFSim simulation set-up.

Table 3 .
Summary of the WFSim simulation set-up.