2.1. Governing equations and boundary conditions

We consider a thermally neutral PBL, with constant pressure gradient $\boldsymbol{{\rm\nabla}}p_{\infty }/{\it\rho}\equiv f_{\infty }\boldsymbol{e}_{1}$ (with $\boldsymbol{e}_{1}$ the unit vector in the streamwise direction). The governing equations are the filtered incompressible Navier–Stokes equations for neutral flows and the continuity equation, i.e. 

where $\widetilde{\boldsymbol{u}}=[\widetilde{u}_{1},\widetilde{u}_{2},\widetilde{u}_{3}]$ is the resolved velocity field, $\widetilde{p}$ is the remaining pressure field (after subtracting $p_{\infty }$ ), ${\bf\tau}_{M}$ is the subgrid-scale model, and the density ${\it\rho}$ is assumed to remain constant. Furthermore, $\boldsymbol{f}$ represents the forces (per unit mass) introduced by the turbines on the flow (see discussion in § 2.2). Since the Reynolds number in ABLs is very high, we neglect the resolved effects of viscous stresses in the LES.

The computational domain is schematically represented in figure 1. In the streamwise and spanwise directions, periodic boundary conditions are used (i.e. respectively on ${\it\Gamma}_{1}$ and ${\it\Gamma}_{2}$ ). At the top boundary ( ${\it\Gamma}_{3}^{+}$ ), symmetry conditions are imposed. At the ground surface ( ${\it\Gamma}_{3}^{-}$ ), impermeability is imposed in combination with Schumann’s (Reference Schumann1975) wall-stress boundary conditions and Monin–Obukhov similarity theory for neutral rough boundary layers. It relates the wall stress $[{\it\tau}_{w1},{\it\tau}_{w2}]$ to the wall-parallel velocity components $[\widetilde{u}_{1},\widetilde{u}_{2}]$ at the first grid point using (Moeng Reference Moeng1984; Bou-Zeid, Meneveau & Parlange Reference Bou-Zeid, Meneveau and Parlange2005)

where $z_{0}$ is the surface roughness of the wall, and $z_{1}$ is the vertical location of the first grid point. Furthermore, the bar on $\overline{\widetilde{u}}_{1}$ and $\overline{\widetilde{u}}_{2}$ represents a local average obtained by filtering the wall-parallel velocity $[\widetilde{u}_{1},\widetilde{u}_{2}]$ in directions parallel to the wall, avoiding an overestimation of the wall stresses (Bou-Zeid et al. Reference Bou-Zeid, Meneveau and Parlange2005). In the current work, we use two successive one-dimensional Gaussian filters with filter width $4{\it\Delta}$ (and ${\it\Delta}$ the grid spacing).

Figure 1. Computational domain ${\it\Omega}$ and boundaries ${\it\Gamma}$ .

In view of the complexity associated with the formulation of the adjoint equations and adjoint subgrid-scale model required for the optimal control (cf. § 3), we opt for a relatively simple subgrid-scale model, i.e. the Smagorinsky (Reference Smagorinsky1963) model with wall damping,

with $\unicode[STIX]{x1D64E}=(\boldsymbol{{\rm\nabla}}\widetilde{\boldsymbol{u}}+(\boldsymbol{{\rm\nabla}}\widetilde{\boldsymbol{u}})^{\text{T}})/2$ the resolved rate-of-strain tensor. The Smagorinsky length $\ell$ ( $=C_{s}{\it\Delta}$ far from the wall) is damped using Mason & Thomson’s (Reference Mason and Thomson1992) wall damping function, i.e.  $\ell ^{-n}=[C_{s}{\it\Delta}]^{-n}+[{\it\kappa}(z+z_{0})]^{-n}$ , with ${\it\kappa}=0.4$ the von Kármán constant, and where we take $n=3$ . Furthermore, ${\it\Delta}=({\it\Delta}_{1}{\it\Delta}_{2}{\it\Delta}_{3})^{1/3}$ is the local grid spacing, and $C_{s}$ is the Smagorinsky coefficient. We employ $C_{s}=0.14$ , consistent with the high-Reynolds-number Lilly value for cubical sharp cutoff filters (Meyers & Sagaut Reference Meyers and Sagaut2006). Note that some other works have used more advanced subgrid-scale models in LES of wind farms, such as the scale-dependent Lagrangian model of Bou-Zeid et al. (Reference Bou-Zeid, Meneveau and Parlange2005) (Calaf et al. Reference Calaf, Meneveau and Meyers2010, Reference Calaf, Parlange and Meneveau2011; Abkar & Porté-Agel Reference Abkar and Porté-Agel2013). In Calaf et al. (Reference Calaf, Meneveau and Meyers2010) a comparison was made between this model and the current Smagorinsky implementation, and the differences in mean velocity and Reynolds stress distributions were found to be small.

2.2. Actuator-disk model

Actuator-disk models add the axial forces exerted by the wind turbines on the flow to the Navier–Stokes equations. Tangential forces are usually neglected (given the tip-speed ratios at which turbines are operated, they are more than an order of magnitude smaller). In a validation study, comparing ADM with and without tangential forces (Wu & Porté-Agel Reference Wu and Porté-Agel2011), it was demonstrated that standard ADMs provide an accurate representation of the overall wake structures behind turbines except for the very near wake ( $x/D<3$ ). Moreover, also Reynolds stresses were found to be accurately predicted, thus yielding a good representation of the interaction of the wind farms with the boundary layer. Later, this was further corroborated by Meyers & Meneveau (Reference Meyers and Meneveau2013) in a detailed analysis of energy fluxes in wind farms, comparing models with and without tangential forces. In the current work, we employ a standard ADM. It corresponds to the version used by Calaf et al. (Reference Calaf, Meneveau and Meyers2010), Meyers & Meneveau (Reference Meyers and Meneveau2010) and Meyers & Meneveau (Reference Meyers and Meneveau2013), and is briefly reviewed below.

The axial force of a turbine $i~(=1,\dots ,N_{t})$ on the flow field can be expressed as

with $C_{T,i}^{\prime }$ the disk-based thrust coefficient, $\widehat{V}_{i}$ the average axial flow velocity at the turbine rotor disk (see further below) and $A={\rm\pi}D^{2}/4$ the rotor-disk surface. The disk-based thrust coefficient $C_{T,i}^{\prime }$ results from integrated lift and drag coefficients over the turbine blades, taking design geometry and flow angles into account (cf. appendix A for a detailed formulation). For an ideal design, and in the absence of any drag forces, $C_{T,i}^{\prime }=2$ . Moreover, below rated wind speed, conventional turbines use generator torque control to keep the turbine at a constant optimal tip-speed ratio independent of wind speed, while the blade pitch is kept constant at its optimal design value. In an ADM, this corresponds to using a constant value for $C_{T,i}^{\prime }$ (see also (A 5) in appendix A).

In an ideal turbine design, the force $F_{i}$ is uniformly distributed over the disk area. Therefore, in an ADM, a uniform force (per unit mass) is distributed over the LES grid cells in the vicinity of the actuator disk using (Calaf et al. Reference Calaf, Meneveau and Meyers2010; Meyers & Meneveau Reference Meyers and Meneveau2010, Reference Meyers and Meneveau2013)

where $\boldsymbol{e}_{\bot }$ is the unit vector perpendicular to the turbine disk, and in (2.2) we employ $\boldsymbol{f}=\sum \boldsymbol{f}^{(i)}$ .

Further, $\mathscr{R}_{i}(\boldsymbol{x})$ is a geometrical smoothing function that distributes the uniform surface force of turbine $i$ over the surrounding LES grid cells. To this end, a Gaussian filter is used, leading to (Meyers & Meneveau Reference Meyers and Meneveau2010)

