RANS modeling 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 experimental data from the literature are used for validation of the model.

1 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 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 RANS an attractive option for parametric studies or for isolating various physical effects (see van der Laan et al., 2021). The wind turbine forces are commonly represented as actuator disks (ADs) 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, ground effects and flow over complex geometries) are solved for directly rather than being prescribed through empirical relations. Disadvantages are that fatigue loading cannot be determined due to the steady nature of the method and 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 in stable, neutral and unstable conditions. The wake was found to recover faster (i.e., approach the freestream velocity at a shorter downstream distance) in the unstable ASL and slower (i.e., approach the freestream velocity at a longer downstream distance) in the stable ASL compared to in a neutral ASL, and this was later confirmed in field experiments by Magnusson and Smedman (1994) and full RANS simulations, including the temperature equation, by Alinot and Masson (2002). Rados et al. (2009) added a parametrized buoyancy term to the k–ε equations based on the MOST expressions, eliminating the need for a temperature equation. The “indirect” method of Rados et al. (2009) was shown by El-Askary et al. (2017) to produce similar wake deficit and turbulence intensity (TI) profiles as the “direct” method that employs a temperature equation.

In all the RANS studies discussed thus far, the simulations suffer from a known imbalance in the k and ε equations; this means that the freestream velocity and turbulence profiles vary horizontally throughout the domain, so that different wake results will be obtained depending on the streamwise position of the simulated turbine. van der Laan et al. (2017) solved this problem via the indirect method by adding analytical terms to the equations to be consistent with the ideal turbulent kinetic energy (TKE) budget under MOST and thus enforce a mean balance at all points. Han et al. (2019) used this approach in the direct method but did not show the extent to which their model is in balance.

Although there seems to be a general consensus that wakes should recover faster in unstable conditions, field measurements by Hansen et al. (2012) and Machefaux et al. (2016) found similar wake deficits for unstable and neutral conditions. This can possibly be attributed to the large uncertainties inherent in such measurement campaigns arising from sensors, post-processing and the unpredictable inflow provided by nature. In contrast, large eddy simulation (LES) offers a controlled environment, where complete statistics of all field variables can be extracted but at a large computational cost compared to RANS. Examples include Churchfield et al. (2012), Abkar and Porté-Agel (2015), Ghaisas et al. (2017) and Xie and Archer (2017), which simulate wakes in both unstable and neutral conditions for a wide variety of cases. All these studies agree with the general consensus and explain it with the increased TI encountered in unstable conditions due to buoyant production of turbulence. Nevertheless, both Alinot and Masson (2002) and Keck et al. (2014) show that a faster wake recovery in unstable conditions is still observed, even when the reference TI is kept fixed (by changing the roughness length); they argue that the enhanced wake recovery must be caused by the increased turbulent length scale associated with the unstable conditions because the turbulent velocity scale is approximately constant for fixed reference TI and wind speed.

The balanced k–ε MOST model by van der Laan et al. (2017) can be combined with the f_P correction, which was originally formulated by Apsley and Castro (1997) and later used by van der Laan (2014) to circumvent the over-diffusiveness of the k–ε model in the wake region. However, the f_P limiter was derived and calibrated for a neutral ASL, where it has been applied in many cases with success, but for non-neutral conditions it has yielded unphysical behavior, especially in unstable cases (van der Laan et al., 2021). Modifications to the MOST k–ε–f_P equations in the unstable regime are therefore suggested in this paper and validated against various field experiments and LESs.

2 Simulation setup

The wakes are simulated with the incompressible, finite-volume flow solver EllipSys3D (Michelsen, 1992; Sørensen, 1995). After continuous development since then, it is now a highly scalable code, which can be run in parallel on large high-performance clusters (HPCs) via a message passing interface (MPI). Thus a typical RANS simulation of a single wake only takes a few minutes to simulate on a contemporary HPC, while a similar case with LES would take several hours to run, even with an order of magnitude more computer resources available. In terms of CPU hours, van der Laan et al. (2015b) estimated LES to be approximately 10³ times more expensive than RANS, and this estimate may even be considered conservative because the LES inflow was created using a Mann-model turbulence box and not with the more expensive precursor method. Furthermore, several LES runs are in principle necessary in order to create an ensemble average, which multiplies the cost of LES. This 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.

2.1 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:

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathit{\zeta }\equiv {\displaystyle \frac{z}{L}},\end{array}$$

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{|}^{-\mathrm{1}}\mathit{\lesssim }\mathrm{0.002}$ m⁻¹ (e.g., Gryning et al., 2007) and tend to occur most often, with observed distributions of the stability ($\mathrm{1}/L$ or ζ) having a peak around zero (Kelly and Gryning, 2010). The most common unstable Obukhov lengths occur at −0.02 m${}^{-\mathrm{1}}\mathit{\lesssim }{L}^{-\mathrm{1}}\mathit{\lesssim }-\mathrm{0.002}$ m⁻¹ (Kelly and Gryning, 2010), but offshore there tends to be a bias towards more unstable conditions, i.e., more negative L⁻¹ 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

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }U={\displaystyle \frac{u_{*}}{\mathit{\kappa }}}\left[\ln \left({\displaystyle \frac{z}{z_{\mathrm{0}}}}\right)-{\mathrm{\Psi }}_m\right],V=W=\mathrm{0},\text{(3)}& {\displaystyle }& {\displaystyle }k={\displaystyle \frac{u_{*}^{\mathrm{2}}}{\sqrt{C_{\mathit{\mu }}}}}{\left({\displaystyle \frac{{\mathrm{\Phi }}_{\mathit{\epsilon }}}{{\mathrm{\Phi }}_m}}\right)}^{\mathrm{1}/\mathrm{2}},\text{(4)}& {\displaystyle }& {\displaystyle }\mathit{\epsilon }={\displaystyle \frac{u_{*}^{\mathrm{3}}}{\mathit{\kappa }z}}{\mathrm{\Phi }}_{\mathit{\epsilon }},\end{array}$$

