RANS modelling of a single wind turbine wake in the unstable surface layer

Unstable atmospheric conditions are often observed during the daytime over land and for significant periods offshore, and are hence relevant for wake studies. A simple k-ε RANS turbulence model for simulation of wind turbine wakes in the unstable surface layer is presented, which is based on Monin-Obukhov similarity theory (MOST). The turbulence model parametrizes buoyant production of turbulent kinetic energy (TKE) without the use of an active temperature equation, and flow balance is ensured throughout the domain by modifications of the turbulence transport equations. Large eddy simulations and 5 experimental data from the literature are used for validation of the model.


Introduction
Wind turbine wakes have been studied for decades using many different methodologies, including wind tunnels, field experiments, analytical engineering models, and numerical simulations. A review of these methodologies is given by Porté-Agel et al. (2020) and it is noteworthy that many of the references therein are from the past decade. The motivation for many of these new 10 studies is the large number of new wind farms emerging each year, where wake effects significantly impact the Annual Energy Production (AEP), as well as wind farm lifetime through increased fatigue.
A sub-category of "numerical simulations" is the Reynolds-Averaged Navier Stokes (RANS) approach, which is a Computational Fluid Dynamics (CFD) method that solves for the mean fields. This means that no time history of the flow is obtained, however the computational resources required for RANS are very small compared to higher-fidelity CFD methods, making 15 RANS an attractive option for parametric studies or for isolating various physical effects (c.f. van der Laan et al., 2021). The wind turbine forces are commonly represented as Actuator Disks (AD) in RANS; several types of AD models are reviewed by van der Laan et al. (2015a). Compared to engineering models, an advantage of RANS is that physical features of the flow (e.g. induction zones, wake interaction, shear layers and ground effects) are solved for directly, rather than being prescribed through empirical relations. Disadvantages are that fatigue loading can not be determined due to the steady nature of the method, and 20 that the solution relies heavily on the turbulence model.
The part of the atmosphere closest to the ground, i.e. the atmospheric surface layer (ASL), can be parametrized with the similarity theory of Monin and Obukhov (1954) (MOST) and used as inflow for RANS simulations of wind turbine wakes. The k-ε turbulence model is usually preferred in RANS wake studies and Crespo et al. (1985) for example utilized this (although in a parabolized RANS setup, which requires less computational resources, but is less accurate) to simulate a single wake clearly motivates the development of the RANS model as a fast, albeit less accurate, alternative to LES.
The different components of the RANS simulations will be discussed in the following sections.

Inflow profile for unstable ASL
Numerous articles have been written about MOST and a historical review is given by Foken (2006). The theory is expressed and applied via the dimensionless stability parameter where z is the height above ground and L is the Obukhov length. Negative ζ corresponds to unstable conditions, while ζ = 0 corresponds to the neutral limit where there is no effect of buoyancy. Neutral conditions are typically defined as |L| −1 0.002 m −1 (e.g. Gryning et al., 2007) and tend to occur most often, with observed distributions of the stability (1/L or ζ) having a peak around zero (Kelly and Gryning, 2010). The most common unstable Obukhov lengths occur at −0.02 m −1 75 L −1 −0.002 m −1 (Kelly and Gryning, 2010); but offshore, there tends to be a bias towards more unstable conditions, i.e. more negative L −1 compared to onshore (Sathe et al., 2013). Various parametrizations have been suggested for wind speed, k, and ε in terms of ζ; in this paper we use the widely accepted forms of Dyer (1974) for U (namely the Ψ m and Φ m functions), and those found in Kaimal and Finnigan (1994) for ε and k (see van der Laan et al., 2017, for more details). Under unstable conditions these are: Although the blade tip of a modern turbine can reach beyond −2L in unstable conditions (e.g. for z tip = 200 m this happens when L −1 −0.01 m −1 , which is not rare), we nevertheless still choose to apply the profiles-and in fact use them all the way up to the top boundary. More realistic inflow profiles for RANS covering the whole Atmospheric Boundary Layer (ABL) are indeed a current research topic, but will not be discussed further in this paper. Maronga and Reuder (2017) reason that MOST is a "pragmatic solution," because the parameters needed for more realistic inflow profiles are often not available.