with $G(\boldsymbol{x})=[6/({\rm\pi}{\it\Delta}_{R}^{2})]^{3/2}\exp (-6\Vert \boldsymbol{x}\Vert ^{2}/{\it\Delta}_{R}^{2})$ the Gaussian filter kernel, with filter width ${\it\Delta}_{R}$ . Further, $\boldsymbol{x}_{i}$ is the coordinate of the turbine rotor centre, ${\it\delta}(x)$ the Dirac delta function, $\text{H}(x)$ the Heaviside function, and where $\boldsymbol{y}^{\prime }=(\boldsymbol{x}^{\prime }-\boldsymbol{x}_{i})-((\boldsymbol{x}^{\prime }-\boldsymbol{x}_{i})\boldsymbol{\cdot }\boldsymbol{e}_{\bot })\boldsymbol{e}_{\bot }$ is the projection of $(\boldsymbol{x}^{\prime }-\boldsymbol{x}_{i})$ on the rotor plane. Similar to earlier studies (Calaf et al. Reference Calaf, Meneveau and Meyers2010; Meyers & Meneveau Reference Meyers and Meneveau2010, Reference Meyers and Meneveau2013), we select ${\it\Delta}_{R}=3{\it\Delta}/2$ , with ${\it\Delta}$ the LES grid resolution. Finally, note that by construction $\int _{{\it\Omega}}\mathscr{R}_{i}(\boldsymbol{x})\,\text{d}\boldsymbol{x}^{\prime }=A$ .

In order to determine the axial disk-averaged velocity $\widehat{V}_{i}$ , first spatial averaging of the velocity is performed over the rotor disk using the geometrical rotor footprint $\mathscr{R}_{i}(\boldsymbol{x})$ , followed by a local time filter. Thus, first the disk-averaged velocity is defined as

and $\widehat{V}_{i}$ is obtained from $V_{i}$ using a first-order time filter, i.e. 

with ${\it\tau}$ the filter time scale. In our simulations, this ordinary differential equation is discretized using an implicit Euler method, such that

with ${\it\alpha}={\rm\Delta}t/({\it\tau}+{\rm\Delta}t)$ . We use ${\it\tau}=5~\text{s}$ ; combined with a time step ${\rm\Delta}t$ of $0.7~\text{s}$ (see § 2.3) this yields ${\it\alpha}\approx 0.12$ .

Finally, the power that is extracted from the boundary layer by all turbines is expressed as

This is not equivalent to the power $P_{ax}$ that is extracted at the turbine axle, which is related to the torque and rotational velocity of the turbine. The drag forces on the turbine blade increase the thrust force, but reduce the torque. Similar to $C_{T,i}^{\prime }$ , a disk-based power coefficient $C_{P,i}^{\prime }$ may be defined that is based on projected forces in the tangential direction. In the absence of drag, $C_{T,i}^{\prime }=C_{P,i}^{\prime }$ (cf. appendix A for details). In the current work, we focus on increasing $P$ by controlling $C_{T,i}^{\prime }$ (cf. § 3), and do not explicitly take $C_{P,i}^{\prime }$ into account. We just presume that $C_{T,i}^{\prime }/C_{P,i}^{\prime }$ is roughly constant, so that the extracted power from the boundary layer is representative of the mechanical power at the turbine axle. Such an approximation does not take into account deleterious effects that increased turbulence levels may have on local blade lift and drag coefficients, in particular, as a result of increased occurrences of stall. However, as further shown in §§ 4 and 5, in the considered optimal control cases, turbulence levels do not increase in front of the turbines, so that the above working assumption is reasonable. A more involved representation that more accurately models the effect of turbulence levels on blade performance is a subject for further research (see also discussion in § 6).

2.3. Discretization and case set-up

Simulations are performed in SP-Wind, an in-house research code that was developed in a series of earlier studies on LES and wind-farm simulations, and flow optimization (see e.g. Meyers & Sagaut Reference Meyers and Sagaut2007; Delport et al. Reference Delport, Baelmans and Meyers2009; Meyers & Meneveau Reference Meyers and Meneveau2010). SP-Wind uses a pseudospectral discretization in the horizontal directions. The nonlinear convective terms and the subgrid-scale stress are dealiased using the $3/2$ rule (Canuto et al. Reference Canuto, Hussaini, Quarteroni and Zang1988). Message passing interface (MPI) is used to run the simulations in parallel mode, and the FFTW library is employed for Fourier transforms (Frigo & Johnson Reference Frigo and Johnson2005). In the vertical direction, a fourth-order energy-conservative finite-difference discretization is used (Verstappen & Veldman Reference Verstappen and Veldman2003), while time integration is performed using a classical four-stage fourth-order Runge–Kutta scheme.

We focus on the simulation and optimal control of an aligned wind farm, i.e. turbines are aligned in rows that are parallel to the wind direction. Details of the case set-up are summarized in table 1 using typical orders of magnitude that are relevant for a wind farm. We take a boundary-layer height of $H=1~\text{km}$ , and use a domain size of $L_{x}\times L_{y}\times H=7~\text{km}\times 3~\text{km}\times 1~\text{km}$ . Fifty turbines with diameter $D=100~\text{m}$ are arranged in a $10\times 5$ matrix, with streamwise spacing $S_{x}=7D$ and spanwise spacing $S_{\text{y}}=6D$ . As routinely done, the set-up can also be non-dimensionalized with turbine hub height or boundary-layer height, and with friction velocity $(\,f_{\infty }H)^{1/2}$ ( $=1~\text{m}~\text{s}^{-1}$ here). The computational grid corresponds to $N_{x}\times N_{y}\times N_{z}=256\times 192\times 80$ . For dealiasing, this is extended to $384\times 288\times 80$ for all operations in real space.

The resulting case resembles earlier aligned wind-farm simulations (cf. case A3 in Calaf et al. (Reference Calaf, Meneveau and Meyers2010), and case 1 in Meyers & Meneveau (Reference Meyers and Meneveau2013)), but we slightly altered the wind-farm parameters, so that the turbine spacings are integer multiples of the rotor diameter, while keeping the ground surface per turbine roughly the same ( $S_{x}\times S_{y}=7.85D\times 5.23D$ in the earlier studies, with $8\times 6$ turbines). Furthermore, we also used a slightly finer mesh spacing. We refer the reader to Calaf et al. (Reference Calaf, Meneveau and Meyers2010) and Meyers & Meneveau (Reference Meyers and Meneveau2013) for the effects of domain size and grid refinement studies. Finally, for the time integration, we use a fixed time step of 0.7 s, which corresponds to a Courant–Friedrichs–Lewy (CFL) number of approximately 0.4.

For conventional simulations, i.e. those discussed in § 2.4, we first start from an initial logarithmic velocity profile to which a set of random perturbations are added. After an initial spin-up period of 16 h (corresponding to approximately 85 through-flow times), during which the velocity profile and turbulence statistics evolve into statistical equilibrium, we accumulate averaged flow properties for a time window of 21 h. Subsequent simulations (for parameter variations, cf. § 2.4) start from an earlier statistically stationary field, and use a spin-up period of 6 h (30 through-flow times).

In figure 2, a snapshot of the instantaneous velocity fields is shown. In the horizontal plane, we observe typical meandering of the turbine wakes. At the same time, patches of high-speed wind can also be seen passing through the spaces between turbine columns.

Figure 2. (a) Snapshot representing an instantaneous streamwise velocity field. (b) Zoom on a subset of four turbines. The horizontal plane is taken at hub height, and turbines are represented by small white disks.

2.4. Boundary-layer response and optimal energy extraction using standard turbine control

In the current subsection, we study the response of a wind-farm boundary layer to static changes in $C_{T}^{\prime }$ , while keeping the $C_{T}^{\prime }$ value the same for all turbines. This will help in defining a reference case, as well as a starting point for the dynamic optimization in § 3.

In order to discuss this further, we first introduce some concepts related to fully developed wind-farm boundary layers with wind turbines situated in the inner layer. In such boundary layers, a double log layer is observed (Frandsen Reference Frandsen1992; Calaf et al. Reference Calaf, Meneveau and Meyers2010), one below the turbine level, characterized by the friction velocity $u_{{\it\tau}l}=({\it\tau}_{w}/{\it\rho})^{1/2}$ , and the other above the turbine level, characterized by the total friction induced by the ground surface and the turbines. Thus with friction velocity (Calaf et al. Reference Calaf, Meneveau and Meyers2010)

with $A={\rm\pi}D^{2}/4$ the turbine rotor area. In a PBL, integration of the momentum balance over the full height of the boundary layer further yields $u_{{\it\tau}h}^{2}=f_{\infty }H$ .

When considering a single turbine in idealized conditions, the optimal operating condition of the turbine corresponds to $C_{T}^{\prime }=2$ (cf. appendix A), corresponding to the Betz limit. However, in an infinite wind-farm boundary layer, the boundary layer responds to the surface roughness induced by the wind farm, and the wind velocity at turbine hub height depends on parameters such as turbine spacing and thrust coefficient. Therefore, when comparing different control cases, it is important to normalize the total power extraction of the farm $P$ by a correct reference value that itself is not dependent on the control and remains constant in real conditions. The logical reference to use for a wind farm in the ABL is the geostrophic wind $G$ in the free atmosphere above the ABL. This approach was followed by Meyers & Meneveau (Reference Meyers and Meneveau2012), for example, when investigating optimal turbine spacing for large wind farms, thus optimizing $P_{ABL}^{+}=P/G^{3}$ .