where

$$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle }{\mathrm{\Psi }}_m=\ln \left[{\displaystyle \frac{\mathrm{1}}{\mathrm{8}}}\left(\mathrm{1}+{\mathrm{\Phi }}_m^{-\mathrm{2}}\right){\left(\mathrm{1}+{\mathrm{\Phi }}_m^{-\mathrm{1}}\right)}^{\mathrm{2}}\right]-\mathrm{2}\arctan \left({\mathrm{\Phi }}_m^{-\mathrm{1}}\right)+{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}},\text{(6)}& {\displaystyle }& {\displaystyle }{\mathrm{\Phi }}_m=(\mathrm{1}-\mathrm{16}\mathit{\zeta }{)}^{-\mathrm{1}/\mathrm{4}},\end{array}$$

and

$$\begin{array}{}\text{(7)}& {\mathrm{\Phi }}_{\mathit{\epsilon }}=\mathrm{1}-\mathit{\zeta }.\end{array}$$

The above relations are valid for $-\mathrm{2}\mathit{\lesssim }\mathit{\zeta }<\mathrm{0}$, so for a fixed L<0, it means that the equations are in principle only valid up to $z\approx -\mathrm{2}L$, where the free convection regime starts (i.e., buoyant production dominates over shear production of TKE). 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^{-\mathrm{1}}\mathit{\lesssim }-\mathrm{0.01}$ m⁻¹, 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.

The roughness length z₀ and friction velocity u_* in Eqs. (2) to (7) can be set using reference values (i.e., defined at z=z_ref) of wind speed (U_ref) and total TI (I_ref) along with the stability parameter (ζ_ref):

$$\begin{array}{}\text{(8)}& {\displaystyle }& {\displaystyle }{z}_{\mathrm{0}}=z_{\mathrm{ref}}\exp \left[-I_{\mathrm{ref}}^{-\mathrm{1}}{C}_{\mathit{\mu }}^{-\mathrm{1}/\mathrm{4}}\mathit{\kappa }\sqrt{{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}}{\left({\displaystyle \frac{{\mathrm{\Phi }}_{\mathit{\epsilon }}\left({\mathit{\zeta }}_{\mathrm{ref}}\right)}{{\mathrm{\Phi }}_m\left({\mathit{\zeta }}_{\mathrm{ref}}\right)}}\right)}^{\mathrm{1}/\mathrm{4}}-{\mathrm{\Psi }}_m\left({\mathit{\zeta }}_{\mathrm{ref}}\right)\right],\text{(9)}& {\displaystyle }& {\displaystyle }{u}_{*}=U_{\mathrm{ref}}{I}_{\mathrm{ref}}{C}_{\mathit{\mu }}^{\mathrm{1}/\mathrm{4}}\sqrt{{\displaystyle \frac{\mathrm{3}}{\mathrm{2}}}}{\left({\displaystyle \frac{{\mathrm{\Phi }}_m\left({\mathit{\zeta }}_{\mathrm{ref}}\right)}{{\mathrm{\Phi }}_{\mathit{\epsilon }}\left({\mathit{\zeta }}_{\mathrm{ref}}\right)}}\right)}^{\mathrm{1}/\mathrm{4}}.\end{array}$$

Note that TI (I) here is not the same as streamwise turbulence intensity (${\mathit{\sigma }}_u/U$); in this paper “TI” will refer to total TI (i.e., $I\equiv \sqrt{\frac{\mathrm{2}}{\mathrm{3}}k}/U$), unless stated otherwise. A typographical (sign) error has been corrected in Eq. (8), compared to the similar expression found in van der Laan et al. (2017). The von Karman constant κ and C_μ parameter are given with the other constants in Table 1.

Table 1 Parameters of the k–ε MOST turbulence model; see Sørensen (1995).

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., $\frac{\mathrm{1}}{A}\iint \frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }U(z{)}^{\mathrm{3}}\mathrm{d}A\approx \mathrm{308.9}$ W m⁻² (±0.6 %). One peculiarity 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 ($\mathrm{d}U/\mathrm{d}z$) because the velocity gradient scales with u_*, which scales with I_ref (Eq. 9); this is a consequence of specifying both TI and hub height velocity.

Figure 1 Analytical MOST profiles. Combinations of low/high TI and neutral/unstable stability. The rotor area of a NREL5MW turbine is shown. Dashed lines are used for ζ to make all profiles visible.

The eddy viscosity profile, ${\mathit{\nu }}_{\mathrm{t}}\left(z\right)=C_{\mathit{\mu }}\frac{k^{\mathrm{2}}}{\mathit{\epsilon }}$, is especially interesting to compare between the cases because ν_t features as 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 can be obtained by increasing either the turbulence strength (I_ref) or how unstable the atmosphere is (−ζ_ref) or both.

The eddy viscosity is sometimes expressed as a product of turbulent velocity and length scales (see Pope, 2000): ν_t=u_tℓ_t, where $u_{\mathrm{t}}=C_{\mathit{\mu }}^{\mathrm{1}/\mathrm{4}}{k}^{\mathrm{1}/\mathrm{2}}$ and ${\mathrm{\ell }}_{\mathrm{t}}=C_{\mathit{\mu }}^{\mathrm{3}/\mathrm{4}}{k}^{\mathrm{3}/\mathrm{2}}{\mathit{\epsilon }}^{-\mathrm{1}}$. These are plotted in Fig. 2 from which one can see that I_ref only affects u_t, while ℓ_t is unchanged. Changing ζ_ref mainly alters ℓ_t, while the hub height u_t and rotor-averaged u_t are unchanged. The increased ν_t (which gives faster wake recovery) due to unstable conditions can therefore not only be attributed to a larger I_ref, but also to an increased turbulent length scale caused solely by ζ_ref.