95
In this paper κ = 0.4 and C µ = 0.03 are used, while z 0 and u * in Eqs.
(2) to (8) Note that TI (I) here is not the same as streamwise turbulence intensity (σ u /U ); in this paper "TI" will refer to total TI (i.e., 100 I ≡ 2 3 k/U ), unless stated otherwise. A typographical (sign) error has been corrected in Eq. (9), compared to the similar expression found in van der Laan et al. (2017).
Examples of four inflow profiles with identical U ref are shown in Fig 1. The stability and TI differ among the cases, but they still have approximately the same averaged power density; i.e., 1 A 1 2 ρU (z) 3 dA ≈ 308.9 W m −2 (±0.6 %). One peculiarity 105 of the unstable MOST profiles is that the TI does not go to zero for z → ∞; this is connected to the ASL assumptions used by MOST (i.e., the mixed and upper layers above the ASL are lacking, where I becomes constant and then vanishes). Another peculiarity (for both neutral and unstable conditions) is that higher TI leads to larger shear (dU/dz), because the velocity gradient scales with u * , which scales with I ref (Eq. 10); this is a consequence of specifying both TI and hub height velocity.
The eddy viscosity profile, ν t (z) = C µ k 2 ε , is especially interesting to compare between the cases, because ν t features as 110 a diffusion coefficient in the Reynolds-Averaged momentum equation and therefore is connected with the entrainment of ambient air into the wake. A faster wake recovery is therefore expected for larger ν t , which as seen in Fig. 1   The eddy viscosity is sometimes expressed as a product of turbulent velocity and length scales, c.f. Pope (2000): where u t = C 1/4 µ k 1/2 and t = C 3/4 µ k 3/2 ε −1 . These are plotted in Fig. 2, from which it is clear that TI only affects u t (l t is 115 unaffected due to u * cancelling when dividing k 3/2 from Eq. (3) by ε from Eq. (4)), while stability mainly alters t .

Wind Turbine representation
A recently developed Actuator Disk (AD) model by Sørensen et al. (2020) is utilized in this paper. The model can be derived from conservation of energy (Bernoulli's equation), conservation of angular momentum (Euler's turbine equation) and an analytical expression for the near-wake azimuthal velocity distribution. The latter is modelled by a vortex extending from the 120 center of the AD to infinity with constant circulation, hence it resembles the classical Joukowsky optimum rotor, c.f. Okulov and Sørensen (2010), and the AD model is therefore referred to as the "Joukowsky-AD". A summary of the model formulation is given in Appendix A2.
The main advantage of the Joukowsky-AD over the widely used "airfoil-AD" (e.g. Sørensen and Kock (1995), Porté-Agel et al. (2011) andvan der Laan et al. (2015a)) is that only a few parameters are necessary: The thrust coefficient C T , tip-speed 125 ratio λ, rotor radius R and freestream reference wind speed U ref (in addition to these, the power coefficient C P is also made an input parameter in our implementation as described by van der Laan et al. (2020)). Nevertheless, it is still able to model non-axissymmetric force distributions and wake rotation, similar to the disk loading of an airfoil-AD. Porté-Agel et al. (2011) argued that these are important features to capture the correct wake behaviour in the near-wake, while van der Laan et al.
(2015c) showed that they are only of minor importance for the far-wake. Wake deficit and rotor loading of the Joukowsky-AD 130 have been found to compare well with several validation cases conducted by Sørensen et al. (2020) and Sørensen and Andersen (2020). This is verified to also be the case for our RANS simulations in Appendix A2.
No nacelle nor tower are included in our simulations, which have been shown to be a good approximation for > 3D downstream of the turbine according to Kasmi and Masson (2008) and Li et al. (2020).