An issue is that the boundary layer considered in the current study is a regular PBL. However, we can use the main working hypothesis that the wind turbines are in the inner layer of the boundary layer (cf. discussion in the introduction), and thus their overall effect on the outer layer is characterized by the friction velocity $u_{{\it\tau}h}$ . By further using basic momentum and energy conservation laws for an Ekman spiral, it is then possible to obtain a simple heuristic relation between $u_{{\it\tau}h}$ and $G$ , i.e. (cf. appendix B for a derivation)

with $\mathscr{D}$ the total turbulent dissipation per unit wind-farm area, and $\mathscr{P}$ the average turbine power extraction per unit farm area (i.e.  $\mathscr{P}=\overline{P}/(L_{x}L_{y})$ , with $\overline{P}$ the time average of $P$ ). Furthermore, $A\approx 12$ is an empirical constant that depends on the outer-layer behaviour of the ABL.

The response of the boundary layer to changes in $C_{T}^{\prime }$ is now investigated. We do not yet consider optimal coordinated control of $C_{T}^{\prime }$ at farm level (cf. § 3 below), but instead keep $C_{T}^{\prime }$ constant in time, and the same for all turbines. Thirteen different cases are considered, with $0.02\leqslant C_{T}^{\prime }\leqslant 3.5$ . Results for the averaged total power extraction are shown in figure 3. We show two normalizations, i.e. one using $\mathscr{P}_{PBL}^{+}=\mathscr{P}/u_{{\it\tau}h}^{3}$ , which is the standard normalization for a PBL, and the other using $\mathscr{P}_{ABL}^{+}=\mathscr{P}/G^{3}$ , as relevant for ABLs. For the second normalization, (2.14) is used to determine $u_{{\it\tau}h}/G$ . These two normalizations reflect the different reactions of a PBL and an ABL to a changing load. In figure 3, it is appreciated that the extracted power depends strongly on the disk-based thrust coefficient. Furthermore, the selected normalization leads to quite different optimal values for $C_{T}^{\prime }$ : using the maximum of $\mathscr{P}/u_{{\it\tau}h}^{3}$ leads to much lower optimal $C_{T}^{\prime }$ values than using $\mathscr{P}_{ABL}^{+}=\mathscr{P}/G^{3}$ . Thus, presuming a constant pressure gradient, independent of the wind-farm load, does not lead to the same optimum as presuming constant geostrophic wind.

Figure 3. Mean power output of an uncontrolled wind farm as a function of $C_{T}^{\prime }$ . Power $\mathscr{P}^{+}$ is normalized by either $u_{{\it\tau}h}$ (red circles), geostrophic wind $G$ (green squares), or driving power (blue triangles). Curves are further normalized by their maximum values of  $\mathscr{P}^{+}$ .

In figure 3, we also consider a third normalization that is based on the total driving power of the PBL, i.e.  $\mathscr{P}_{DRP}^{+}=\mathscr{P}/(\,f_{\infty }U_{b}H)$ , with $U_{b}$ the total bulk velocity. Thus, such a normalization presumes a constant driving power in the boundary layer, which is independent of the wind-farm load. It is observed that this third normalization leads to an optimal $C_{T}^{\prime }$ value that is close to the one found for $\mathscr{P}/G^{3}$ , with $C_{T}^{\prime }\approx 1.33$ .

3.1. Receding-horizon approach

We employ a receding-horizon optimal control approach for the control of wind-farm boundary layer interaction. This essentially follows the standard paradigm of model-predictive control (Rawlings & Mayne Reference Rawlings and Mayne2008), but where the model in our case consists of the full LES equations (2.1) and (2.2), and where we do not have a state estimation problem, i.e. in our simulations the flow state is perfectly known at each time step. In the context of DNS-based and LES-based optimal control, a similar setting was employed by Chang & Collis (Reference Chang and Collis1999) and Bewley et al. (Reference Bewley, Moin and Temam2001).

In a receding-horizon optimal control approach, time is split into a number of control windows with length $T$ , also called the time horizon – a schematic overview is presented in figure 5. Starting with the first time horizon, an optimization problem is formulated (cf. § 3.2) in which the control parameters are optimized as a function of time. To that end, an iterative gradient-based optimization approach is used (cf. § 3.3) requiring several LES, combined with adjoint simulations for the determination of the gradients (cf. § 3.4). Once a set of optimal controls is found for the interval $[0,T]$ , they are effectively used as control inputs to advance the system over a time window $T_{A}$ (see figure 5). Subsequently, a new optimization problem is formulated that optimizes the controls for the time window $[T_{A},T_{A}+T]$ , and so forth.

Figure 5. Receding horizon optimal control approach.

In standard receding-horizon optimal control, $T_{A}$ often just corresponds to the control time step, so that with every time step a control optimization problem is solved, leading to an optimization time horizon that smoothly moves forward with the control inputs. In the context of DNS-based or LES-based optimal control, this is not done, as every optimization problem itself requires very large computational resources. Bewley et al. (Reference Bewley, Moin and Temam2001) used $T_{A}=T$ to limit computational cost, though this led to non-smooth transitions between time windows. They also looked at one case with $T_{A}=T/2$ ; and in a similar study, Collis et al. (Reference Collis, Chang, Kellogg and Prabhu2000) also explored $T_{A}=3T/4$ and $T_{A}=T/4$ . In the current work, we choose $T_{A}=T/2$ as an ad hoc balance between computational cost and control smoothness, and we will refer to time windows $T_{A}$ as the control windows.

3.2. Optimization problem formulation

We consider the optimal control of a wind-farm boundary layer. The control parameters correspond to all disk-based turbine thrust coefficients ${\bf\varphi}\equiv [C_{T,1}^{\prime }(t),C_{T,2}^{\prime }(t),\dots ,C_{T,N_{t}}^{\prime }(t)]$ (with $N_{t}$ the total number of turbines). The state variables in the optimal control problem are $\boldsymbol{q}\equiv [\widetilde{\boldsymbol{u}}(\boldsymbol{x},t),\widetilde{p}(\boldsymbol{x},t),\widehat{\boldsymbol{V}}(t)]$ , i.e. corresponding to the LES velocity field, pressure field and the time-filtered turbine-disk velocity fields $\widehat{\boldsymbol{V}}\equiv [\widehat{V}_{1},\dots ,\widehat{V}_{N_{t}}]$ .

The optimal control problem is formulated as a minimization problem in which we employ the following cost functional:

The first term corresponds to the amount of energy extracted from the boundary layer by the wind turbines over the optimization time horizon $T$ . The second term is a penalization term ( $0\leqslant {\it\gamma}<1$ ) that is related to penalization of the total dissipation $\mathscr{D}L_{x}L_{y}$ , but formulated in an alternative way that simplifies implementation. Details are discussed in § 5, where results for ${\it\gamma}>0$ are presented. In § 4, ${\it\gamma}=0$ is discussed.

The optimal control problem under consideration is a PDE-constrained optimization problem that corresponds to

In (3.3), the wall-stress model is explicitly added to the momentum equation using the Dirac delta function ${\it\delta}(\boldsymbol{x}-z_{1}\boldsymbol{e}_{3})$ with $z_{1}$ the location of the first grid point near the wall, and where ${\bf\tau}_{\boldsymbol{w}}=[{\it\tau}_{w1},{\it\tau}_{w2},0]$ , with ${\it\tau}_{w1}$ and ${\it\tau}_{w2}$ defined by (2.3) and (2.4). For the sake of further use, the state constraints (3.3) are written in short-hand notation as $\boldsymbol{B}({\bf\varphi},\boldsymbol{q})=0$ , representing LES momentum and continuity equations, and the time filter of the disk velocity.

Finally, we remark that we add some additional box constraints on the controls, i.e.  $0\leqslant C_{T,i}^{\prime }(t)\leqslant 4$ . These are trivial to add, and are not formally included here so as not to further complicate the equations. See § 3.5 for further discussion on these constraints.

3.3. Optimization method

