Evaluation of tilt control for wind-turbine arrays in the atmospheric boundary layer

Wake redirection is a promising approach designed to mitigate turbine-wake interactions which have a negative impact on the performance and lifetime of wind farms. It has recently been found that substantial power gains can be obtained by tilting the rotors of spanwise-periodic wind-turbine arrays. Rotor tilt is associated to the generation of coherent streamwise vortices which deflect wakes towards the ground and, by exploiting the vertical wind shear, replace them with higher-momentum fluid (high-speed streaks). The objective of this work is to evaluate power gains that can be obtained by tilting rotors in spanwise-periodic wind-turbine arrays immersed in the atmospheric boundary layer and, in particular, to analyze the influence of the rotor size on power gains in the case where the turbines emerge from the atmospheric surface layer. We show that, for the case of wind-aligned arrays, large power gains can be obtained for positive tilt angles of the order of 30{\deg}. Power gains are substantially enhanced by operating tilted-rotor turbines at thrust coefficients higher than in the reference configuration. These power gains initially increase with the rotor size reaching a maximum for rotor diameters of the order of five boundary layer momentum thicknesses (for the considered cases) and decrease for larger sizes.Maximum power gains are obtained for wind-turbine spanwise spacings which are very similar to those of large-scale and very large scale streaky motions which are naturally amplified in turbulent boundary layers. These results are all congruent with the findings of previous investigations of passive control of canonical boundary layers for drag-reduction applications where high-speed streaks replaced wakes of spanwise-periodic rows of wall-mounted roughness elements.


Introduction
An unavoidable byproduct of wind-turbine operation is the generation of wakes containing the low-speed, highly turbulent fluid from which mechanical power has been extracted. As in wind farms the streamwise distance between wind turbines is typically much shorter than the distance required for the wake to diffuse and recover the incoming wind speed, wind-farm interior turbines experience large reductions of their power production and significant unsteady loads when they are shadowed by wakes of upstream turbines (see e.g. Stevens and Meneveau, 2017;Porté-Agel et al., 2019, for recent reviews). Among the significant number of techniques proposed to minimize the negative effects of turbine-wake interactions, the approach where wakes are deflected away from downstream turbines by yawing upstream rotors has recently attracted much interest. It has indeed been shown that power gains in downstream turbines associated to wake deflections can exceed power losses experienced by yawed turbines (Dahlberg and Medici, 2003;Jiménez et al., 2010).
In addition to yaw control, where wakes are deflected horizontally, it has been more recently proposed to deflect wakes in the vertical direction by acting on the rotor-tilt angle (see e.g. Guntur et al., 2012;Fleming et al., 2014). Best performances were obtained for positive tilt angles, where the wake is deflected towards the ground. This type of control, with positive tilt angles, is well suited for turbines with downwind-oriented rotor whose interest is being currently reconsidered because of their favorable distributions of blade bending loads, which are critical in the design of nextgeneration light and flexible blades, and their resilience in off-grid/extreme wind conditions (Loth et al., 2017;Kiyoki et al., 2017).
Also for control based on tilt, global power gains are achievable depite the power losses experienced by tiltedrotors upstream turbines (Fleming et al., 2015;Annoni et al., 2017;Bay et al., 2019;Cossu, 2020). In a companion paper (Cossu, 2020, referred to as C2020 in the following), we have recently shown that by tilting upstream rotors of spanwiseperiodic turbine rows, turbine wakes (low-speed fluid) can be replaced with high-speed streaks, i.e. streamwise elongated regions where the wind speed exceeds the mean, similarly to what observed in passive control approaches designed for drag-reduction applications (Fransson et al., 2006;Pujals et al., 2010b). An important aspect of the tilt-control approach is that it exploits the beneficial effect of the vertical wind shear, which is related to the lift-up effect by which streaks are amplified (Moffatt, 1967;Landahl, 1980;Schmid and Henningson, 2001) with a non-modal amplification mechanism (Böberg and Brosa, 1988;Butler and Farrell, 1992;Gustavsson, 1991;Trefethen et al., 1993;Cossu et al., 2009;Cossu and Hwang, 2017). Such a beneficial effect is not exploited by yaw control where the wakes are deflected in the horizontal plane and therefore tilt control can potentially lead to power gains larger than those associated to yaw control.
In C2020 it was found that significant improvements of the global power production can be obtained by operating tilted-rotor turbines at higher induction (higher thrust coefficient) to compensate for the reduction of the normal velocity component. It was also shown that those power gains could be further improved when considering higher values of the D/H ratio of the rotor diameter D to the boundary layer height H, at least for D/H < 0.36, in accordance with theoretical predictions of streak amplifications in turbulent wall-bounded shear flows Pujals et al., 2009;Hwang and Cossu, 2010;Cossu and Hwang, 2017) and with previous experimental results where streaks were forced with roughness elements instead of tilted-rotor wind turbines (White, 2002;Fransson et al., 2004;Pujals et al., 2010a). Global power gains of the order of 40% were attained for the largest considered ratio D/H = 0.36.
The results of C2020, nevertheless, call for confirmation in the high-D/H range because they do not necessarily extrapolate to atmospheric boundary layers since they were obtained in a pressure-driven boundary layer (PBL) where the effects of Coriolis acceleration and capping inversion were neglected, and furthermore with rotors taller than the logarithmic layer region mimicking the atmospheric surface layer. The scope of the present investigation is therefore to evaluate the level of power gains that can be obtained by tilting upstream rotors in wind-turbine arrays immersed in atmospheric boundary layers (ABL) where the effects of Coriolis acceleration and capping inversion are fully taken into account.
The questions in which we are interested are the following: What are the typical power gains that can be obtained by tilt-control in the ABL? How do they compare to those found in the PBL? How do these power gains depend on the relative rotor size, especially in the case of large rotors? To answer these questions, we use large-eddy simulations (LES) to simulate neutral atmospheric boundary layers with capping inversions in the presence of wind-turbine arrays in wind-aligned configurations which are the worst case for turbine-wake interactions. The turbines are modeled with the actuator-disk method and are assumed to operate at constant thrust coefficient in order to obtain generic results which do not depend of the specificity of particular control laws chosen for turbines operation. The effect of tilt angle, induction factor and rotor size on the global power production will be studied with respect to reference configurations where the turbines are operated in standard mode.
The paper is organized as follows: the problem formulation is introduced in §2, results are presented in §3 and summarized and further discussed in §4. Additional details on used numerical methods are provided in the appendix.