RANS
A homogeneous, flat lower surface is assumed for all cases in this paper. The inner part of the mesh surrounding the AD is called the "wake domain" and is shown for a typical case in Fig. 3. In this area, a horizontal resolution of ∆x = ∆y = D/10 is used (based on the grid study in Appendix A1), while grid stretching is used in the vertical direction with ∆z = z 0 at the first cell and ∆z = D/10 at the cell at z/D = 3. The wake domain is however only a small part of the full domain: The full domain extends an additional x in = 5 km to the west, north, east and south, respectively, while the top of the grid is at z/D = 25. Grid 140 stretching is used in all directions outside of the wake domain to circumvent an excessively large number of cells (the case shown in Fig. 3 has ≈ 2.1 · 10 6 cells in total). The choice of having grids with size on the order of ∼ 10 km is made to avoid tunnel-like blockage effects, and to have fully developed inflow profiles at the turbine position. The numerical solution strategy of the incompressible RANS equations in EllipSys3D is thoroughly discussed in other publications (Michelsen, 1992;Sørensen, 1995;Sørensen et al., 2007), so only the main features are discussed here. The

145
SIMPLE method is used with a modified Rhie-Chow algorithm, following Réthoré (2009) and Troldborg et al. (2015), to avoid the numerical wiggles induced by the discrete actuator disk body forces.
As mentioned in the introduction, the flow variables in an empty domain with MOST inflow can be kept in balance by modifying the k-and ε-equations, as suggested by van der Laan et al. (2017): where The source term, S k , and the C ε,3 parameter constitute the two modifications compared to the usual k-ε equations (similar 160 corrections exist for the stable ASL, but are not discussed in this paper). Viscous terms have been neglected in the above equations, which is a good approximation in atmospheric flow applications due to the Reynolds number being very large (Wyngaard, 2010). The Coriolis force is also neglected, hence no veer is present in the simulations. Definitions of the f P -correction (which was in fact set to f P = 1 in the work of van der Laan et al. (2017)) and buoyant production, B, are deferred to the Section 3. The parameters used in the above equations are summarized in Table 1. Finally, S k and C ε,3 differ slightly from those printed in 165 van der Laan et al. (2017), with the only difference being that here we choose Φ h Φmσ θ → 1; this "modeller's choice" for turbulent Prandtl number (σ θ ) avoids the inconsistency mentioned in that paper, and makes the model independent of the temperature 3 Modification of the k-ε-f P model in the unstable ASL The background eddy viscosity shown in Fig. 1 is perturbed in the turbine presence and especially so in the wake region, see 170 Fig. 4 for an example with neutral inflow. As explained in Section 2.1, ν t is very important for the wake development and the f P -correction effectively attenuates the ν t perturbation in the interface between the wake and freestream, known as the shear layer, and in the region around the AD to improve wake predictions. This attenuation can also be viewed as a modification of the turbulence scales, u t and t .
At z = z hub The cause of the ν t perturbation in the first place is the large velocity gradients across the AD and the shear layer, which 175 enhances TKE shear production, but other terms of Eq. (11) are also highly active in these regions, and it is this complex interplay together with the f P -formulation that in the end determine the wake recovery. The effect of the buoyancy term in this interplay is discussed first, then afterwards the role of f P in the unstable ASL.

Buoyant production term
The buoyant production of TKE is B ≡ g θ0 w θ and the heat flux is typically obtained using a temperature equation and a 180 flux-gradient relationship. In this work, we pursue an alternative way and investigate two simple parametrizations: The "2017 model" is the one utilized by van der Laan et al. (2017)     wake recovery is seen. It can be noted that in the near-and far-wake of the cstB model both B and S k terms are small compared to the other TKE terms, c.f. Fig 5, and as such it effectively resembles the neutral model, but with the one difference that it has a larger turbulent length scale, c.f. Fig. 2, which explains the faster wake recovery seen in Fig. 6. The B and S k terms must nevertheless still be retained to enforce the freestream balance of the k-and ε-equations throughout the domain.
The faster wake recovery of the unstable cstB model compared to the similar neutral case is also seen in Fig. 7 along with 210 the shear parameter and turbulence time scale, which are both used in the f P -correction to be discussed in the next section.