In the current work, we do not solve the PDE-constrained optimization problem written in its standard form (3.2) and (3.3), where the PDE is explicitly formulated as a constraint. Though it is possible and sometimes beneficial to do this for smaller problems (see e.g. Hinze & Kunisch (Reference Hinze and Kunisch2001) for a discussion), the size of the space–time state space in our optimal control problem (of the order of 1 billion degrees of freedom) does not allow such an approach. Instead, we reformulate the problem in a reduced form, with a reduced cost functional, i.e. 

where $\boldsymbol{q}({\bf\varphi})$ is the solution to the state equations given the control inputs ${\bf\varphi}$ , implicitly defined by $\boldsymbol{B}({\bf\varphi},\boldsymbol{q}({\bf\varphi}))\equiv 0$ . Thus, in its reduced form the problem is unconstrained, but at every step of the optimization algorithm the state constraints need to be explicitly satisfied. The size of the optimization space in this reduced formulation corresponds to the number of degrees of freedom on ${\bf\varphi}$ , which is approximately $2\times 10^{4}$ in this study.

To solve (3.4) the same approach is followed as first used by Bewley et al. (Reference Bewley, Moin and Temam2001) for DNS-based optimal control, i.e. the combination of a Polak–Ribière conjugate-gradient method and the Brent line-search algorithm (Press et al. Reference Press, Teukolsky, Vetterling and Flannery1996; Luenberger Reference Luenberger2005; Nocedal & Wright Reference Nocedal and Wright2006). It is an iterative method for solving unconstrained optimization problems. Given an intermediate estimate of the optimum ${\bf\varphi}^{(k)}$ , a search direction ${\it\delta}{\bf\varphi}^{(k)}$ is determined using the Polak–Ribière conjugate-gradient direction

where $\boldsymbol{{\rm\nabla}}\!\widetilde{\mathscr{J}}^{(k)}$ is the gradient of the cost functional (cf. § 3.4 for its determination based on the adjoint equations), and ${\it\beta}_{k}$ is given by

Using the search direction ${\it\delta}{\bf\varphi}^{(k)}$ , a new estimate of the optimum is obtained from

where ${\it\alpha}$ is the result of a line search that minimizes $\widetilde{\mathscr{J}}({\bf\varphi}^{(k)})$ in the direction ${\it\delta}{\bf\varphi}^{(k)}$ . To that end, an iterative gradient-free line-search method is used that is based on the mnbrak and Brent algorithms (Press et al. Reference Press, Teukolsky, Vetterling and Flannery1996). Details on the implementation used in our work can be found in Delport et al. (Reference Delport, Baelmans and Meyers2009).

3.4. Gradient of the reduced cost functional and adjoint equations

An important element in the conjugate-gradient algorithm discussed above is the determination of the gradient of the reduced cost functional $\boldsymbol{{\rm\nabla}}\!\widetilde{\mathscr{J}}$ for a given set of controls ${\bf\varphi}$ . The use of a simple finite difference approach is not feasible if the design space ${\bf\varphi}$ is large, since this requires an evaluation of the state equations for every possible dimension in ${\bf\varphi}$ . Instead, a mathematically equivalent formulation for the determination of the gradient can be used that requires once the solution of an additional set of partial differential equations, i.e. the adjoint equations, with a cost that is roughly equivalent to that of the original state equations.

The derivation of the gradient of the reduced cost functional in terms of the adjoint solution is straightforward, but lengthy. It follows standard approaches as discussed in Tröltzsch (Reference Tröltzsch2010) and Borzi & Schultz (Reference Borzi and Schultz2012), for example. We refer the reader to appendix C for a detailed derivation. Here, we summarize the most important results that are used in our optimization algorithm. First of all, the gradient of the reduced cost functional may be expressed as (Borzi & Schultz Reference Borzi and Schultz2012)

with $[\partial \boldsymbol{B}/\partial {\bf\varphi}]^{\ast }$ the adjoint of $\partial \boldsymbol{B}/\partial {\bf\varphi}$ (cf. § C.1) and $\boldsymbol{q}^{\ast }=({\bf\xi},{\rm\pi},{\bf\chi})$ the solution of the adjoint equations (see further below). For the current cost functional (3.1), this leads to

with $\textstyle \pmb{\mathscr{R}}\equiv [\mathscr{R}_{1},\dots ,\mathscr{R}_{N_{t}}]$ , and where $\circ$ is used to denote the entry-wise product (or Hadamard product), and $\widehat{\boldsymbol{V}}^{\circ 2}$ is the entry-wise square of $\widehat{\boldsymbol{V}}$ . Furthermore, ${\bf\xi}(\boldsymbol{x},t)$ is the adjoint velocity field that is obtained by solving the adjoint equations.

The derivation of the adjoint equations for the standard transient Navier–Stokes equations is well documented in the literature (see e.g. Choi, Hinze & Kunisch (Reference Choi, Hinze and Kunisch1999), Bewley et al. (Reference Bewley, Moin and Temam2001), Wei & Freund (Reference Wei and Freund2006) and Delport et al. (Reference Delport, Baelmans and Meyers2009), among others). The main extension here is the addition of adjoints for the Smagorinsky model and wall-stress model, and the adjoint of the ADM. The reader is referred to appendix C for details of their derivation. The resulting adjoint equations correspond to

and with

and ${\bf\tau}_{\boldsymbol{w}}^{\ast }=[{\it\tau}_{w,1}^{\ast },{\it\tau}_{w,2}^{\ast },0]$ . Further

where $\unicode[STIX]{x1D64E}^{\ast }=(\boldsymbol{{\rm\nabla}}{\bf\xi}+(\boldsymbol{{\rm\nabla}}{\bf\xi})^{\text{T}})/2$ .

The spatial boundary conditions of the adjoint equations are equivalent to those of the forward equations. In the $x_{1}$ and $x_{2}$ directions, periodic boundary conditions are required. In the normal direction, impermeability is required at the top and bottom walls, and symmetry boundary conditions at the top wall. For the ‘initial conditions’, it is important to realize that the adjoint equations are solved backwards in time (cf. the sign of the time derivatives), so that the ‘initial’ conditions should be provided at $t=T$ . They correspond to

The adjoint equations (3.10) show some similarity to the flow equations of the forward problem, e.g. time derivatives and convective terms can be recognized (though with different signs), continuity looks the same, and there is also an adjoint pressure variable. Therefore, much of the discretization of the forward problem can be reused, with the same pseudospectral discretization in the horizontal directions, in combination with a fourth-order energy-conservative discretization in the vertical direction. For the time integration, a fourth-order Runge–Kutta method is also used. Note that the adjoint equations follow from a linearization of the governing equations around a state $(\widetilde{\boldsymbol{u}},\widetilde{p},\widehat{\boldsymbol{V}})$ . In the adjoint equations, this state is also required (cf. (3.10)). To this end, the nonlinear forward problem is solved first, and the full space–time state is stored on disk. Subsequently, it is used during the solution of the adjoint equations.

Finally, for the discretization of the filter equation, we choose a time discretization that is equivalent to the discrete adjoint of the discrete forward filter equation (2.11). This corresponds to

3.5. Computational set-up

The geometrical set-up, grid, time step, etc., remain the same as in § 2.3 (cf. table 1). The optimal control is started from a statistically stationary field of an ‘uncontrolled’ wind farm with $C_{T,i}^{\prime }=1.33$ , and ${\bf\varphi}^{(0)}=1.33$ is used for the starting value of the optimization algorithm.

As already introduced in § 3.2, we use box constraints on the controls, i.e.  $0\leqslant C_{T,i}^{\prime }(t)\leqslant 4$ . The lower constraint prevents the turbine from starting to operate as a fan, even if the optimization algorithm would ask for this. For the upper boundary, we do not a priori want to limit $C_{T}^{\prime }$ to the Betz limit ( $C_{T}^{\prime }=2$ ), but at the same time we cannot leave it free, as $C_{T}^{\prime }\rightarrow \infty$ is not very practicable from a turbine construction point of view. Therefore, we selected an ad hoc limit of $C_{T}^{\prime }=4$ , which corresponds to a wind turbine that is constructed with double blade chord lengths compared to Betz-optimal blade design (cf. § A.2 for details). Moreover, we have further investigated an extra case with optimal control over one control window without box constraints on $C_{T}^{\prime }$ (and using ${\it\gamma}=0$ ). For this case, we found that $C_{T}^{\prime }$ fluctuates between $-19$ and $+24$ , but compared to the case with box constraints this leads to no significant additional energy extraction (i.e. 17.7 % extra energy for the case without versus 17.66 % for the case with box constraints, averaged over the first control window).