Problem formulation
We consider the flow developing on wind-turbine arrays immersed in the atmospheric boundary layer. The flow is computed by means of large-eddy simulations implemented in SOWFA, the NREL's Simulator for On/Offshore Wind Farm Applications which solves the filtered Navier-Stokes equations (see Churchfield et al., 2012, and app. A for additional details). The horizontal component of the Coriolis acceleration is included in the equations, compressibility effects are accounted for by means of the Boussinesq approximation and Smagorinsky (1963) model is used to approximate subgridscale stresses. It is assumed that Monin-Obukhov similarity theory for turbulent boundary layers above rough surfaces applies near the ground by implementing appropriate stress boundary conditions as in Schumann (1975) while slip boundary conditions are enforced at the upper plane of the solution domain. Periodic boundary conditions are applied in the spanwise direction with L y periodicity (the spanwise extension of the domain).
Preliminary 'precursor' simulations of the atmospheric boundary layer in the absence of wind turbines are run with periodic boundary conditions enforced also in the streamwise direction (with L x streamwise periodicity). The precursor simulations are used to generate realistic inflow wind conditions (Keating et al., 2004;Tabor and Baba-Ahmadi, 2010;Churchfield et al., 2012) for the simulations with wind turbines by storing the temporal evolution of flow variables on a vertical plane which is then used as inflow boundary condition for the simulations in the presence of wind turbines.
The actuator disk model (ADM) is used to approximate the forces exerted by wind turbines on the fluid and the power produced by the turbines. This model has been shown to provide reliable results for the characteristics of turbines wakes except in the wake formation region (Wu and Porté-Agel, 2011). To obtain general results, not depending on the specific control law assumed for the considered turbines, the total force exerted by each turbine on the fluid is assumed to be of the form F = −C T ρu 2 n Ae n /2, as done by e.g. Calaf et al. (2010), Goit and Meyers (2015) and Munters and Meyers (2017), where C T is the disk-based thrust coefficient, e n is the unit vector normal to the rotor disk, u n is the wind velocity component normal to the rotor averaged over the disk surface of area A = πD 2 /4 and D is the rotor diameter. The force is assumed to be uniformly distributed on the rotor and wake rotation effects are neglected because their modeling would reintroduce a dependence on specific turbine design and control. Turbines are assumed to (always) operate in Region II at constant C T .
The power produced by each turbine is P = C P ρu 3 n A/2 where C P = χC T with the coefficient χ = 0.9 accounting for wing-tip power losses and LES (coarse) grid effects (Martinez et al., 2016;Shapiro et al., 2019). Power gains that will be computed in the following, however, being ratios of powers computed with the same χ, do not depend on the specific value chosen for χ. The value C T = 2 corresponds to the optimal Betz value maximizing the power output of an isolated ideal turbine (Burton et al., 2001;Munters and Meyers, 2017).
We will consider the effect of rotor tilt in arrays composed of three rows of wind turbines aligned with the mean wind and spaced by λ x = 7D in the streamwise direction as Annoni et al. (2017), who considered a single column of three aligned turbines immersed in the ABL, and C2020 who studied spanwise-periodic arrays immersed in the pressure boundary layer (PBL). The arrays are periodic in the spanwise direction with spanwise turbine spacing λ y = 4D as in C2020, a value inspired by the spanwise spacing used for roughness elements used as 'streak generators' in previous boundary layer control studies as those of Fransson et al. (2004Fransson et al. ( , 2005Fransson et al. ( , 2006, Hollands andCossu (2009), Pujals et al. (2010b) and Pujals et al. (2010a). Tilt control will be applied to the two upwind rows of turbines by changing their rotortilt angle ϕ and thrust coefficient C T while turbines in the third, downwind, row are left at reference conditions.

Precursor simulations of the ABL
In the precursor simulations we consider neutral atmospheric boundary layers (ABL) at latitude 41 o N , whose thickness is kept almost steady by the presence of a capping inversion centered at z = H and of thickness ∆z CI , inside which there is a positive vertical potential temperature gradient (dθ/dz) CI . The potential temperature θ neut below the capping inversion is constant (300K), while a constant positive vertical gradient (dθ/dz) G = 0.03K/m is enforced in the geostrophic region above. The flow is driven by a pressure gradient enforcing mean westerly winds of 8m/s at  z = 100m. Three boundary layers are considered with H = 750m, H = 500m and H = 350m respectively, as detailed in Table 1. The precursor simulations are run up to t 1 = 20000s to extinguish initial transients and obtain a well-developed ABL. The driving mean pressure gradient and the flow variables on the west boundary (x = 0) are then stored from t 1 to t 2 = 30000s, i.e. for more than two and a half hours, to be used as driving term and inflow boundary conditions in the simulations with wind turbines. The flow is simulated in two sets of domains, the first extending 3km x 3km in the streamwise (x, west-east) and spanwise (y, north-south) directions respectively and the second of 6km x 3km extension all discretized with cells of 15m size in the horizontal directions and size ranging from 7m (near the ground) to 21m (in the freestream above the capping inversion) in the vertical direction.
The mean-wind profiles obtained for the three selected H in the 3km x 3km domains are nearly indistinguishable up to z ≈ 250m, where they are well fitted by the logarithmic law (see Fig. 1). In this region the wind direction is essentially directed along the x axis (from west to east) with a small y (north-south) component which is increasingly southerly for increasing heights. Identical results are found in the 6km x 3km domains (not shown).

Effect of tilt angle and C T on power gains
We first consider the effect of tilt on wind turbines with rotor diameter D = 126m and hub height z h = 89m (the same dimensions as the NREL 5-MW turbine model defined by Jonkman et al., 2009) in the 3km x 3km domain which accommodates three rows of six turbines each (see Fig. 3). The simulations with wind turbines are performed in the same numerical domain, with the same grid and parameters of the precursor simulations for 10000s using the driving mean pressure gradient and the west-plane (x = 0) inflow boundary conditions recorded in the precursor simulation. Statistics are accumulated for the last 6000s of the simulation. First, a reference case is simulated in the ABL with H = 750m where all turbine rotors have the usual −5 o tilt angle used to avoid blade-tower collisions for standard upwindaligned rotors (Jonkman et al., 2009) and are operated at constant C T = 1.5 which is consistent with those observed in real wind farms (Wu and Porté-Agel, 2015;Stevens et al., 2015;Munters, 2017). Then, a set of controlled cases are considered where the rotors of the two most upwind turbine rows are tilted all by the same angle ϕ (as in Annoni et al., 2017, andC2020).
The ratios of the global power produced in the controlled case to the global power produced in the reference case are reported in Fig. 2a for the selected tilt angles ϕ = 20 o , 30 o and 40 o . From Fig. 2a it is seen that, for the present case where all turbines are operated at the same C T = 1.5, best power gains above 15% are obtained for ϕ = 30 0 . These gains are of the same order of those found by Annoni et al. (2017) for three aligned turbines in the ABL (using the specific NREL5 load distribution and control law) and those found by C2020 in the pressure boundary layer (PBL) for the same array configuration and operation.
Following C2020, the simulations are repeated operating tilted-rotor turbines at the higher C T = 3, which is still in the range attainable by real turbines (see e.g. Munters and Meyers, 2017). From Fig. 2a we see that in this case the power gains have almost doubled, attaining values above 30% for ϕ = 30 o , similarly to what previously found in the PBL.
From Fig. 2b, where relative powers extracted by each row of turbines are reported, it can be seen that the effect of tilt is to moderately increase the power extracted by the middle row of turbines and greatly increase the power extracted by the most downwind row. The effect of tilt, indeed, is to create a pair of counter-rotating streamwise vortices (Dahlberg and Medici, 2003;Howland et al., 2016;Bastankhah and Porté-Agel, 2016) which are horizontally-staked (see Fig. 4b) and are associated with a strong (vertical) downwash in the wake of the tilted rotors (see Fig. 4b and the blue regions of negative vertical velocity clearly discernible in Fig. 3d). The downwash induced by the first and second row of turbines is almost additive (see Fig. 3d). From panels a and b of Fig. 3 and from Fig. 4 it can be seen that the forced downwash is strong enough to expel the low-speed region (the wake) first towards the ground and then on the sides, and replace it with high-speed streaks by exploiting the positive wind shear (Fleming et al., 2015;Annoni et al., 2017;Cossu, 2020).
These results, obtained in the ABL, are consistent with those found in the pressure boundary layer (PBL), as could be expected because for the considered D = 126m turbines and the H = 750m ABL (D/H = 0.126) the rotors are mostly immersed in the atmospheric surface layer, where the velocity profile is logarithmic (see Fig. 1a) and which is well approximated by the inner region of the PBL. In the considered case, the slight wind veering with height has for effect only a slight bending of the wake region and of the highspeed region in the cross-stream plane (see Fig. 4) with no apparent significant influence on the coherent lift-up process by which streaks are amplified Pujals et al., 2009;Cossu and Hwang, 2017;Cossu, 2020).

Influence of the relative rotor size on power gains
In C2020, for the case of the pressure boundary layer (PBL), increasing power gains were found for increasing D/H ratios of the rotor diameter to the boundary layer thickness, at least up to D/H = 0.36. As already mentioned, however, for the largest values of D/H considered in C2020 the results obtained in the PBL are not necessarily a good approximation of those that one would find in the ABL because rotors emerged from the logarithmic region where PBL and ABL mean wind profiles are similar. Furthermore, in the PBL considered in C2020 a slip boundary condition at z = H was enforced, which could also have artificially affected streak growths for the largest considered D/H values. We therefore investigate if the results found in the PBL hold also in  the ABL. In the present ABL setting, where H is given by the capping inversion height, we will also be able to access values of D/H larger than in C2020. We begin by considering the effect of reducing H from H = 750m to H = 500m and then to H = 350m for the same array of turbines with D = 126m considered in §3.2 (where tilted-rotor turbines are operated at C T = 3 while the others are operated at the reference C T = 1.5). Reducing H at constant D results in increasing the D/H ratio. From Fig. 5a, where the results are reported, it is seen that for all considered H, the maximum power gains are obtained for ϕ = 30 o and that for this tilt angle they do increase with D/H but begin to saturate for the largest considered D/H values.
To better explore the high-D/H regime, additional simulations are performed with larger turbines keeping constant the relative spacing of 4D and 7D in the spanwise and streamwise directions respectively and the operation mode (C T = 3 for tilted rotors, C T = 1.5 for the others). First, an array of turbines with D = 180m and z h = 115m (inspired by the DTU-10 MW model) is considered in the same 3km x 3km domains already used (three rows with four turbines each are accommodated in the domain). Then, arrays of even larger turbines with D = 250m, z h = 150m and D = 360m, z h = 230m are accommodated in the 6km x 3km domains (three rows with 3 and 2 turbines respectively in each row in the simulated domain). Also for the D = 180m turbines (as shown in Fig. 5b) and for the D = 250m and D = 360m turbines (not shown) the best power gains are obtained for ϕ = 30 o , a value which is therefore robust to the change of D and H. However, from Fig. 5b it is also seen that for D = 180m turbines, the best power gains are obtained in the H = 500m ABL (D/H = 0.36) and not in the H = 350m one corresponding to the largest value D/H = 0.51. Similar results are found for D = 250m and D = 360m turbines.
To better appreciate the influence of H and D on power gains, the P/P ref values obtained for ϕ = 30 o are reported as a function of D/H in Fig. 6 for all the considered D-H combinations and are compared to those obtained in the PBL (reproduced from data reported in C2020). From Fig. 6 it is seen that, for all the considered boundary layers, power gains initially increase with turbine diameter (and therefore D/H) before reaching a maximum value (D/H ≈ 0.24 − 0.51 depending on the considered H) and eventually decrease when D is further increased.
The saturation of power gains at sufficiently large D/H can not be attributed to lateral wake/streaks displacements associated to strong wind-veer effects because, for all considered ABLs, saturation is reached for the D = 180 turbines, which remain in the region where mean velocity profiles are identical for the three considered H and wind veer is negligible (see Fig. 1). The saturation can probably be associated to the existence of an optimal spanwise wavelength where the amplification of the streamwise streaks induced by the streamwise vortices generated by the tilted rotors is maximum, as theoretical models of optimal energy amplifications predict Cossu et al., 2009;Hwang and Cossu, 2010) and as is observed in experiments of forced streaks amplifications in turbulent boundary layers (Pujals et al., 2010b). Following this line of thought, the maximum of P/P ref , obtained for D/H ≈ 0.24 − 0.51, i.e. for a spanwise wavelength of λ y = 4D ≈ 1 − 2.5H, would therefore correspond to the spanwise wavelength of maximum amplification of the streamwise streaks leading to the largest mean velocity increase in correspondence to downstream rotors. From Fig. 6 it is also seen that P/P ref data are quite dispersed when plotted as a function of D/H and, most importantly, maximum power gains are attained at values of D/H which strongly vary from one boundary layer to another. This could have been expected, since the considered turbulent boundary layers are not self-similar. Previous investigations (see e.g. Corbett and Bottaro, 2000;Cossu et al., 2009) had, however, shown that the most relevant lengthscale selecting optimal spanwise wavelengths λ y associated to maximum streaks amplification in non-self-similar boundary layers is the boundary layer momentum thickness where U ∞ is the streamwise velocity in the freestream above the boundary layer. We therefore replot in Fig. 7 the P/P ref data as a function of D/δ 2 . From Fig. 7 it is seen that this new rescaling (on δ 2 instead of H) greatly reduces the data dispersion and that maximum gains are obtained for D/δ 2 ≈ 4.3 − 5. This enhanced scaling based on δ 2 provides further evidence that, also in the ABL, streaks amplification is the main driving mechanism of the observed power gains obtained by tilting rotors.

Conclusions
In this study we have evaluated the gains of global extracted wind power that can be obtained by tilting rotors in wind-turbine arrays immersed in neutral atmospheric boundary layers with capping inversions. In particular, we have considered spanwise-periodic arrays with three turbine rows where rotors of the two upwind rows are all tilted by the same angle ϕ while rotors of the downwind row are left in reference conditions.
It is found that significant power gains can be obtained. For all the considered combinations of turbine rotors diameters (D = 126m, D = 180m, D = 250m and D = 360m), boundary layer heights (H = 750m, H = 500m and H = 350m), and C T of operation, the best power gains are obtained for tilt angles ϕ ≈ 30 o . Power gains are substantially improved when tilted turbines are operated at induction rates (corresponding to C T = 3) higher than the reference one (corresponding to C T = 1.5).
The influence of rotor size on power gains has also been investigated for ratios D/H up to 0.72, larger than those previously considered by Cossu (2020) in the PBL. Results show that power gains do initially increase with rotor size (starting from small D/H values) but they eventually saturate and decrease for sufficiently large ratios D/H.
Following the rationale of Cossu (2020), relating tiltinduced power gains to the amplification of large-scale streaks forced by tilted rotors, the existence of an optimal rotor diameter maximizing tilt-induced power gains has been connected to the existence of an optimal spanwise wavelength maximizing streak amplifications. Credit to this interpretation is given by the fact that the range of optimal spanwise turbine spacings λ y = 4D ≈ 1 − 2.5H corresponds to the spacing of large-scale streaks associated to large-scale and very-large-scale streaky motions (LSM and VLSM) in neutral turbulent boundary layers (see e.g. Hutchins et al., 2012;Shah and Bou-Zeid, 2014;Fang and Porté-Agel, 2015) whose natural amplification is based on the same mechanism Cossu et al., 2009;Hwang and Cossu, 2010;Cossu and Hwang, 2017). Further credit to this interpretation is given by the fact that power gain data obtained for all the considered D-H turbine-ABL combinations do scale better on D/δ 2 than on D/H ratios, where δ 2 is the boundary layer (streamwise) momentum thickness, as do streak-energy amplifications in non-self-similar boundary layers (Corbett and Bottaro, 2000;Cossu et al., 2009).
Additional final evidence is, nonetheless, needed to confirm the relation between streaks amplification and tiltinduced power gains in the ABL by comparing streamwise vortices forced by tilted rotors to the optimal perturbations leading to maximum large-scale streaks amplification in the ABL and in particular their optimal spanwise spacing. Optimal perturbations and energy amplifications in the ABL are, however, currently unknown (they have, actually, been computed in the Ekman boundary layer by Foster, 1997, but not for ABL profiles with logarithmic profiles in the surface layer). This is the subject of current intensive research effort.
It is also important to emphasize that a non-negligible part of the obtained power gains is associated to the oper-ation of tilted-rotor turbines at higher induction rates. In the present study and in the previous related one of Cossu (2020) only a few basic combinations of tilt angles ϕ and (steady) thrust coefficients C T have been considered for tilted-rotor turbines. It is believed that there is large room for further power-gain improvements that would come from the optimization of the ϕ − C T distribution among controlled turbines. Further work is also needed to evaluate power gains attainable for typical wind roses instead of the single windaligned case considered here.
Moreover, additional investigations, based on more realistic turbine models and control laws, are needed to assess the level of practically attainable power gains and their dependence of the chosen control laws. More realistic turbine models are especially needed in high-D/H regimes where rotors are large enough to directly interact with the capping inversion and the geostrophic free-stream. Such regimes would not only be encountered for future-generation very large turbines but also in current-generation wind turbines immersed in stable ABLs where boundary layer heights H are small. In this case it is known that wake rotation effects are important (see e.g. Englberger et al., 2020) and wind veer and Coriolis acceleration are likely to strongly influence the streak amplification process, especially trough a modification of the trajectory of the pair of counter-rotating streamwise vortices generated by the tilt.