Turbulence closure with f P in non-neutral conditions
As stated in the introduction, k-ε models tend to predict faster wake recovery compared to experiments and LES. This can be corrected by using f P = 1 in the ν t definition, Eq. (16), which clearly affects the veloctiy deficit as shown in Fig. 4 (2014), can be summarized as: The shear parameter, σ, is large in the region surrounding the rotor and in the shear layer, see Fig. 7, which explains the drop 220 of ν t in these regions.
It was recognized by van der Laan et al. (2020), that the freestream shear parameter,σ, has to be adjusted for MOST inflow in order to have f P = 1 in the freestream: This can simply be derived by inserting the freestream profiles of U , k and ε (see Sect. 2.1) into the shear parameter definition,

225
Eq. (24); the form of Eq. (26) has been used in all previous papers utilizing the 2017 model for wake modelling.
A more subtle modification arises recognizing that the f 0 parameter is also stability-dependent, i.e., This is actually the form suggested by Apsley and Leschziner (1998), but they were not considering stability effects, i.e. no variation ofσ nor C R with stability; it has not been used in previous applications of the 2017 model. We note Eq. (27) is 230 consistent with the neutral limit, sinceσ 2 → C −1 µ for ζ → 0. The Rotta constant was calibrated to C R = 4.5 for wind turbine wakes in the neutral ASL in the work of van der Laan (2014) and we therefore require C R → 4.5 in the neutral limit (ζ → 0).

One form that satisfies this is
whereB/ε is the freestream buoyant production to dissipation ratio and C B is a new parameter to be calibrated. The effect of 235 modifications 1 and 2 is shown in Fig. 8. It is clear that the two modifications increase f P in regions where σ/σ > 1, i.e. in the wake. This is necessary because larger σ/σ are encountered with non-neutral inflow, sinceσ decreases and k/ε increases in unstable conditions (note from Eq. (24) that σ ∼ k/ε); see Fig. 7.
When both modifications are used, faster wake recovery for a given I ref occurs in unstable conditions, as shown in Fig. 9; this was also the case when f P was 'turned off' (set to 1), c.f. Fig. 6.  The cstB model with the f P -modifications described in the previous section is tested for the cases summarized in Table 2. Each case is simulated with a range of C B parameters, C B = {0.0, 5.0, 10.0}, while keeping C R = 4.5 fixed. The latter has been calibrated for a suite of neutral EllipSys3D LES's by van der Laan (2014), but in practice if it was calibrated with another LES code, a different, optimal C R might have been obtained. In the same way, it cannot be expected that a universally valid C B 245 exists, when we compare with results from a range of different codes and experiments. Hence, no optimal C B will be obtained in this section, but rather the qualitative effect of C B is shown.
The numerical setup for each case follows that described in Section 2.3, e.g. the cell size and extent of the wake region are scaled with the rotor diameter.

Case
Type