For the optimal receding-horizon control cases considered here, we select an optimization time horizon $T=280~\text{s}$ . This corresponds approximately to 0.4 times the through-flow time, or the average convection time for the flow to pass four rows of turbines. In this way, we avoid any interaction between the optimal control approach and the periodicity of the streamwise boundary conditions. Once the optimization of the controls is converged, they are used for $T_{A}=140~\text{s}$ , before the next optimization problem is started. This process is repeated for a total of 25 control windows, totalling 3500 s (approximately 1 h) of wind-farm operation.

Finally, in order to limit computational costs, we do not fully converge the conjugate-gradient algorithm or the line-search algorithms. Instead, we stop after five conjugate-gradient iterations, and use a maximum of eight forward function evaluations per line search. This leads to a maximum of 45 PDE simulations per control window (40 forward and five adjoint), or 1125 PDE simulations in total, where one PDE simulation takes approximately 90 min of wall time on 32 processors.

In figure 6, we show a typical convergence history of an optimization (using ${\it\gamma}=0$ ) for three different control windows. On the $x$ axis the number of subsequent PDE simulations during the iterative conjugate-gradient optimization algorithm is shown. Closed symbols refer to standard LES, while open symbols refer to adjoint simulations. It is appreciated that the cost function does not decrease monotonically with the number of simulations. This is related to the line-search algorithm (cf. § 3.3), which sometimes overshoots the optimal step length along a conjugate-gradient search direction. In figure 6 it is appreciated that cost functionals decrease significantly during the optimization, but optimization is stopped after five conjugate-gradient iterations, so they are not formally converged to $\boldsymbol{{\rm\nabla}}\!\widetilde{\mathscr{J}}=0$ .

Figure 6. Typical convergence history of the conjugate-gradient optimization for three different control windows $[(n-1)T_{A},(n-1)T_{A}+T]$ : ▪, control window $n=1$ ; ▴, control window $n=13$ ; ●, control window $n=20$ . Open symbols (○, ▫, ▵) correspond to adjoint simulations required for gradient evaluations, and are plotted at the same cost-functional level as the previous forward simulation.

4.1. Adjoint fields and optimal controls

In figure 7, snapshots of the instantaneous adjoint fields are shown. Unlike the flow field, the adjoint equations evolve backwards in time, and propagate in the upstream direction. The adjoint field depends strongly on the definition of the cost functional, and the forward flow state around which the equations are linearized. The fields that are shown in figure 7 belong to the first adjoint equations in the optimization sequence of the first optimal control time horizon (see figure 5): at this point, the equations are linearized around a flow state that is obtained for initially constant controls at all turbines with $C_{T,i}^{\prime }(t)=1.33~(i=1,\dots ,N_{t})$ . The initial condition for the adjoint equations (with ${\it\gamma}=0$ ) corresponds to ${\bf\xi}(\boldsymbol{x},T)=0$ . This is also visible in the first snapshot at $T-t=14$ (figure 7 a), where the field is still largely zero.

Figure 7. Contours of instantaneous streamwise adjoint fields: (a)  $T-t=14~\text{s}$ ; (b)  $T-t=70~\text{s}$ ; (c)  $T-t=174~\text{s}$ ; (d) zoom of (b). Horizontal planes are taken at hub height.

The adjoint equations are driven by the cost function of the optimization problem (i.e.  $\partial \mathscr{J}/\partial \boldsymbol{q}$ in (C 13)). They essentially express where possible changes in the cost function $\mathscr{J}$ may be originating from (thus the equations evolve backwards in time and in the reverse flow direction). When looking at the first two snapshots in figure 7 (at $T-t=14$ and $70$ ), we see that changes to the cost functional in the last time interval of the optimal control time horizon originate from a tube upstream of the turbine rotors. At this point, only changes to the flow velocity in this region affect the later energy extraction at the turbine rotor. When looking earlier in time ( $T-t=174$ ), the ‘tube’ observed in figure 7(b) has ‘hit’ the previous row of turbines, so that upstream turbines have the potential to influence the energy extraction of their downstream neighbours. Looking at figure 7(c) (at $T-t=174$ ), it is observed that the adjoint field has become fully turbulent in the whole domain. This shows that it is the full interaction with the boundary layer that influences the wind-farm energy extraction.

In figure 8 the behaviour of the optimal thrust coefficient is shown for one of the turbines in the controlled wind farm. Approximately 1 h of time is shown, corresponding to subsequent optimal control in 25 control windows. It is appreciated that $C_{T}^{\prime }$ is strongly changing in time, but remains limited by the box constraints used during the optimization (cf. § 3.5). A further zoom in the figure reveals that the changes in $C_{T}^{\prime }$ are well resolved in time; no additional smoothing of the gradients used in the optimization was required for this. Moreover, typical time scales with which the controls change remain above 10 s.

Figure 8. Time evolution of the thrust coefficient of one of the turbines in the farm.

4.2. Optimized power output

In the current subsection, the energy balance and power extraction are discussed in detail. In figure 9 the total wind-farm power extraction per unit wind-farm area $\mathscr{P}_{{\it\Omega}}$ ( $=P/(L_{x}L_{y})$ ) and the total gains and losses per unit area are shown. To this end, the total kinetic energy equation is horizontally averaged and integrated over the boundary layer height, leading to

where we use $\langle \cdots \rangle _{{\it\Gamma}}$ to denote horizontal averages, and the subscript ${\it\Gamma}$ to indicate that averages are taken over a finite-domain horizontal plane ${\it\Gamma}\subset {\it\Omega}$ . In contrast to $\mathscr{P}$ and $\mathscr{D}$ in § 2.4, $\mathscr{P}_{{\it\Omega}}$ and $\mathscr{D}_{{\it\Omega}}$ are not averaged in time, so that they still fluctuate significantly from one time step to another (i.e. the horizontal extent of our domain is not large enough to obtain statistical convergence based on horizontal averaging only).

Figure 9. Time evolution of (a) total wind-farm power output and (b) gains and losses: – – –, driving power by pressure force; —— (grey), rate of change of kinetic energy; —— (black), farm power; - $\cdot$ - $\cdot$ -, dissipation.

In figure 9(a) the power extraction is shown before control ( $t<0$ ) and after the start of the optimal model predictive control. It is seen that overall the power extraction increases, but starts to fluctuate significantly more than before coordinated optimal control. On average, a gain in energy extraction of 16 % is achieved (integrated between $t=0$ and $t=3500~\text{s}$ ). In figure 9(b) gains and losses to the boundary layer are shown. Here, the wind-farm power is plotted as a loss term, together with the dissipation $\mathscr{D}_{{\it\Omega}}$ . We observe that the increase in power extraction is mainly balanced by an overall deceleration of the flow. In addition, also the dissipation $\mathscr{D}_{{\it\Omega}}$ increases (recall that dissipation is not penalized, i.e.  ${\it\gamma}=0$ ).

In figure 9(b) the driving power $f_{\infty }U_{b}H$ is also plotted. It is observed that the driving power slowly decreases during the optimal control. This is directly related to the fact that $f_{\infty }$ is kept constant during the optimization, while the flow slowly decelerates. Returning to the discussion in § 2.4, we remark here that $f_{\infty }U_{b}H$ is not explicitly kept constant during optimization. This would require adding a non-trivial state constraint to (3.2), requiring the solution of an additional set of adjoint equations for every gradient evaluation. However, it is appreciated that changes in the driving power remain small. Moreover, we also remark that the average level of $C_{T}^{\prime }$ in the optimal control moves to higher values, and thus away from the constant $C_{T}^{\prime }$ optimum observed in figure 3 for the constant-forcing case. This clearly indicates that the optimal controls found here do not lead to a statistically stationary optimal situation, but instead purely exploit transient boundary-layer deceleration.

4.3. Energy balance in the turbine region

The energy balances in the wind farm are further investigated. To smooth the results, they are additionally integrated per control window $T_{A}$ . Next to the balance integrated over the whole height of the computational domain ${\it\Omega}$ (cf. (4.1)), we also look at the balance in the turbine region, i.e. integrated from $z_{h}-D/2$ to $z_{h}+D/2$ , and we denote this horizontal slab of the computational domain with ${\it\Omega}_{D}$ (note that, for numerical evaluation, we take the integration bounds slightly wider, to ensure that all filtered turbine forces (cf. (2.7)) are included in ${\it\Omega}_{D}$ ).