Figure 2 Turbulence scales in the freestream. Same labels as in Fig. 1.

2.2 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 modeled by a vortex extending from the center of the AD to infinity with constant circulation; hence, it resembles the classical Joukowsky optimum rotor (see 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; van der Laan et al., 2015a) is that only a few parameters are necessary: the thrust coefficient C_T, tip-speed 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-axisymmetric 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 behavior 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 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 or tower is included in our simulations, which have been shown to be a good approximation for >3 D downstream of the turbine according to Kasmi and Masson (2008) and Li et al. (2020).

2.3 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 $\mathrm{\Delta }x=\mathrm{\Delta }y=D/\mathrm{10}$ is used (based on the grid study in Appendix A1), while grid stretching is used in the vertical direction with Δz=z₀ at the first cell and $\mathrm{\Delta }z=D/\mathrm{10}$ at the cell at $z/D=\mathrm{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=\mathrm{25}$. Grid 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 $\approx \mathrm{0.45}\times {\mathrm{10}}^{\mathrm{6}}$ cells in the wake domain and $\approx \mathrm{2.1}\times {\mathrm{10}}^{\mathrm{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, while the drawback is that the majority of the cells are actually outside the region of interest, i.e., the wake domain.

Figure 3 Top and side views of the wake domain, the size of which is $\mathit{\{}{l}_x,l_y,l_z\mathit{\}}/D=\mathit{\{}\mathrm{16},\mathrm{4},\mathrm{3}\mathit{\}}$. The total grid size is $\mathit{\{}{L}_x,L_y,L_z\mathit{\}}=\mathit{\{}{l}_x+\mathrm{2}{x}_{\mathrm{in}},l_y+\mathrm{2}{x}_{\mathrm{in}},l_z+\mathrm{22}D\mathit{\}}$ and is too large to be shown here because x_in=5 km. Only every fourth cell is plotted.

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 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):

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }{U}_j{\displaystyle \frac{\partial k}{\partial x_j}}={\mathcal{D}}_k+\mathcal{P}-\mathit{\epsilon }+\mathcal{B}-S_k,\text{(11)}& {\displaystyle }& {\displaystyle }{U}_j{\displaystyle \frac{\partial \mathit{\epsilon }}{\partial x_j}}={\mathcal{D}}_{\mathit{\epsilon }}+\left(C_{\mathit{\epsilon }\mathrm{1}}\mathcal{P}-C_{\mathit{\epsilon }\mathrm{2}}\mathit{\epsilon }+C_{\mathit{\epsilon }\mathrm{3}}\mathcal{B}\right){\displaystyle \frac{\mathit{\epsilon }}{k}},\end{array}$$

where

$$\begin{array}{}\text{(12)}& {\displaystyle }& {\displaystyle }{\mathcal{D}}_{\mathit{\{}k,\mathit{\epsilon }\mathit{\}}}={\displaystyle \frac{\partial }{\partial x_j}}\left({\displaystyle \frac{{\mathit{\nu }}_{\mathrm{t}}}{{\mathit{\sigma }}_{\mathit{\{}k,\mathit{\epsilon }\mathit{\}}}}}{\displaystyle \frac{\partial \mathit{\{}k,\mathit{\epsilon }\mathit{\}}}{\partial x_j}}\right),\text{(13)}& {\displaystyle }& {\displaystyle }\mathcal{P}=-\stackrel{\mathrm{‾}}{u_i^{\prime }{u}_j^{\prime }}{\displaystyle \frac{\partial U_i}{\partial x_j}},\text{(14)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{‾}}{u_i^{\prime }{u}_j^{\prime }}={\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}k{\mathit{\delta }}_{ij}-{\mathit{\nu }}_{\mathrm{t}}\left({\displaystyle \frac{\partial U_i}{\partial x_j}}+{\displaystyle \frac{\partial U_j}{\partial x_i}}\right),\text{(15)}& {\displaystyle }& {\displaystyle }{\mathit{\nu }}_{\mathrm{t}}=C_{\mathit{\mu }}{f}_{\mathrm{P}}{\displaystyle \frac{k^{\mathrm{2}}}{\mathit{\epsilon }}},\text{(16)}& {\displaystyle }& {\displaystyle }{S}_k={\displaystyle \frac{u_{*}^{\mathrm{3}}}{\mathit{\kappa }L}}\left[{\mathit{\zeta }}^{-\mathrm{1}}\left({\mathrm{\Phi }}_m-{\mathrm{\Phi }}_{\mathit{\epsilon }}\right)-\mathrm{1}-{\displaystyle \frac{{\mathit{\kappa }}^{\mathrm{2}}}{\mathrm{4}{\mathit{\sigma }}_k\sqrt{C_{\mathit{\mu }}}}}{\mathrm{\Phi }}_m^{\mathrm{13}/\mathrm{2}}{\mathrm{\Phi }}_{\mathit{\epsilon }}^{-\mathrm{3}/\mathrm{2}}{f}_{\mathrm{un}}\right],\text{(17)}& {\displaystyle }& {\displaystyle }{f}_{\mathrm{un}}=(\mathrm{2}-\mathit{\zeta })+\mathrm{8}\left(\mathrm{1}-\mathrm{12}\mathit{\zeta }+\mathrm{7}{\mathit{\zeta }}^{\mathrm{2}}\right)-{\displaystyle \frac{\mathrm{1}}{\mathrm{16}}}\left(\mathrm{3}-\mathrm{54}\mathit{\zeta }+\mathrm{35}{\mathit{\zeta }}^{\mathrm{2}}\right),\text{(18)}& {\displaystyle }& {\displaystyle }{C}_{\mathit{\epsilon }\mathrm{3}}={\displaystyle \frac{\mathrm{1}}{\mathit{\zeta }}}\left(C_{\mathit{\epsilon }\mathrm{1}}{\mathrm{\Phi }}_m-C_{\mathit{\epsilon }\mathrm{2}}{\mathrm{\Phi }}_{\mathit{\epsilon }}+\left[C_{\mathit{\epsilon }\mathrm{2}}-C_{\mathit{\epsilon }\mathrm{1}}\right]{\mathrm{\Phi }}_{\mathit{\epsilon }}^{-\mathrm{1}/\mathrm{2}}{\mathrm{\Phi }}_m^{\mathrm{5}/\mathrm{2}}(\mathrm{1}-\mathrm{12}\mathit{\zeta })\right).\end{array}$$

The source term, S_k, and the C_ε,3 parameter constitute the two modifications compared to the usual k–ε equations (similar 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, ℬ, are deferred to the Sect. 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 van der Laan et al. (2017), with the only difference being that here we choose $\frac{{\mathrm{\Phi }}_h}{{\mathrm{\Phi }}_m{\mathit{\sigma }}_{\mathit{\theta }}}\to \mathrm{1}$; this modeler's choice for turbulent Prandtl number (σ_θ) avoids the inconsistency mentioned in that paper and makes the model independent of the temperature similarity function Φ_h.

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 Fig. 4 for an example with neutral inflow. As explained in Sect. 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.

Figure 4 Streamwise velocity (a, b) and kinematic eddy viscosity (c, d) are normalized by their freestream values, i.e., U_ref and ν_t,ref, respectively. The neutral I_ref=6 % case from Fig. 1 with a single NREL5MW turbine is shown.

Figure 5 TKE budgets of the 2017 model (a–c) and cstB model (d–f). The profiles are extracted at $y/D=\mathrm{0}$, i.e., in the center of the wake. ${\mathit{\zeta }}_{\mathrm{ref}}=-\mathrm{0.5}$ and I_ref=12 % for both rows.

The cause of the ν_t perturbation in the first place is the large velocity gradients across the AD and the shear layer, which enhances TKE shear production, but other terms of Eq. (10) are also highly active in these regions, and it is this complex interplay together with the f_P formulation that in the end determines the wake recovery. The effect of the buoyancy term in this interplay is discussed first and then afterwards the role of f_P in the unstable ASL.

3.1 Buoyant production term

The buoyant production of TKE is $\mathcal{B}\equiv \frac{g}{{\mathit{\theta }}_{\mathrm{0}}}\stackrel{\mathrm{‾}}{w^{\prime }{\mathit{\theta }}^{\prime }}$, and the heat flux is typically obtained using a temperature equation and a flux–gradient relationship. In this work, we pursue an alternative way and investigate two simple parametrizations:

$$\begin{array}{}\text{(19)}& {\displaystyle }& {\displaystyle }\mathcal{B}=-{\mathit{\nu }}_{\mathrm{t}}\left[{\left({\displaystyle \frac{\partial U}{\partial z}}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{\partial V}{\partial z}}\right)}^{\mathrm{2}}\right]{\displaystyle \frac{\mathit{\zeta }{\mathrm{\Phi }}_h}{{\mathit{\sigma }}_{\mathit{\theta }}{\mathrm{\Phi }}_m^{\mathrm{2}}}}\left(\mathrm{2017}\mathrm{model}\right),\text{(20)}& {\displaystyle }& {\displaystyle }\mathcal{B}=-{\displaystyle \frac{u_{*}^{\mathrm{3}}}{\mathit{\kappa }L}}\left(\mathrm{cstB}\mathrm{model}\right).\end{array}$$

The “2017 model” is the one utilized by van der Laan et al. (2017), van der Laan et al. (2020), Doubrawa et al. (2020) and van der Laan et al. (2021). This model does not require a temperature equation for closure but instead utilizes the temperature similarity function, Φ_h, and Prandtl number, σ_θ. The “cstB” model is as the name suggests simply a constant source term throughout the domain, and again no temperature equation is necessary. The k and ε equations are the same for the two models, except some minor changes to S_k, Eq. (16), and C_ε3, Eq. (18), are needed for the 2017 model; see Sect. 2.3. To isolate the effect of the ℬ parametrizations, they are tested here with f_P=1. A NREL5MW turbine (Jonkman et al., 2009) with U_ref=8 m s⁻¹ is used for all plots in this section.

Figure 6 Disk-averaged streamwise velocity, 〈U_AD〉, for low/high TI and neutral/unstable. ${\mathit{\zeta }}_{\mathrm{ref}}=-\mathrm{0.5}$ for unstable.

The upstream ($x/D=-\mathrm{10}$) budget in Fig. 5 shows the inconsistency of the 2017 model mentioned by van der Laan et al. (2017): the buoyant production goes to $\mathcal{B}=\frac{-u_{*}^{\mathrm{3}}}{\mathit{\kappa }L}{\mathrm{\Phi }}_m$, although $\mathcal{B}\to \frac{-u_{*}^{\mathrm{3}}}{\mathit{\kappa }L}$ is expected in the freestream, which can be derived from the ASL definition $\mathcal{B}\equiv \frac{g}{{\mathit{\theta }}_{\mathrm{0}}}{\stackrel{\mathrm{‾}}{{\mathit{\theta }}^{\prime }{w}^{\prime }}}_s$ and the Obukhov length definition. The cstB model on the other hand by definition complies with the freestream ASL limit of ℬ. Additionally, it can be noted that the cstB model has $\mathcal{B}/\mathcal{P}\approx \mathrm{1}$ at z_ref because ${\mathit{\zeta }}_{\mathrm{ref}}=-\mathrm{0.5}$ was used (see Fig. 5.23 of Stull, 1988).

A clear distinction between the two parametrizations is seen both in the near-wake ($x/D=\mathrm{1}$) and far-wake ($x/D=\mathrm{5}$) TKE budgets: in the top shear layer of the 2017 model simulation, large buoyant production is produced by the large velocity gradients in this region. This is observed in neither direct RANS simulations by El-Askary et al. (2017) nor in wind tunnel experiments by Zhang et al. (2013) or Hancock and Zhang (2015), so this may be deemed as an unphysical artifact. Indeed the 2017 model is derived for a homogeneous ASL, and applying it to a wind turbine wake violates this assumption. The cstB model on the other hand effectively uncouples the buoyant production and wake dynamics. This assumption can partly be justified by the aforementioned studies, which show that temperature changes very little in the wake from the ambient conditions and that the heat flux actually decreases in the wake.

Another deficiency of the unstable 2017 model is illustrated in Fig. 6: for a given I_ref, it unphysically predicts slower wake recovery than in neutral as also noted by van der Laan et al. (2021). This is remedied in the cstB model, where a slightly faster wake recovery is seen. It can be noted that in the near and far wake of the cstB model both ℬ and S_k terms are small compared to the other TKE terms (see Fig. 5), and as such it effectively resembles the neutral model, but with the one difference that it has a larger turbulent length scale (see Fig. 2), which explains the faster wake recovery seen in Fig. 6. The ℬ 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 neutral model (as demonstrated in Fig. 6) is also seen in the top row of Fig. 7. The second row shows the shear parameter, σ, which is a non-dimensional scalar metric describing the normalized velocity gradients; see Eq. (23). Both the magnitude and contours of σ are similar for the considered cases; i.e., it is large in the region around the AD and in the shear layers. The freestream shear parameter, $\stackrel{\mathrm{̃}}{\mathit{\sigma }}$, is however not equal for the two cases (see inlets in third row), which means that $\mathit{\sigma }/\stackrel{\mathrm{̃}}{\mathit{\sigma }}$ differs a lot between the two cases. This observation is important because $\mathit{\sigma }/\stackrel{\mathrm{̃}}{\mathit{\sigma }}$ is the main parameter used in the f_P correction to be discussed in the next section.

Figure 7 Normalized streamwise velocity (a, b), shear parameter (c, d) and normalized shear parameter (e, f). The cstB model is used for the unstable case (b, d, f), which has ${\mathit{\zeta }}_{\mathrm{ref}}=-\mathrm{0.5}$.

3.2 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. (15), which clearly affects the velocity deficit as shown in Fig. 4. The form of f_P used for wakes in the neutral ASL by van der Laan (2014) can be summarized as

$$\begin{array}{}\text{(21)}& {\displaystyle }& {\displaystyle }{f}_{\mathrm{P}}={\displaystyle \frac{\mathrm{2}{f}_{\mathrm{0}}}{\mathrm{1}+\sqrt{\mathrm{1}+\mathrm{4}{f}_{\mathrm{0}}\left(f_{\mathrm{0}}-\mathrm{1}\right)(\mathit{\sigma }/\stackrel{\mathrm{̃}}{\mathit{\sigma }}{)}^{\mathrm{2}}}}},\text{(22)}& {\displaystyle }& {\displaystyle }{f}_{\mathrm{0}}=\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{\left(C_{\mathrm{R}}-\mathrm{1}\right)}}\left(\mathrm{neutral}\right),\text{(23)}& {\displaystyle }& {\displaystyle }\mathit{\sigma }={\displaystyle \frac{k}{\mathit{\epsilon }}}\sqrt{{\left({\displaystyle \frac{\partial U_i}{\partial x_j}}\right)}^{\mathrm{2}}},\text{(24)}& {\displaystyle }& {\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\sigma }}={\displaystyle \frac{\mathrm{1}}{\sqrt{C_{\mathit{\mu }}}}}\left(\mathrm{neutral}\right).\end{array}$$