Appendix A: Methods
The filtered Navier-Stokes equations, including the effect of Coriolis acceleration, with Boussinesq fluid model and the Smagorinsky (1963) model of turbulent subgrid scale motions are solved with SOWFA (Churchfield et al., 2012), an open-source code based on the OpenFOAM finite-volumes code (OpenCFD, 2011). The PIMPLE scheme is used for time advancement. Schumann (1975) stress boundary conditions are enforced at the near-ground horizontal boundary.
Numerical simulations have been performed in numerical domains of streamwise length L x (3km and 6km) and width L y (3km). Three domains heights L z have been considered as detailed in Table 1. Totally, six numerical domains have therefore been used in the simulations. All the domains are discretized with cells of constant 15m size in the horizontal directions and height-increasing size ranging from the 7m to 21m in the vertical direction. The solutions are advanced with ∆t = 0.8s time steps which keeps reasonable the amount of data stored in precursor simulations.
In order to enforce the turbine response to depend only on C T and not on the turbine controller, for the sake of simplifying the interpretation of the results, an in-house ADMC model has been derived from the original actuator disk (ADM) turbine model implemented in SOWFA (Martinez-Tossas and Leonardi, 2013) by distributing the body force uniformly in the radial direction while keeping the same dis-cretization points on the disk and setting its magnitude to F = −C T ρu 2 n Ae n /2 (Calaf et al., 2010;Goit and Meyers, 2015;Munters and Meyers, 2017) and setting to zero body force components parallel to the rotor plane, thus setting to zero the wake rotation. A Gaussian projection of discretized body forces with a smoothing parameter ε = 20m is also used to avoid numerical instabilities (Martinez-Tossas and Leonardi, 2013).