The following notation is used for the horizontal and time average (here for $\widetilde{\boldsymbol{u}}$ )

using $\overline{\,\cdot \,}^{T_{A}}$ to denote the time average over a time window with length $T_{A}$ . Moreover, define $\widetilde{\boldsymbol{u}}^{\prime \prime }=\widetilde{\boldsymbol{u}}-\overline{\langle \widetilde{\boldsymbol{u}}\rangle }_{{\it\Gamma}}^{T_{A}}$ , $e=\widetilde{\boldsymbol{u}}\boldsymbol{\cdot }\widetilde{\boldsymbol{u}}/2$ and $e^{\prime \prime }=e-\overline{\langle e\rangle }_{{\it\Gamma}}^{T_{A}}=\widetilde{\boldsymbol{u}}^{\prime \prime }\boldsymbol{\cdot }\widetilde{\boldsymbol{u}}^{\prime \prime }/2+\widetilde{\boldsymbol{u}}^{\prime \prime }\boldsymbol{\cdot }\overline{\langle \widetilde{\boldsymbol{u}}\rangle }_{{\it\Gamma}}^{T_{A}}$ . The energy balance for ${\it\Omega}_{D}$ is then expressed as

Here we have also further introduced the following notation: turbulent (and dispersive) transport $\overline{T}^{\prime \prime }$ , driving power $\mathscr{F}_{{\it\Omega}_{D},T_{A}}$ , wind-farm power extraction $\mathscr{P}_{{\it\Omega},T_{A}}$ and dissipation $\mathscr{D}_{{\it\Omega}_{D},T_{A}}$ , all per unit farm area. The terms of (4.3) are plotted in figure 10(b). For reference, in figure 10(a) the same terms are shown integrated for the whole domain: this corresponds to the balance shown in figure 9(b), but in addition averaged in time per control window.

Figure 10. Gains and losses per unit wind-farm area averaged over control windows for (a) the whole computation domain ${\it\Omega}$ and (b) the disk region ${\it\Omega}_{D}$ : ♦, rate of change of kinetic energy ${\rm\Delta}E_{{\it\Omega}}/T_{A}$ ; ▪, driving power $\mathscr{F}_{{\it\Omega},T_{A}}$ ; ▴, wind-farm power extraction $\mathscr{P}_{{\it\Omega},T_{A}}$ ; ●, dissipation $\mathscr{D}_{{\it\Omega},T_{A}}$ ; $+$ , total transport $\overline{T}^{\prime \prime }(z+D/2)-\overline{T}^{\prime \prime }(z-D/2)$ .

Looking at figure 10(b), it is observed that the turbine region is in equilibrium for $t\geqslant 5T_{A}$ , i.e. in this region ${\rm\Delta}E_{{\it\Omega}_{D}}(t)$ is approximately zero. The farm power extraction is the largest sink in this region, followed by the dissipation, which amounts to 37 % of the dissipation in the whole boundary layer. The main source of energy is the turbulent transport $\overline{T}^{\prime \prime }(z+D/2)$ at the top boundary of the region.

The results in figure 10 suggest that the inner region of the boundary layer has found a new statistical equilibrium after roughly five control windows (which is roughly equivalent to one through-flow time), and only the outer layer is decelerating. Moreover, the energy is transported from the outer region towards the turbine region by increased turbulent and dispersive transport. This is further investigated next.

4.4. Flow profiles in the controlled wind farm

In the current subsection, we investigate time-averaged and horizontally averaged profiles. First of all, we look at the streamwise velocity profile averaged over five different time windows, i.e. corresponding to the windows $[5(n-1)T_{A},5nT_{A}]$ with $n=1,\dots ,5$ . Results are shown in figure 11, together with the velocity profile for the uncontrolled case. For all cases, we observe two distinct logarithmic regions, one above and one below the turbines. This is consistent with observations by Cal et al. (Reference Cal, Lebron, Castillo, Kang and Meneveau2010) and Calaf et al. (Reference Calaf, Meneveau and Meyers2010) for uncontrolled wind farms. The control in the current work does not change these features. It is further observed in figure 11 that the velocity profile for the controlled cases are lower than for the uncontrolled case. When looking at the inner layer ( $z/H<0.15$ ), it is seen that the velocity profiles of the middle three averages cluster around a new equilibrium position (see inset in figure). The first average $[0,5T_{A}]$ is found to be somewhat higher (closer to the uncontrolled case), while the last average $[20T_{A},25T_{A}]$ , is somewhat lower than the three previous averages. For the current case, we pushed the optimal control a bit further, i.e. up to 33 control windows (results not shown in the plots), and this further confirms that the inner layer starts to decelerate again after about 20 to $25T_{A}$ . When looking higher up in the wind-farm boundary layer and close to the top, it is appreciated that the velocity profiles keep decreasing in all control windows. This is consistent with the observations in figure 10. Therefore, for analysis of higher-order moments, we will mainly show averages over the middle time window $[5T_{A},20T_{A}]$ , which extends over 35 min of wind-farm operation, during which we presume that the flow is in statistical equilibrium in the inner layer.

Figure 11. Streamwise mean velocities: ——, uncontrolled case; - $\cdot$ - $\cdot$ -, optimal control case averaged over the time interval $[0,5T_{A}]$ ; – – – (see also inset), averaged over later intervals.

In figure 12 the total, dispersive and Reynolds stresses are shown for the different cases. Recall that these stresses are constructed based on horizontal averages over a domain that is horizontally periodic, but not homogeneous. We first look at the total shear stress in figure 12(a), averaged over different windows in the controlled case, and for the uncontrolled case. For the latter it is seen that $-\langle \overline{\widetilde{u}_{1}^{\prime \prime }\widetilde{u}_{3}^{\prime \prime }}\rangle =(1-z/H)f_{\infty }H$ above the turbines, as can be expected from the plane-averaged momentum balance in a conventional channel-flow boundary layer. Below the turbine level, the turbine forces and subgrid-scale stresses (close to the wall) take over the role of the turbulent shear stress in the total balance.

Figure 12. Vertical profiles of (a) total stresses and (b–e) total stresses, Reynolds and dispersive components. In panel (a): —— (black), uncontrolled case; - $\cdot$ - $\cdot$ - (blue), – – – (black) and —— (green), controlled case, respectively averaged over time windows $[0,5T_{A}]$ , $[5T_{A},20T_{A}]$ and $[20T_{A},25T_{A}]$ . In panels (b–e): —— (black), —— (orange) and —— (cyan), respectively the total stresses, Reynolds stresses and dispersive stresses for the uncontrolled case; – – – (black), – – – (orange) and – – – (cyan), respectively the total stresses, Reynolds stresses and dispersive stresses for the controlled case averaged over $[5T_{A},20T_{A}]$ .

Figure 13. Vertical profiles of horizontally averaged mean-flow kinetic energy flux. Uncontrolled case: —— (black), total kinetic energy flux; —— (orange), flux due to Reynolds stress; —— (cyan), flux due to dispersive stress. Optimal control case averaged over time window $[5T_{A},20T_{A}]$ : – – – (black), total kinetic energy flux; – – – (orange), flux due to Reynolds stress; – – – (cyan), flux due to dispersive stress.

Looking at the total shear stresses $\langle \overline{\widetilde{u}_{1}^{\prime \prime }\widetilde{u}_{3}^{\prime \prime }}\rangle$ for the controlled cases in figure 12(a), the picture changes. Now the boundary layer is no longer in equilibrium, so that $u_{{\it\tau}h}\neq f_{\infty }H$ , but rather

with $-\text{d}U_{b}/\text{d}t$ the deceleration of the bulk flow. Since $\text{d}U_{b}/\text{d}t$ is roughly constant for $t\in [5T_{A},20T_{A}]$ (cf. figure 10), we can expect $u_{{\it\tau}h}$ also to be roughly constant in this region, but higher than $f_{\infty }H$ . This is observed in figure 12(a) for the average over the interval $[5T_{A},20T_{A}]$ . Note that for this case the deceleration mainly takes place in the outer layer of the boundary layer, while the inner layer ( $z/H<0.15$ ) is in a new equilibrium (cf. discussion above). For the average over $[0,5T_{A}]$ , this is not yet the case, and this is also apparent from figure 12(a).