As shown in Fig. 7, $\mathit{\sigma }/\stackrel{\mathrm{̃}}{\mathit{\sigma }}$ is large in the region surrounding the rotor and in the shear layer (see Fig. 7), which leads to f_P<1 and hence a drop of ν_t in these regions. This is a desirable feature because it corrects the over-diffusiveness of the standard k–ε model.

It was recognized by van der Laan et al. (2020) that the freestream shear parameter, $\stackrel{\mathrm{̃}}{\mathit{\sigma }}$, has to be adjusted for MOST inflow in order to have f_P=1 in the freestream (the k–ε–f_P model should reduce to the standard k–ε model in the freestream):

$$\begin{array}{}\text{(25)}& \stackrel{\mathrm{̃}}{\mathit{\sigma }}={\displaystyle \frac{\mathrm{1}}{\sqrt{C_{\mathit{\mu }}}}}\sqrt{{\displaystyle \frac{{\mathrm{\Phi }}_m}{{\mathrm{\Phi }}_{\mathit{\epsilon }}}}}\left(\mathrm{general}\right).\end{array}$$

This can simply be derived by inserting the freestream profiles of U, k and ε (see Sect. 2.1) into the shear parameter definition, Eq. (23); the form of Eq. (25) has been used in all previous papers utilizing the 2017 model for wake modeling.