SWiFT case 250
A large wake benchmark study was conducted by Doubrawa et al. (2020) to compare various simulation methodologies and codes against LIDAR measurements in different atmospheric conditions. The measurements were carried out for a Vestas V27 turbine at the Scaled Wind Technology Facility (SWiFT) in Lubbock, Texas, USA, which is an area of flat terrain.
The inflow parameters of the SWiFT row in Table 2 were obtained from the ensemble average of five 10 min-averages from a met. mast located 2.5D upstream of the turbine. Note, that the stability parameter was measured to ζ = −0.089 at This conversion could actually be slightly different in the unstable ASL, because the ratios of velocity variance change with stability (e.g. Chougule et al., 2018), but unfortunately only the vertical velocity variance follows MOST, so that no general surface layer formula can be constructed (Panofsky and Dutton, 1984;Wyngaard, 2010). The operational state parameters 260 {C T , P , Ω} were taken from the OpenFAST steady-state curves, which were supplied for the benchmark.
For the unstable SWiFT case, the wake profile was only measured at 3D downstream and the results are shown in Fig. 10.
Three different LES codes were used in the benchmark and the filled area in Fig. 10 represents the spread of the LES results.
It can be seen that all three LES's underpredict velocity deficit compared to the LIDAR measurements, which highlights the fundamental problem of comparing measurements with numerical models: Even highly computational expensive simulations 265 do not always match experimental results. This must either be due to experimental errors in the provided input data or because the idealizations used for the LES's are too simple to capture the wake behaviour.
RANS can generally not be expected to perform better than a well-performed LES and if it does, it is likely due to fortunate error cancellations. Therefore, from a theoretical point of view one could argue that the performance of RANS should mainly be assessed with regards to how well it matches the LES results. Both RANS and LES use many of the same idealizations 270 (uniform roughness, flat terrain, homogeneous inflow, etc.) and indeed our RANS results in Fig. 10 are also closer to the LES results, than to the experimental results.
For the SWiFT case the 2017 model seemingly performs better than the cstB model, but this is probably due to fortuitous model error and/or some unaccounted mesoscale effects. Model error was expected, because the neutral EllipSys3D RANS simulation in the SWiFT study (surprisingly) did not compare well with experimental results (Doubrawa et al., 2020), despite

NTK41 case
A Nordtank NTK41 500kW wind turbine was installed at what is now the Risø campus of the Technical University of Denmark in 1992, and was used for many research studies before its decommissioning in 2021. Among these studies, the "NTK41 280 testcase" of this paper is based on LIDAR measurments and LES's conducted by Machefaux et al. (2016). They used two different models for their LES's; the results included in Fig. 11 (along with our RANS results) are from their more advanced model, which they called the "LES-ABL" or "extended model." The inflow parameters ("NTK41" row in Table 2) and LIDAR measurements were ensemble-averaged over 20 10-min averages.
The NTK41 turbine is a stall-regulated wind turbine and is therefore operated at constant rotational speed independent of 285 the inflow wind speed (Hansen, 2015), in this case at Ω = 27.1 rpm. The thrust coefficient for the unstable case of Machefaux et al. (2016) was measured with strain-gauges to be C T,meas = 0.71, while their LES gave C T,LES = 0.83. Looking up the thrust curve of the NTK41 turbine at U ref = 6.8 m/s also gives C T,curve = 0.83, so this will be used in the present RANS simulations. They argued that the lower thrust coefficient of their measurement could be explained by the large uncertainty of the strain gauges. Finally, the measured power was P meas = 120 kW, while P LES = 127 kW and P curve = 125 kW, where 290 the latter will be used to set C P for the AD model of the RANS simulations. Figure 11 shows that the cstB model matches the LES and experimental data better than the 2017 model, although a still faster wake recovery is seen for both of these reference data. Compared to more conventional LES setups (e.g. V80-Abkar, V80-Keck and NREL5MW cases), the LES model used by Machefaux et al. (2016) is simplified by using a modified Mann box for its inflow, and a slip condition at the bottom wall; together these add some uncertainty to the LES results. The experimental 295 wake data can also be expected to have large uncertainties and/or biases connected to them, but the sources and sizes of these were not discussed in detail by Machefaux et al. (2016).

V80-Abkar case
Abkar and Porté-Agel (2014) investigated the effect of atmospheric stability using LES and used the results to modify the analytical Bastankhah wake model (Abkar and Porté-Agel, 2015). They studied the wake of a single Vestas V80 turbine (known 300 from e.g. the Horns Rev 1, North Hoyle and Princess Amalia wind farms), which has been used in many previous wake studies.
The turbine was modelled in their studies by an airfoil-AD and operated at Ω = 16.1 rpm. Neither C T nor P were mentioned in the two papers, so in the following RANS simulations the values deduced from the power and thrust coefficient curves of Relative to the LES, the velocity deficits shown in Fig. 12 are overpredicted by the 2017 model for all three downstream 305 distances, while the cstB model corrects this especially well in the far-wake. Besides the velocity deficit, the TI based on the freestream velocity was also available for this case and is plotted in the lower row of Fig. 12 with the RANS results. The wake TI is overpredicted by RANS, which is also typically seen for neutral RANS simulations (van der Laan, 2014).