In figure 12(b–e) we further decompose the total stresses (both the shear stress and the normal stresses) into dispersive stresses ( $\langle \overline{\widetilde{u}}_{i}^{\prime \prime }\overline{\widetilde{u}}_{j}^{\prime \prime }\rangle$ ) and plane-averaged Reynolds stresses ( $\langle \overline{\widetilde{u}_{i}^{\prime }\widetilde{u}_{j}^{\prime }}\rangle$ ), averaged over the time window $[5T_{A},20T_{A}]$ , i.e. (e.g. following Calaf et al. Reference Calaf, Meneveau and Meyers2010) $\langle \overline{\widetilde{u}_{i}^{\prime \prime }\widetilde{u}_{j}^{\prime \prime }}\rangle =\langle \overline{\widetilde{u}}_{i}^{\prime \prime }\overline{\widetilde{u}}_{j}^{\prime \prime }\rangle +\langle \overline{\widetilde{u}_{i}^{\prime }\widetilde{u}_{j}^{\prime }}\rangle$ , with $\widetilde{u}_{i}^{\prime }=\widetilde{u}_{i}-\overline{\widetilde{u}}_{i}$ . First of all, in figure 12(b), it is very interesting to notice that the increase in total stress is caused by an increase of dispersive shear stresses, which are roughly doubled, while the plane-averaged Reynolds stresses decrease. The same trends are observed in figure 12(c) for the streamwise stresses: dispersive stresses increase significantly, but Reynolds stresses decrease. For the spanwise and vertical stresses, we observe that dispersive stresses also increase significantly, while Reynolds stresses remain largely unchanged, except in the turbine-tip region, where they slightly increase.

Looking at fluxes of horizontally averaged mean-flow energy in figure 13, the same trends are observed. It is seen that the energy flux at the top of the farm ( $z=z_{h}+D/2$ ) is considerably higher for the controlled case than for the uncontrolled case. Below the farm, the inverse is observed. Here the energy flux towards the flow below the farm is decreased compared to the uncontrolled case. Thus, as a result of the optimal control, the energy flux towards the farm increases, while at the same time this energy is better captured by the turbines. We further also looked at total kinetic energy fluxes (not shown here), which contain additional elements of triple dispersive correlations, and turbulent and pressure diffusion, but these are minor effects, and the differences from the fluxes in figure 13 are small.

Figure 14. (a–d) Contours of mean streamwise velocity field averaged over control windows and turbine elements: (a,b) in a horizontal plane through the hub; (c,d) in a vertical plane through the turbine. (e,f) Contours of $\overline{\widetilde{u}}_{1}^{\prime \prime }\overline{\widetilde{u}}_{3}^{\prime \prime }$ in a vertical plane through the turbine. (a,c,e) Uncontrolled case; (b,d, f) the optimal control case averaged over time window $[5T_{A},20T_{A}]$ .

In order to further understand the increase of dispersive stresses and decrease of Reynolds stresses observed above, we first look at elements of the dispersive stress in figure 14 averaged in time and over the $10\times 5$ turbine subdomains (each with horizontal size $S_{x}\times S_{y}$ ). First of all, in figure 14(a,b) the mean streamwise velocity is shown in a horizontal plane at hub level. It is appreciated that the inflow velocities of the turbines in the controlled and uncontrolled cases are more or less equal. However, the wake velocity in the controlled case is much lower than for the uncontrolled case, which is clearly related to the fact that more energy is extracted in the controlled case. In spite of the lower wake velocity in the controlled case, the wake recovers faster than the uncontrolled wake. This is also visible in figure 14(c,d), where the streamwise velocity is shown in a vertical streamwise plane through the turbine centre. Here it is appreciated that, over the rotor height, the turbine inflow in the controlled case is even a bit higher than for the uncontrolled case. In figure 14(e, f) a vertical streamwise plane of $-\overline{\widetilde{u}}_{1}^{\prime \prime }\overline{\widetilde{u}}_{3}^{\prime \prime }$ is shown, which contributes to the dispersive shear stress (i.e.  $-\langle \overline{\widetilde{u}}_{1}^{\prime \prime }\overline{\widetilde{u}}_{3}^{\prime \prime }\rangle$ ) when averaged over horizontal planes. It is observed that positive $-\overline{\widetilde{u}}_{1}^{\prime \prime }\overline{\widetilde{u}}_{3}^{\prime \prime }$ regions increase in strength for the controlled case, while negative regions are not much altered. Given the distribution of the horizontal velocity in figure 14(c,d), it is clear that these positive $-\overline{\widetilde{u}}_{1}^{\prime \prime }\overline{\widetilde{u}}_{3}^{\prime \prime }$ regions are associated with a correlation of low mean streamwise velocity with upward mean motion. In the high-speed channels between the turbines, we also observe (not shown here) increased positive $-\overline{\widetilde{u}}_{1}^{\prime \prime }\overline{\widetilde{u}}_{3}^{\prime \prime }$ correlations in the controlled case, which is here associated with high mean streamwise velocity transported downwards by mean negative vertical motion.

In figure 15, the Reynolds stresses are further investigated in a turbine subdomain. In figure 15(a,b) the Reynolds shear stresses are shown in a horizontal plane through the turbine rotor tip (this is the region where the Reynolds stresses are largest in figure 12). We observe that, in the controlled case, the Reynolds shear stresses increase significantly above the turbine-wake region, which can explain the faster wake recovery observed in figure 14(b). On average, the increased effect of Reynolds shear stresses is compensated by lower shear stresses in the regions between the turbine rows when comparing controlled and uncontrolled cases, leading to lower horizontally averaged values (cf. figure 12 b). Also, in a vertical plane (figure 15 c,d), it is appreciated that Reynolds shear stresses increase significantly in the wake region for the controlled case. However, they do not increase in front of the turbine.

Figure 15. Contours of Reynolds stresses averaged over control windows and turbine elements: (a,b) horizontal plane at the turbine-tip level; (c–j)  $x$ – $z$ planes through the rotor centre. (a,c,e,g,i) Uncontrolled case; (b,d, f,h, j) the optimal control case averaged over time window $[5T_{A},20T_{A}]$ .

In figure 15(e–j), the normal Reynolds stresses are shown in a vertical plane. All normal stresses increase significantly in the wake region for the controlled case. When looking $1D$ upstream of the turbine, we observe that the streamwise stresses $\overline{\widetilde{u}_{1}^{\prime }\widetilde{u}_{1}^{\prime }}$ decrease in the controlled case compared to the uncontrolled case. This is possibly beneficial for reducing turbine loading and local blade stall (cf. also discussion at end of § 2.2). The spanwise stresses $\overline{\widetilde{u}_{2}^{\prime }\widetilde{u}_{2}^{\prime }}$ remain roughly unchanged in front of the turbine, while the vertical stresses $\overline{\widetilde{u}_{3}^{\prime }\widetilde{u}_{3}^{\prime }}$ slightly increase. Nevertheless, overall, compared to the uncontrolled case, the turbulent kinetic energy $\overline{\widetilde{u}_{i}^{\prime }\widetilde{u}_{i}^{\prime }}/2$ in the controlled case decreases by almost 9 % in front of the turbines (measured $1D$ upstream). We should remark that our simulations do not resolve turbulent fluctuations of the size of the turbine blade chord. It has been demonstrated experimentally that the energy exchange is dominated by turbulent scales with size of the rotor diameter (Hamilton et al. Reference Hamilton, Kang, Meneveau and Cal2012). However, smaller scales can be very relevant for the local blade performance, e.g. having an effect on local blade stall. In a classical turbulent energy cascade, it is expected that such smaller scales follow the larger scales (that are here resolved), but this has to be further established using better turbine representations (e.g. using actuator line models) and finer resolutions. This is subject of further research (cf. discussion in § 6).

Given the fact that the turbulence levels decrease in front of the turbines, the increased fluctuations in power output observed in figure 9 all result from fluctuations in the control $C_{T}^{\prime }$ . It is possible to make a Reynolds decomposition of the extracted power, using (2.12), and defining $C_{T,i}^{\prime }\equiv \overline{C_{T,i}^{\prime }}+{\it\Delta}[C_{T,i}^{\prime }]$ and $\widehat{V}_{i}^{2}V_{i}\equiv \overline{\widehat{V}_{i}^{2}V_{i}}+{\it\Delta}[\widehat{V}_{i}^{2}V_{i}]$ . Thus,