A more subtle modification arises recognizing that the f₀ parameter is also stability-dependent, i.e.,

$$\begin{array}{}\text{(26)}& f_{\mathrm{0}}=\mathrm{1}+{\displaystyle \frac{C_{\mathit{\mu }}{\stackrel{\mathrm{̃}}{\mathit{\sigma }}}^{\mathrm{2}}}{C_{\mathrm{R}}-\mathrm{1}}}\left(\mathrm{modification}\mathrm{1}\right).\end{array}$$

This is actually the form suggested by Apsley and Leschziner (1998), but they were not considering stability effects, i.e., no variation of $\stackrel{\mathrm{̃}}{\mathit{\sigma }}$ or C_R with stability; it has not been used in previous applications of the 2017 model. We note Eq. (26) is consistent with the neutral limit, since ${\stackrel{\mathrm{̃}}{\mathit{\sigma }}}^{\mathrm{2}}\to C_{\mathit{\mu }}^{-\mathrm{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

$$\begin{array}{}\text{(27)}& C_{\mathrm{R}}=\mathrm{4.5}+C_{\mathrm{B}}{\displaystyle \frac{\stackrel{\mathrm{̃}}{B}}{\stackrel{\mathrm{̃}}{\mathit{\epsilon }}}}\left(\mathrm{modification}\mathrm{2}\right),\end{array}$$

where $\stackrel{\mathrm{̃}}{B}/\stackrel{\mathrm{̃}}{\mathit{\epsilon }}$ is the freestream buoyant production to dissipation ratio and C_B is a new parameter to be calibrated. The effect of the stability-dependent $\stackrel{\mathrm{̃}}{\mathit{\sigma }}$, modification 1 and modification 2 is shown in Fig. 8. When plotted with $\mathit{\sigma }/\stackrel{\mathrm{̃}}{\mathit{\sigma }}$ on the abscissa, the stability-dependent $\stackrel{\mathrm{̃}}{\mathit{\sigma }}$, Eq. (25), has no effect compared to the neutral model (compare “Neutral” and “Unstable” in Fig. 8); this is clear from the definition of f_P in Eq. (21). The two modifications increase f_P when $\mathit{\sigma }/\stackrel{\mathrm{̃}}{\mathit{\sigma }}>\mathrm{1}$, i.e., in the wake. This is necessary to compensate for the larger $\mathit{\sigma }/\stackrel{\mathrm{̃}}{\mathit{\sigma }}$ encountered with non-neutral inflow; see Fig. 7.

Figure 8 Analytical form of the f_P correction as function of normalized shear parameter (f_P is independent of TI but depends on ζ; the unstable curves use ${\mathit{\zeta }}_{\mathrm{ref}}=-\mathrm{0.5}$ in this plot). The effect of stability-dependent $\stackrel{\mathrm{̃}}{\mathit{\sigma }}$ (Eq. 25), modification 1 (Eq. 26) and modification 2 (Eq. 27 with C_B=10.0) is shown.

When both modifications are used, faster wake recovery for a given I_ref occurs in unstable conditions, as shown in Fig. 9, i.e., similar behavior as when f_P was “turned off” (i.e., f_P=1); see Fig. 6.

Figure 9 Same as Fig. 6, but with f_P included. $C_{\mathrm{B}}=\mathit{\{}\mathrm{0.0},\mathrm{5.0},\mathrm{10.0}\mathit{\}}$ are tested for the unstable cstB model.

Table 2 Overview of test cases. SWiFT: Doubrawa et al. (2020). NTK41: Machefaux et al. (2016). V80-Abkar: Abkar and Porté-Agel (2015). V80-Keck: Keck et al. (2014). NREL5MW: Churchfield et al. (2012). The air density used for all cases is ρ=1.225 kg m⁻³.

4 Validation with experiments and LES

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_{\mathrm{B}}=\mathit{\{}\mathrm{0.0},\mathrm{5.0},\mathrm{10.0}\mathit{\}}$, while keeping the first term in Eq. (27) fixed to 4.5. The latter has been calibrated for a suite of neutral EllipSys3D LESs by van der Laan (2014), but in practice if it was calibrated with another LES code, a different, optimal value might have been obtained. In the same way, it cannot be expected that a universally valid C_B 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 Sect. 2.3; e.g., the cell size and extent of the wake region are scaled with the rotor diameter.

4.1 SWiFT case

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 Farm 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 meteorological mast located 2.5 D upstream of the turbine. Note that the stability parameter was measured to $\mathit{\zeta }=-\mathrm{0.089}$ at z=10 m, which at hub height corresponds to ${\mathit{\zeta }}_{\mathrm{ref}}=-\mathrm{0.29}$. Also, the streamwise turbulence intensity was measured at hub height to $I_{u,\mathrm{ref}}=\mathrm{12.6}$ %, which is converted to the total turbulence intensity as $I_{\mathrm{ref}}\approx \mathrm{0.8}{I}_{u,\mathrm{ref}}$ (van der Laan et al., 2015b). 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 {C_T, P and Ω} 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 3 D 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 LESs underpredict velocity deficit compared to the lidar measurements, which highlights the fundamental problem of comparing measurements with numerical models: even highly computational expensive simulations do not always match experimental results. This must be either due to experimental errors in the provided input data or because the idealizations used for the LESs are too simple to capture the wake behavior.

Figure 10 The unstable SWiFT case. The fixed frame of reference experimental and LES results were digitized from Doubrawa et al. (2020).

Figure 11 The unstable NTK41 case, where LES and experimental results were digitized from Machefaux et al. (2016). Profiles extracted at z=z_ref.

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 cancelations. 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 (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 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 EllipSys3D having been validated by van der Laan (2014) for numerous neutral cases. The very low wind veer of the stable SWiFT case (see Doubrawa et al., 2020) indicates that mesoscale effects were present during the experimental campaign.

4.2 NTK41 case

A Nordtank NTK41 500 kW 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 test case of this paper is based on lidar measurements and LESs conducted by Machefaux et al. (2016). They used two different models for their LESs; 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.

Figure 12 The unstable V80-Abkar case, where LES results were digitized from Abkar and Porté-Agel (2015). Profiles extracted at z=z_ref. Note that total TI (d–f) is based on U_ref and not the local velocity.

The NTK41 turbine is a stall-regulated wind turbine and is therefore operated at constant rotational speed independent of 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_{\mathrm{T},\mathrm{meas}}=\mathrm{0.71}$, while their LES gave $C_{\mathrm{T},\mathrm{LES}}=\mathrm{0.83}$. Looking up the thrust curve of the NTK41 turbine at U_ref=6.8 m s⁻¹ also gives $C_{\mathrm{T},\mathrm{curve}}=\mathrm{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 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 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).

4.3 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 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 modeled in their studies by an airfoil AD and operated at Ω=16.1 rpm. Neither C_T nor P was mentioned in the two papers, so in the following RANS simulations the values deduced from the power and thrust coefficient curves of Hansen et al. (2012) evaluated at U_ref=8 m s⁻¹ were used: C_T=0.81 and P=696 kW.

Relative to the LES, the velocity deficits shown in Fig. 12 are overpredicted by the 2017 model for all three downstream 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). The slower wake recovery of the 2017 model is consistent with less turbulence development in the near-AD region, which might explain its seemingly better TI prediction; however, we prioritize predicting a correct velocity deficit, since this is used for AEP calculations.

Figure 13 The unstable V80-Keck case, where LES results were digitized from Keck et al. (2014). Profiles extracted at z=z_ref. Note that streamwise TI (d–f) is based on U_ref and ${\mathit{\sigma }}_u\approx \sqrt{\frac{\mathrm{2}}{\mathrm{3}}k}/\mathrm{0.8}$ for the RANS simulations.

4.4 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 Churchfield et al. (2012), i.e., using a precursor simulation for the inflow and modeling the turbine with actuator lines (ALs). More specifically the unstable North Hoyle row A case is considered here; it features four V80 turbines spaced 11 D apart, using the inflow parameters described in the “V80-Keck” row of Table 2. Wake data are available at $x/D=\mathit{\{}\mathrm{4},\mathrm{5},\mathrm{6}\mathit{\}}$ 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 (see Table 2). The streamwise TI given by Keck et al. (2014) is converted to total TI with $I_{\mathrm{ref}}\approx \mathrm{0.8}{I}_{u,\mathrm{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 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 the fact 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.

4.5 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 7 D 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 ≥2 D upstream of the second turbine.

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. (2009), except that the thrust of Jonkman et al. (2009) also includes gravity and therefore cannot be used to obtain the aerodynamic thrust.

The velocity deficit of the new cstB model compares well with the LES data, especially so for C_B=5.0; see Fig. 14.

Figure 14 The unstable NREL5MW case, where LES results were digitized from Churchfield et al. (2012). Profiles extracted at z=z_ref.

5 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. For example comparisons with more detailed LES data (temperature profiles, shear production, buoyancy production, Reynolds stress tensor, anisotropy tensor, derived kinematic viscosity, etc.) would be useful.

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 analogous to C_R in the original f_P formulation. Both modifications 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 studies from the literature. Based on these comparisons, we recommend C_B=5.0 to be used but also acknowledge that each reference case was 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.

Testing the cstB model behavior in more complicated scenarios, i.e., aligned row cases (see Appendix A3), full wind farms, complex terrain, AEP calculations, etc., is a natural next step to map the applicability and limitations of the model. Also the extension to stable conditions for the cstB model is straightforward, but as the surface layer height in such conditions is small compared to modern turbine dimensions, it is questionable if the cstB model will be viable or if full ABL models are more appropriate. Yet another question to be answered is the effect of freestream turbulence anisotropy, which changes with stability and cannot be modeled with turbulence models based on the standard Boussinesq hypothesis.

Appendix A: Simulation details

A1 Grid study

Earlier studies by van der Laan et al. (2015b) have shown that a domain resolution of eight cells per 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 Fig. 3, and the Joukowsky AD (see next section) is used. Disk averages of velocity and TI are evaluated 1 D downstream of the turbine to verify grid independence in Fig. A1.

Figure A1 Grid-independence study. “Ref” means reference, i.e., the finest resolution available. Metrics are evaluated at $x/D=\mathrm{1}$. The cores used for the increasing domain resolutions are {45, 54, 63, 54, 60, 57}, respectively (non-constant because the domains are decomposed in different number of blocks).

Based on this small grid study, a domain resolution of 10 cells per diameter and an AD resolution of $(N_{\mathit{\phi }},N_{\mathrm{r}})=(\mathrm{32},\mathrm{32})$ are 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=\mathrm{5}$ the differences are only approximately 0.2 % for the velocity metric and 2 % for the TI metric). A simulation with this choice can be executed in about 1 wall clock minute on 63 cores (AMD EPYC 7351 processors are used).

A2 The Joukowsky AD method

In summary, the surface force distributions (unit: N m⁻²) on the AD are calculated in each iteration as

$$\begin{array}{}\text{(A1)}& {\displaystyle }\begin{array}{rl}{f}_{\mathrm{n},ij}=& {\displaystyle \frac{\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{C}_{\mathrm{T}}A{U}_{\mathrm{ref}}^{\mathrm{2}}}{F_{\mathrm{n}}^{\mathrm{0}}}}{f}_{\mathrm{n},ij}^{\mathrm{0}},f_{\mathrm{n},ij}^{\mathrm{0}}=\mathrm{4}\mathit{\rho }{q}_{\mathrm{0}}{\displaystyle \frac{g\left({\mathit{\chi }}_i\right)F\left({\mathit{\chi }}_i\right)}{{\mathit{\chi }}_i}}\\ & \left(\mathit{\lambda }{\mathit{\chi }}_i+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{q}_{\mathrm{0}}{\displaystyle \frac{g\left({\mathit{\chi }}_i\right)F\left({\mathit{\chi }}_i\right)}{{\mathit{\chi }}_i}}\right){\displaystyle \frac{U_{ij}^{\mathrm{2}}}{{\left(\mathrm{1}+\sqrt{\mathrm{1}-C_{\mathrm{T}}}\right)}^{\mathrm{2}}}},\end{array}\text{(A2)}& {\displaystyle }\begin{array}{rl}{f}_{\mathit{\theta },ij}=& {\displaystyle \frac{\frac{\mathrm{1}}{\mathrm{2}}\mathit{\rho }{C}_{\mathrm{P}}A{U}_{\mathrm{ref}}^{\mathrm{3}}}{P^{\mathrm{0}}}}{f}_{\mathit{\theta },ij}^{\mathrm{0}},f_{\mathit{\theta },ij}^{\mathrm{0}}=\mathrm{2}\mathit{\rho }{q}_{\mathrm{0}}{\displaystyle \frac{g\left({\mathit{\chi }}_i\right)F\left({\mathit{\chi }}_i\right)}{{\mathit{\chi }}_i}}\\ & {\displaystyle \frac{U_{ij}^{\mathrm{2}}}{\mathrm{1}+\sqrt{\mathrm{1}-C_{\mathrm{T}}}}},\end{array}\text{(A3)}& {\displaystyle }\begin{array}{rl}{q}_{\mathrm{0}}& ={\displaystyle \frac{\sqrt{\mathrm{16}{\mathit{\lambda }}^{\mathrm{2}}{a}_{\mathrm{2}}^{\mathrm{2}}+\mathrm{8}{a}_{\mathrm{1}}{C}_{\mathrm{T}}}-\mathrm{4}\mathit{\lambda }{a}_{\mathrm{2}}}{\mathrm{4}{a}_{\mathrm{1}}}},a_{\mathrm{1}}=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\displaystyle \frac{g^{\mathrm{2}}{F}^{\mathrm{2}}}{\mathit{\chi }}}\mathrm{d}\mathit{\chi },\\ & a_{\mathrm{2}}=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}gF\mathit{\chi }\mathrm{d}\mathit{\chi },\end{array}\text{(A4)}& {\displaystyle }g=\mathrm{1}-\exp \left[-a{\left({\displaystyle \frac{\mathit{\chi }}{\stackrel{\mathrm{‾}}{\mathit{\delta }}}}\right)}^b\right],\text{(A5)}& {\displaystyle }F={\displaystyle \frac{\mathrm{2}}{\mathit{\pi }}}\arccos \left[\exp \left(-{\displaystyle \frac{N_{\mathrm{b}}}{\mathrm{2}}}\sqrt{\mathrm{1}+{\mathit{\lambda }}^{\mathrm{2}}}\cdot (\mathrm{1}-\mathit{\chi })\right)\right].\end{array}$$

Figure A2 The normal blade loading, F_n [N m⁻¹], and tangential blade loading, F_t [N m⁻¹], are normalized by the density, ρ, rotor radius, R, and freestream hub height velocity, U_ref. The blade loadings for the Joukowsky AD and airfoil AD have been obtained by azimuthal averages, while $F_{\mathrm{n}}=\frac{\mathrm{2}\mathit{\pi }rT/A}{N_{\mathrm{b}}}$ is prescribed a priori for the uniform AD.

Figure A3 Disk average of velocity, 〈U_AD〉, and disk average of turbulence intensity, 〈I_AD〉.

Figure A4 Aligned row of 10 turbines with inflow similar to the Case 5 of van der Laan et al. (2021).

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. (2014) and Troldborg et al. (2015). $F_{\mathrm{n}}^{\mathrm{0}}\equiv \sum _i\left(\sum _j\left[f_{\mathrm{n},ij}^{\mathrm{0}}{A}_{ij}\right]\right)$ is the total normal force of the unscaled distribution, $P^{\mathrm{0}}\equiv U_{\mathrm{ref}}\mathit{\lambda }\sum _i\left({\mathit{\chi }}_i\sum _j\left[f_{\mathit{\theta },ij}^{\mathrm{0}}{A}_{ij}\right]\right)$ is the total power of the unscaled distribution, $\mathit{\lambda }\equiv \mathrm{\Omega }R/U_{\mathrm{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, ${\mathit{\chi }}_i\equiv 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₀ 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 $\stackrel{\mathrm{‾}}{\mathit{\delta }}=\mathrm{0.25}$) and Prandtl's tip correction, respectively; see Eqs. (A4) and (A5).

The normal and tangential loadings for the same case as used in Sect. A1 are compared between the uniform AD, airfoil AD and Joukowsky AD, the latter with two different $\stackrel{\mathrm{‾}}{\mathit{\delta }}$'s, in Fig. A2. Clearly, the Joukowsky AD with $\stackrel{\mathrm{‾}}{\mathit{\delta }}=\mathrm{0.25}$ produces similar loadings as the airfoil AD, e.g., qualitatively correct root behavior, tip behavior and constant tangential loading region, which makes it superior to 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_{\mathrm{ref}}=\mathrm{1}-\mathrm{0.5}\left(\mathrm{1}-\sqrt{\mathrm{1}-C_{\mathrm{T}}}\right)\approx \mathrm{0.74}$ at the rotor plane. This is not exactly observed in Fig. A3, but contrary to ideal 1D momentum theory our CFD simulation also includes atmospheric turbulence and shear.

A3 Aligned row case

Even though the focus of the present paper is on single wakes, wake simulations will in practice most often involve wake interaction. Contrary to engineering models, there is no need for empirical wake superposition methods in RANS, since the wake interaction automatically results from solving the RANS equations.

The Case 5 from the study of van der Laan et al. (2021), which consists of an aligned row of 10 NREL5MW turbines spaced 7 D apart, is simulated with a similar setup, and the normalized, disk-averaged velocity recovery is shown in Fig. A4. For reference, a neutral k–ε–f_P simulation with the same inflow speed and TI is simulated. The 2017 model performs poorly in the sense that it actually recovers slower than the neutral reference and that no equilibrium wake–wake interaction seems to occur. On the other hand, the cstB model predicts faster wake recovery for the first three to four turbines, while it goes to the same wake deficit as the neutral model in the fully developed part of the wind farm. More studies and validations with LES/experiments for similar wake interaction cases should be undertaken in the future.