V80-Keck case
This case is based on a LES from Keck et al. (2014), where the SOWFA solver was used with a similar setup as in the work of 310 Churchfield et al. (2012), i.e. using a precursor simulation for the inflow and modelling the turbine with Actuator Lines (AL).
More specifically the "unstable North Hoyle row A" case is considered here; it features four V80 turbines spaced 11D apart, using the inflow parameters described in the "V80-Keck" row of Table 2. Wake data is available at x/D = {4, 5, 6} downstream of the first turbine and the induction effect of the downstream turbines on this data should therefore be minimal, hence they are omitted from the RANS simulations. The inflow wind speed is U ref = 8 m/s, the same as in the V80-Abkar case, thus the same operational state of the wind turbine is utilized, c.f. Table 2. The streamwise TI given by Keck et al. (2014) is converted to total TI with I ref ≈ 0.8I u,ref , similar to the method used in the SWiFT case. Figure 13 shows that the cstB model improves the velocity deficit prediction over the old 2017 model, when comparing with the LES results, which were digitized from Fig. 11 in Keck et al. (2014) (note that a typo is present in that figure, i.e. label should be "[D]" instead of "[R]"). Streamwise TI was also available at the same downstream distances and in the RANS 320 simulations it was obtained by converting from the total TI as described above. In both the streamwise TI and velocity deficit LES data, a misalignment of the wake center can be observed, which Keck et al. (2014) explains with that only 10 min of LES data was averaged. This is especially visible in the streamwise TI plots, but nevertheless it seems that the cstB model predicts streamwise TI in the right range except for an overprediction in the wake center, which was also seen in the previous case.

NREL5MW case
The last case is based on the LES studies by Churchfield et al. (2012), more specifically their "U-L case" (see inflow parameters in the "NREL5MW row" of Table 2). They model two NREL5MW turbines with the Actuator Line (AL) methodology coupled to the aeroelastic FAST solver and the turbines are spaced 7D apart. For the RANS simulations of the present study, we omit the second turbine and only compare with the first wake of the LES study. To avoid biases from the induction zone of the second turbine, we only consider wake results ≥ 2D upstream of the second turbine.

330
The steady-state power, thrust coefficient and rotational speed were not given in the paper, so therefore the steady-state curves from the DTU in-house aeroelastic solver, HAWCStab2, were used. These are similar to the curves shown by Jonkman et al.  Churchfield et al. (2012). Profiles extracted at z = z ref .

Conclusions
We have proposed a simple k-ε RANS model, the "cstB" model, for simulation of wind turbine wakes in the unstable surface layer. The model does not require an additional temperature equation and instead bases the TKE buoyancy production on MOST and the assumption that it is decoupled from the wake dynamics, which means that the buoyant production of TKE is constant throughout the domain, even in the wake region, hence the name "cstB". Wind tunnel studies and simulations have hinted that the latter assumption is reasonable, but a more thorough investigation would be beneficial for developing simple, non-neutral wake models.
Originally developed to account for the general over-diffusive nature of k-ε models in wind turbine wakes under neutral conditions, here the f P -correction is combined with the new cstB model by making two non-neutral modifications. These introduce a new parameter, C B ; it is a free parameter analagous to C R in the original f P -formulation. Both modifications 345 are consistent in the sense that the new, non-neutral f P -formulation becomes equal to the original neutral f P form for ζ → 0.
By using this updated f P -model with the cstB model, a faster wake recovery is obtained for unstable conditions over neutral conditions, when TI is fixed, as was also the case when no f P -model was applied.
The cstB model with the modified f P -function was generally found to perform better than the previous model of van der Laan et al. (2017) with the old f P -formulation, in terms of velocity deficit profiles from five different reference cases found in 350 studies from the literature. Based on these comparisons, we recommend C B = 5.0 to be used, but also acknowledge that each reference case were originally conducted with different numerical and experimental setups, and that further studies are needed to conclude on a more certain C B value, which could also be slightly code-dependent, as has been seen for C R in the original f P -model.