In the uncontrolled case, the second term on the right-hand side is zero. In the controlled case, we find that $C_{T}^{\prime }$ is slightly anticorrelated with $\widehat{V}_{i}^{2}V_{i}$ , leading to a negative value for the second term, with an observed magnitude that is approximately 6 % of the total extracted power $P$ (the first term on the right-hand side is 106 % of $P$ ). Consequently, the second term is a source in the turbulent kinetic energy equation. This may explain the locally increased turbulence levels observed in the turbine wakes above.

5.1. Adjoint solutions and optimal controls

In figure 16, a snapshot of an instantaneous adjoint field is shown for ${\it\gamma}=2/3$ . Different from the adjoint fields for the unpenalized case (cf. figure 7), in the current case, the initial condition for the adjoint equations differs from zero, and corresponds to ${\bf\xi}(\boldsymbol{x},T)=2/3\boldsymbol{u}(\boldsymbol{x},T)$ instead. This is visible in figure 16 early in the adjoint simulation.

Figure 16. Contours of instantaneous streamwise adjoint field for ${\it\gamma}=2/3$ , obtained from the first gradient calculation in control window 1: (a)  $T-t=14~\text{s}$ : (b)  $T-t=174~\text{s}$ . Horizontal planes are taken at the hub height.

Figure 17 shows the behaviour of the thrust coefficient for one of the turbines in the controlled wind farm. Both ${\it\gamma}=1/2$ and ${\it\gamma}=3/2$ cases show strong response to the turbulent flow field. However, in comparison to figure 8, the changes are less extreme. It can also be seen from this figure that $C_{T}^{\prime }$ mostly stays within the lower and upper limits, i.e. 0 and 4, although occasionally it still hits the upper bound.

Figure 17. Time evolution of the thrust coefficient of one of the turbines in the farm: (a)  ${\it\gamma}=1/2$ ; (b)  ${\it\gamma}=2/3$ .

5.2. Discussion of energy balances

Figure 18(a,b) show the time series of the total instantaneous gains and losses in the boundary layer for ${\it\gamma}=1/2$ and ${\it\gamma}=2/3$ , respectively. It can be appreciated for both cases that the overall deceleration of the boundary layer remains limited, and the optimal control case with ${\it\gamma}=2/3$ even slightly accelerates for $t>1000$ . Also, compared to figure 9, it is observed that now dissipation does not increase much in the controlled case.

Figure 18. Time evolution of gains and losses, for (a)  ${\it\gamma}=1/2$ and (b)  ${\it\gamma}=2/3$ : – – – (grey), driving power by pressure force; —— (grey), rate of change of kinetic energy; —— (black), farm power; - $\cdot$ - $\cdot$ - (black), dissipation.

In order to assess the precise gains and losses for the different cases, we assemble the time-integrated gains and losses in tables 2 and 3. Time averaging is performed from $[0,20T_{A}]$ , i.e. the time window for which the inner layer of the ${\it\gamma}=0$ case remains in equilibrium. The uncontrolled case is averaged over the same time window. First of all, in table 2, gains and losses are normalized with the total power extracted from the system by the wind farm and dissipation (i.e.  $\mathscr{P}+\mathscr{D}$ ). Looking at the results in the table, it is observed that the respective contributions of $\mathscr{P}$ and $\mathscr{D}$ to the total energy extraction from the boundary layer shift only slightly for the different cases. For the ${\it\gamma}=0$ case, which does not penalize dissipation, the ratio $\mathscr{P}/\mathscr{D}\approx 0.65$ deteriorates compared to the uncontrolled case for which $\mathscr{P}/\mathscr{D}\approx 0.68$ . For ${\it\gamma}=1/2$ and ${\it\gamma}=2/3$ , the ratios slightly improve, i.e.  $\mathscr{P}/\mathscr{D}\approx 0.72$ for both cases. We remark that for ${\it\gamma}=2/3$ the ratio is slightly lower than for ${\it\gamma}=1/2$ (cf. table 2), even though ${\it\gamma}=2/3$ penalizes dissipation more. However, given the limited averaging time, this difference is probably not significant. Moreover, the cost functional (5.1) does not directly optimize the ratio $\mathscr{P}/\mathscr{D}$ .

In table 3, the relative gains (changes) compared to the uncontrolled reference case are listed. It can be seen that all three cases increase wind-farm power extraction compared to the uncontrolled reference. Also dissipation increases for the three cases, but much less for ${\it\gamma}=1/2$ and ${\it\gamma}=2/3$ . Any increase in $\mathscr{P}+\mathscr{D}$ (compared to the reference) leads to an equal increase in the sum of power input and energy balance – the latter is also listed in table 3. Only for the optimal control case with ${\it\gamma}=2/3$ does the total power extracted from the system remain close to that of the uncontrolled reference. The gain in wind-farm power extraction in this case is limited to 6 %. Note that, if we were to assume a system that remains perfectly in equilibrium, then a change of $\mathscr{P}/\mathscr{D}$ from 0.68 to 0.72 (cf. above) would be equivalent to an increase of wind-farm power extraction of 5.8 %.

Finally, we remark that the order of magnitude of statistical errors in the discussion above is approximately 1 % – this is appreciated from $\text{d}E/\text{d}t\neq 0$ for the uncontrolled case in table 2.

5.3. Averaged flow profiles

Similar to the case without penalization, it is observed here that the inner region of the boundary layer is in statistical equilibrium after a short transition period. To assess averaged flow profiles in the current section, we will simply use the same averaging window $[5T_{A},20T_{A}]$ as proposed for ${\it\gamma}=0$ .

Average velocity profiles are shown in figure 19. Only the uncontrolled case and the averages over $[5T_{A},20T_{A}]$ are shown for ${\it\gamma}=1/2$ and ${\it\gamma}=2/3$ . We observe that the wind at turbine level decelerates a bit compared to the uncontrolled case, but now the outer layer slightly accelerates. This results from the fact that overall the driving power in these cases remains approximately constant (cf. figure 18), so that, with constant driving force, also the bulk velocity remains constant.

Figure 19. Streamwise mean velocities: ——, uncontrolled case; – – –, ${\it\gamma}=1/2$ ; - $\cdot$ - $\cdot$ -, ${\it\gamma}=2/3$ . The optimal control cases are averaged over the time window $[5T_{A},20T_{A}]$ .

In figure 20, we further look at total, Reynolds and dispersive stresses. In contrast to the unpenalized case, the peaks of total shear stresses in figure 20(a) (just above the wind farm) remain close to the uncontrolled case. For ${\it\gamma}=1/2$ it is a little higher, and for ${\it\gamma}=2/3$ a little lower, than the uncontrolled case, but this difference is statistically not significant in view of the limited temporal averaging time. Similar to before, dispersive stresses increase while Reynolds stresses decrease. When looking at the normal stresses in figure 20(b–d), some further differences are observed compared to optimal control with ${\it\gamma}=0$ . First of all, similar to ${\it\gamma}=0$ (cf. figure 12), all dispersive stresses increase compared to the uncontrolled case (though not as much as for ${\it\gamma}=0$ ). However, in contrast to ${\it\gamma}=0$ , the streamwise total stresses decrease compared to the uncontrolled case as a result of a significant decrease in streamwise Reynolds stresses. Looking at the vertical transport in the penalized cases (figure not shown here), we observe that the total transport remains largely unchanged, but is carried more by dispersive stresses (i.e. by mean flow transport) and less by turbulent stresses.

Figure 20. Vertical profiles of total, Reynolds and dispersive stresses: —— (black), —— (orange) and —— (cyan) respectively total stresses, Reynolds stresses and dispersive stresses for the uncontrolled case; – – – (black), – – – (orange) and – – – (cyan) respectively total stresses, Reynolds stresses and dispersive stresses for the controlled case with ${\it\gamma}=1/2$ averaged over $[5T_{A},20T_{A}]$ ; - $\cdot$ - $\cdot$ - (black), - $\cdot$ - $\cdot$ - (orange) and - $\cdot$ - $\cdot$ - (cyan) same for case with ${\it\gamma}=2/3$ .

When further analysing the spatial distribution of mean velocity profiles and Reynolds stresses (not further shown here), we observe again important changes between the controlled cases and the uncontrolled cases. Similar to ${\it\gamma}=0$ , the turbine inflow improves, with a slightly higher mean velocity, and lower turbulence intensity. However, the wake deficit is much less pronounced as less energy is extracted. Now, Reynolds stresses only slightly increase in the wake regions, but they significantly decrease between the turbine rows.