1. Introduction
Improving wind farm design and control is at the core of wind energy research aimed at facilitating greater wind energy penetration. Design problems focus on maximizing wind farm power production and reducing the levelized cost of energy by optimizing wind farm layouts [
1,
2,
3,
4] and wind turbine set points for yaw, tilt, and thrust [
5,
6]. Active control attempts to actuate turbines dynamically to reduce farm level power fluctuations [
7], track a power output reference signal to provide power grid services [
8,
9,
10,
11], or maximize power production [
12,
13].
Design and control studies rely on engineering models that can be easily computed and optimized because high-fidelity simulations are too computationally expensive to be practical. Traditional engineering design models focus on wake growth rates, the shape of the velocity deficit in the transverse directions downstream of the turbine (the wake profile), and wake interactions. However, the details and capabilities of these models vary considerably. For example, the wake profile shape can be modeled as a top-hat [
14,
15], Gaussian [
16,
17], or a more complicated step-wise function [
18]. Wake areas are often modeled as growing linearly with downstream distance [
14,
15], but have also been represented with other functional forms [
19]. Velocity deficits arising from wake interactions may be computed as a linear sum [
20,
21] or as a sum of squares [
14,
15]. Related modeling approaches include the effects of the atmospheric boundary layer (ABL) on the wind turbine array through one-way [
22] or two-way coupling [
23,
24] with ABL models [
19,
25,
26]. Design-related models may also include certain dynamic effects such as yaw and tilt actuation [
27,
28,
29] or the turbulent meandering of wakes [
30]. These design models are typically either static or account for limited temporal variations and therefore cannot capture changing wind farm conditions or power output levels.
Active and closed-loop control applications, on the other hand, require time-dependent models that can capture the effects of control actions. Dynamic modeling approaches also vary considerably in fidelity, modeling philosophy, and underlying assumptions. The simplest models are discrete-time extensions of steady-state wake models [
31] and continuous-time one-dimensional models [
8,
32]. Higher fidelity, yet more computationally-expensive, models include two-dimensional Reynolds-averaged Navier–Stokes (RANS) simulations [
33] and very coarse large-eddy simulations (LES) [
34].
Design and control-oriented models also differ in their implementation. In design problems, the goal is to estimate the effect of wind turbines on the fluid flow prior to implementation or construction of the farm. Therefore, steady-state model coefficients cannot be directly measured and are instead estimated empirically [
20,
21,
35] based on site-specific data, or more infrequently, via physics-based methods [
24]. In closed-loop control applications, parameter estimation is enabled through real-time data regarding power output or from flow field sensors that provide a means of fitting model coefficients to local site-specific attributes based on the current flow state [
4,
36,
37].
In other ways, dynamic models have a range of complications that are not relevant to steady-state models. For example, the temporal effect of large-scale turbulent structures and the non-linearity of advection affect inter-turbine travel times. These interactions impact the instantaneous power output of individual turbines within a farm, but may not change the time-averaged power output. Therefore, these effects have not been extensively discussed in studies that focus on steady-state applications.
In this paper, we build on our previously-proposed continuous-time one-dimensional dynamic wake model [
32] to provide a unified modeling framework for steady-state design and active control applications. We demonstrate how the selection of modeling parameters, including the wake mixing and expansion model, the wake profile shape, the advection scheme, and the inflow model, can improve model fidelity. For the dynamic model, we introduce a filtering approach to incorporate spanwise inhomogeneity in the inlet conditions to improve the inflow model, as well as state and parameter estimation to improve power output prediction. We then demonstrate that these estimation techniques can correct model errors, such as those introduced when using a linear advection model.
The remainder of this paper is organized as follows. We introduce the model proposed in [
8] in
Section 2, along with a new definition of the wake profile that improves the model accuracy. In
Section 3, we discuss how to parameterize the model for design applications that employ the steady-state version of the model, where improved performance is achieved by linking the wake expansion rates to top-down atmospheric models [
26]. In
Section 4, we apply the dynamic model to a range of problems starting with the prediction of the power trajectory of a wind farm at start-up. We then introduce measurement-based state and parameter estimation approaches to account for spanwise inhomogeneities in the inlet and modeling errors. Finally, we summarize and discuss the conclusions in
Section 5.
2. Wake Model
This paper builds on the dynamic wake model introduced in [
8] and extended to include the effect of yaw in [
32]. The model presented here differs from previous versions in the following ways. Velocity deficits are linearly superposed; the wake growth rate is linked to the friction velocity; and the wake profile model incorporates a smooth transition from a top-hat to a Gaussian wake profile.
We begin by considering a wind farm with N turbines each with diameter . The coordinate system is defined such that x is the streamwise direction, y is the spanwise direction, and z is the vertical direction. The center of the turbine’s rotor is located at . The inflow velocity to the farm may vary in time and space ; however, in some cases, the time-averaged or time- and space-averaged velocity may be the more appropriate, and the arguments of the averaged quantities will be omitted.
The streamwise velocity deficit behind the
turbine
is governed by the one-dimensional partial differential equation (PDE) (no Einstein summation convention over
n is implied):
The advection velocity
can be specified is several ways. Here, we focus on two cases: the linear advection model of the form
, which was used in [
37] for ease of implementation, and the nonlinear form
. The function
is a normalized Gaussian function with characteristic width
R:
The function
defines the wake expansion rate as:
which depends on the wake area
that is defined in
Section 2.1.
The source term
is specified by the initial velocity deficit
, which is obtained from an inviscid model [
14,
38]. In the case where the thrust is pointing in the direction of the upstream velocity, the appropriate inviscid model that relies on actuator disk theory [
38] gives an initial velocity deficit of:
where
is the thrust coefficient and
is the upstream velocity at turbine
n. The initial velocity deficits can also be described in terms of the local thrust coefficient
and the disk-averaged velocity across the turbine rotor
.
For the dynamic model, we follow the convention of LES [
25,
39] and use
and
in the model. Assuming that the resulting velocity field is consistent under Galilean transformations, we adopt a linear superposition of velocity deficits. The resulting velocity field is found by distributing the velocity deficits using the wake profile
defined in
Section 2.2 and subtracting the result from the inlet flow velocity
:
The linear superposition approach simplifies the implementation within a model-based control framework and agrees well with the data, as shown in
Section 3. This disk-averaged velocity for turbine
n,
, is found by sampling the velocity field according to:
where
is the squared radial coordinate of the
turbine.
For the steady-state case, it is easier to use
and
. The corresponding steady-state solution of the dynamic model [
32] is:
Linear superposition of the other wakes can be used to find the upstream velocity of each turbine:
2.1. Turbulent Wake Mixing and Expansion
The wake model in
Section 2 captures the dynamics of wake advection, inter-turbine wake interactions (through velocity deficit superposition), and the effect of wake growth on the magnitude of the velocity deficit in the wake. A functional form for the wake area
, defined using the normalized wake diameter
, that properly captures the effect of turbulent mixing of the wake with the surrounding turbulent atmospheric boundary layer is needed to specify the model fully. We now derive the required form for the wake area through the self-similarity assumption applied to the Reynolds-averaged mean momentum equation associated with a turbine in the ABL. For convenience, we dispense with the turbine indexing in this section, and the origin of the coordinate system is placed at the center of the turbine rotor.
After linearizing the advective term, assuming the pressure gradient vanishes in the region sufficiently downstream from the turbine, neglecting the streamwise Reynolds stress, and using the eddy viscosity model, the Reynolds-averaged mean momentum equation for the streamwise velocity
u in the region away from the turbine forcing becomes:
where
is the eddy viscosity, and we have assumed a constant inflow velocity
for the turbine of interest. The far wake of a turbine has been shown empirically to have a self-similar profile [
16]; therefore, the velocity may be written as:
where
is a self-similar function of
and
is proportional to the diameter of the wake.
For a free wake, the eddy viscosity scales with the wake width and the velocity deficit
[
40]. In the atmospheric boundary layer, however, the appropriate velocity scale is the friction velocity
, giving the following eddy viscosity:
Substituting (
10) and (
11) into (
9) and expanding the derivatives, we obtain:
where primes denote derivatives with respect to
. Rearranging (
12) in dimensionless form leads to:
The expression (
13) shows that in order for the profile to be self-similar, the following two conditions must hold in the wake:
Integrating (14), we find that
, and the wake diameter expands linearly with downstream distance, as in the Jensen model [
14,
15]. This yields a normalized wake diameter of:
where the wake expansion coefficient:
is proportional to the ratio of the friction velocity and hub-height velocity, validating the empirical suggestion of Frandsen [
35]. The parameter
will remain as the only free parameter available to the model.
While the solution to (
15) only requires
, we can find the scaling using momentum flux conservation in the far field limit using:
where
C is a constant that depends on the initial velocity deficit in the wake
and the integral
. Therefore, the velocity deficit scales inversely with the square of the downstream distance
. Solving the steady-state version of the PDE (
1) also provides this scaling.
While the linear growth rate is valid in the wake region of the turbine, it is ill-posed upstream of the turbine, where it can vanish or become negative. Therefore, we instead use the following modified function [
8] that smoothly approximates the linear expansion in the far wake while avoiding becoming less than unity close to the turbine:
Defined this way, the normalized wake diameter specifies the wake expansion rate (
3) needed to define the velocity deficit PDE (
1) fully.
2.2. Wake Profile
We now propose a new functional form for the wake profile
in (
5). Experimental observations have revealed that the shape of the wake profile as a function of spanwise distance slowly transitions from a top-hat profile (cf. the left panel of
Figure 1) immediately following the rotor disk to a self-similar Gaussian profile farther downstream in the far wake [
16,
20,
28]. For simplicity, the Jensen model traditionally used top-hat profiles for the entire downstream evolution of the wake [
14,
15], but such profiles poorly describe the far wake profile that dominates the interaction of upstream wakes on downstream turbines. Other models use the Gaussian profile for the entire wake [
16], which overestimates the centerline deficit near the turbine, or model the transition from the top-hat as a decaying potential core [
28].
Here, we propose a super-Gaussian function that more accurately approximates the known transition from a top-hat profile to a Gaussian profile. As in
Section 2.1, we remove the turbine indexing for convenience and set the origin of the coordinate at the center of the turbine rotor. The proposed wake profile function is of the form:
where
is a parameter defining the width of the wake relative to the wake diameter
and
is a parameter that characterizes the smoothness of the wake. In order to transition from a top-hat profile at the rotor disk to a Gaussian in the far wake, we require that:
Furthermore,
should be specified such that the profile becomes nearly Gaussian at
[
28]. To satisfy these conditions, we define:
where
at
.
The wake function must also be defined to conserve mass, which means that for a single turbine subject to the inflow velocity
:
Alternatively, the integral can be expressed as:
the solution of which [
41] (p. 342 #3.478.1) is:
Therefore, the scaling value
is:
with
and
.
The resulting profiles for the velocity in the wake, assuming
, and
, are shown in
Figure 1. These profiles begin as a top-hat with the expected velocity of
, at the rotor disk, where
a is the induction factor. Around
, the velocity is a smoothed top-hat with
in the middle of the wake. Further downstream, the deficit decays, transitioning to a Gaussian by mimicking the potential core region noted in [
28].
This functional form captures the mean velocity in both the near wake, where the velocity deficit of an actuator disk has a top-hat profile, and the far wake, where the velocity deficit has a Gaussian profile. Small-scale features, such as trailing root and tip vortices, and turbine-specific thrust profiles are not captured by this functional form. In cases where these specifics are deemed important, the functional form can only be assumed to be valid for . With the wake profile specified, the model presented in this section is fully defined with a single free parameter that specifies the wake recovery rate inside the farm.
3. The Steady-State Model for Design Applications
For design problems, we require a method to specify the unknown coefficient
for the wake expansion coefficient
(17). In typical models, it is obtained using empirical tables based on surface and atmospheric conditions. The recently-proposed coupled wake boundary layer (CWBL) model [
23,
24,
42] eliminates the need for empirical tables by coupling the wake expansion coefficient of the Jensen model to the top-down model [
19,
25], which provides information about the turbulence characteristics of the wind turbine array boundary layer.
The CWBL determines the single free parameter of the Jensen wake model, the fully-developed wake expansion coefficient , by matching the power prediction of the top-down and Jensen wake model. For a given wake expansion coefficient, the top-down wake model planform thrust coefficient is determined by finding the wake area of the farm; i.e. , where is the planform area of the farm and is the effective wake area coverage. The wake expansion coefficient in the developing region is then assumed to transition from to using a heuristic equation.
Here, we use a related approach, but make several important modifications. First, we replace the Jensen wake model with the steady-state solution of the wake model described above. The remaining differences in the approaches are highlighted in the remainder of this section. Most notably, we dispense with the rather strong assumption of the standard top-down model that the disk-averaged velocity is equal to the planar-averaged hub-height velocity
to calculate the planform-averaged thrust coefficient and power production. We represent entrance effects using the developing top-down model [
26,
43] and replace the heuristic wake growth parameter development with the classic equation for the developing internal boundary layer.
We detail our approach by first introducing the fully-developed top-down model [
25]. This model assumes constant stress regions above and below the wind turbine layer, with friction velocities of
and
and roughness heights of
and
, respectively. The planar-averaged momentum balance results in the following relationship:
where
is the planar-averaged velocity at hub height
. Within the wind turbine wake region
, where
R is the rotor radius, the wind turbine wakes are assumed to generate an additional eddy viscosity
. The resulting eddy viscosity in the wake region is
below the hub height
and
above the hub height
, where
.
Applying momentum conservation and velocity continuity across the layers, the roughness height of the wind farm is:
where
. The friction velocities are:
where
is the free-stream friction velocity and
is the height of the boundary layer. The planar-averaged hub-height velocity is:
In the CWBL model, the disk-averaged velocity is then assumed to be approximately equal to the planar-averaged hub-height velocity
to calculate the power production of the farm.
In contrast to the CWBL model, we assume that the top-down model only provides useful predictions of the planar-average hub-height velocity, but not of the individual turbine power, which depends on the more local velocity at each turbine rather than the planar-averaged velocity. As a result, the proper matching condition to find the free parameter
of the present model, which scales the wake expansion coefficients as a function of friction velocity, is to require consistency between the planar averaged velocities from the top-down and the wake model. In practice, we find
that minimizes the difference between the modeled planar-averaged velocities:
where the superscripts
and
correspond to the top-down and wake models, respectively.
In order to match the planar-averaged velocity from the wake model
corresponding to each turbine, the wind farm’s planform area must be divided into regions that “belong” to each turbine. Here, we propose to use Voronoi tessellation to identify these areas. Difficulties arise at the edge of the farm, where the edge turbines are closest to infinite sections spanning the entire planform area. We address this issue by performing an initial Voronoi tessellation and reflecting each edge turbine across the adjacent line segments of the tessellation to create ghost points. Using the ghost points, Voronoi tessellation is used again to find the planform area of each turbine. The area sections obtained through this algorithm are shown in
Figure 2 and
Figure 3 for the Horns Rev wind farm and a randomly-generated wind farm, respectively. The left panels show the original Voronoi tessellation with the ghost point projections. The right panels show the final tessellation of the wind farm.
We model entrance effects using the developing top-down model [
26,
43], where the boundary transitions from a single log-law profile with surface roughness height
to the fully-developed top-down model described above. In the transition region, an internal boundary layer grows from the height of the wind farm at
until it reaches the final boundary layer height at
. The classical internal boundary layer growth model dictates that:
until
. The friction velocities (30) and (31) and planar-averaged velocity (32) evolve with streamwise distance
x, where
is replaced with the internal boundary layer height
.
In our application of the model, the planform thrust coefficient needed for the top-down model is calculated using the force determined by the wake model:
where
is the planform area of the Voronoi cell around the
turbine. This planform thrust coefficient also provides the surface roughness
. Although this value should strictly represent the roughness height of the fully-developed regime, in smaller wind farms, it is difficult to find the roughness height from matching. Therefore, we assume that within each cell, the planform thrust coefficient represents the average coefficient from upstream turbines. Furthermore, the starting location of the internal boundary layer is assumed to coincide with the most upstream turbine, the wake of which, defined by the cone with diameter
, intersects the rotor of the turbine.
Using the friction velocities at each turbine, we can determine the wake expansion coefficients that the effective friction velocity at each turbine is the average of the low and high friction velocities
, and the advection velocity is the planar-averaged velocity
. As a result, the wake expansion coefficient (17), which is the input to the wake model, can be written as:
For a single turbine, this reduces to the expected
. This improvement replaces the heuristic development function with a physically-based function and eliminates the difficulty of identifying a fully-developed regime to determine the effective wake area coverage that is needed in the CWBL approach.
We validate our proposed steady state model by applying it to the Horns Rev wind farm and compared the results to LES, as well as the Jensen and CWBL models. The Horns Rev wind farm is composed of 10 rows of eight turbines with rotor diameters of
m and hub heights of
m. The turbines are arranged such that the spacings are
and
in the streamwise and spanwise directions when the wind direction is at
. A friction velocity of
m/s, boundary layer height of 500 m, and surface roughness of
m are selected to match the LES data. The thrust coefficient is
(
). Details about the LES setup and CWBL and Jensen model implementations can be found in [
24,
44].
The wind farm power production predicted by the model is compared to the LES, Jensen, and CWBL models in
Figure 4. The current model shows good agreement with the LES data for all angles, with particularly good agreement for the
wind angle and in the early rows of the farm for all cases. The results are comparable to the CWBL model, but our approach outperforms the Jensen model.
Figure 5 illustrates that the flow field obtained is more realistic than that generated by models using simple top-hat or Gaussian profiles (the latter are shown in [
21]), due to the use of the wake profile function in Equation (20). In particular, the figures show that the wake profiles transitions from top-hats to Gaussians, and the expansion rate of downstream turbines is larger than the first turbines in each row, as seen in the LES data. As a result, at closer streamwise spacings, we expect the model to perform well compared to actuator disk model simulations; however, actuator line simulations with radially-varying induction factors and small-scale tip vortices will not be fully described by the model at close streamwise spacings.