A1 Grid study
Earlier studies by van der Laan et al. (2015b) have shown that a domain resolution of eight cells pr. diameter is sufficient to obtain grid independence for wakes in the neutral ASL. A range of domain and AD resolutions are here tested for the new cstB model and f P -modifications with the TI = 12 % and C B = 5.0 case also used in Fig. 9. The domain size is described in   Based on this small grid study, a domain resolution of 10 cells pr. diameter and an AD resolution of (N ϕ , N r ) = (32, 32) is chosen for the current investigation. The difference between this resolution and the reference resolution is less than 2 % for the velocity metric and less than 6 % for the TI metric (the differences decrease with downstream position, e.g. at x/D = 5 the differences are only approximately 0.2 % for the velocity metric and 2 % for the TI metric). A simulation with this choice can 365 be executed in about one wallclock minute on 63 cores (AMD EPYC 7351 processors are used).

A2 The Joukowsky AD method
In summary, the surface force distributions (unit: N/m 2 ) on the AD are calculated in each iteration as: Here, f n,ij and f θ,ij are the normal and azimuthal surface force distributions at the (i, j)'th AD element (i: radial direction, j: azimuthal direction), which are applied in the CFD domain using the methodology described by Réthoré et al. 375 (2014) and Troldborg et al. (2015). F 0 n ≡ i j f 0 n,ij A ij is the total normal force of the unscaled distribution, P 0 ≡ U ref λ i χ i j f 0 θ,ij A ij is the total power of the unscaled distribution, λ ≡ ΩR/U ref is the tip-speed ratio, A is the area of the AD, A ij is the area of the (i, j)'th AD element, χ i ≡ r i /R is the local normalized radius, U ij is the normal velocity at the (i, j)'th AD element, N b = 3 is the number of blades, q 0 is the normalized circulation, g is the root correction and F is the tip correction. The latter two are obtained with Delery's root correction (with parameters a = 2.335, b = 4.0 and δ = 0.25) and 380 Prandtl's tip correction, respectively, c.f. Eq. (A4) and (A5).
The normal and tangential loadings for the same case as used in Section A1 are compared between the uniform-AD, airfoil-AD and Joukowsky-AD, the latter with two different δ's, in Fig. A2. Clearly, the Joukowsky-AD with δ = 0.25 produces similar loadings as the airfoil-AD, e.g. qualitatively correct root behaviour, tip behaviour and constant tangential loading region, which 385 makes it superior to the uniform-AD. is prescribed a-priori for the uniform-AD.
Finally, in Fig. A3 the velocity and TI disk averages follow the same trend for all four AD methods, but with a slightly larger velocity deficit for the airfoil-AD, possibly because of its also slightly larger blade loadings, see Fig. A2. The thrust coefficient of the uniform-AD and Joukowsky-AD is C T = 0.77, which from 1D momentum theory should give U/U ref = 1 − 0.5 1 − √ 1 − C T ≈ 0.74 at the rotor plane. This is not exactly observed in Fig. A3, but contrary to ideal 1D momentum 390 theory our CFD simulation also includes atmospheric turbulence and shear. Code availability. EllipSys3D and pywake_ellipsys are propritary software of DTU. Information about the latter can however be freely accessed at https://topfarm.pages.windenergy.dtu.dk/cuttingedge/pywake/pywake_ellipsys/index.html.
Data availability. The RANS results were generated with DTU's proprietary software, but the data presented can be made available by contacting the corresponding author. Interested parties are also welcome to hand-digitize the results and use them as reference in other