Wind turbine wake simulation with explicit algebraic Reynolds stress modeling

Reynolds-averaged Navier–Stokes (RANS) simulations of wind turbine wakes are usually conducted with two-equation turbulence models based on the Boussinesq hypothesis; these are simple and robust but lack the capability of predicting various turbulence phenomena. Using the explicit algebraic Reynolds stress model (EARSM) of Wallin and Johansson (2000) can alleviate some of these deficiencies while still being numerically robust and only slightly more computationally expensive than the traditional two-equation models. The model implementation is verified with the homogeneous shear flow, half-channel flow, and square duct flow cases, and subsequently full three-dimensional wake simulations are run and analyzed. The results are compared with reference large-eddy simulation (LES) data, which show that the EARSM especially improves the prediction of turbulence anisotropy and turbulence intensity but that it also predicts less Gaussian wake profile shapes.

1 Introduction

As wind farms increase in size and number of turbines, increasingly more attention should be given to the study of wind turbine wakes, as they can account for a relatively large power production decrease. Simple engineering models, such as the classic Jensen (1983) model, the “new-classic” Bastankhah and Porté-Agel (2014) model, or the more recent Ishihara model (Ishihara and Qian, 2018), can be used to model the flow through a wind farm on a regular laptop in a matter of seconds. However, they all share some common weaknesses: they are based on rather strict flow assumptions, need empirically tuned parameters, and require a superposition model for overlapping wakes.

A wind turbine wake is a complex three-dimensional swirling flow, and its development is governed by turbulent mixing, which is strongly influenced by density stratification in the atmospheric boundary layer (ABL) and the interaction between the ABL and the wake itself. To model a more physically correct wind farm flow we therefore solve the Navier–Stokes equations, which in essence are a reformulation of Newton's second law, along with conservation of mass. The process of discretizing and solving these equations on computers is known as computational fluid dynamics (CFD). Unfortunately, the Reynolds number of atmospheric flows is so large (10⁶–10⁸ or more; see, e.g., Wyngaard, 2010) that it is unfeasible to conduct direct numerical simulation (DNS) for even a single-wind-turbine wake. Instead large-eddy simulation (LES), i.e., the simulation of the spatially filtered Navier–Stokes equations, has been performed by numerous groups in the last decade; see review by Breton et al. (2017). Even though the rotor geometry is reduced to an actuator disk (AD) or an actuator line (AL) model, these simulations require a vast number of CPU hours on modern high-performance computing (HPC) clusters, again even for a single-wind-turbine wake. A much faster simulation can be conducted by solving the Reynolds-averaged Navier–Stokes (RANS) equations instead; however, in contrast to LES, which needs only to parameterize turbulent fluctuations at the smallest (unresolved) scales, RANS simulation relies heavily on the quality of the turbulence model because it must model fluctuations across all simulated scales. For example, the standard two-equation k–ε model of Launder and Spalding (1974) was shown by Réthoré (2009) to perform poorly for simulation of wind turbine wakes in terms of wake recovery prediction and also generating unphysical Reynolds stresses in the vicinity of the turbine (i.e., unrealizable turbulence). This is not unexpected since it was originally developed for simple free shear and boundary layer flows, e.g., flat plate, pipe, plane jet, and cavity flows (Launder and Spalding, 1974). The standard k–ε model is considered to be a linear eddy viscosity model (EVM) because the anisotropy tensor, $a_{ij}\equiv \frac{\stackrel{\mathrm{‾}}{u_i^{\prime }{u}_j^{\prime }}}{k}-\frac{\mathrm{2}}{\mathrm{3}}{\mathit{\delta }}_{ij}$ (where $\stackrel{\mathrm{‾}}{u_i^{\prime }{u}_j^{\prime }}$ is the Reynolds stress tensor; k is the turbulent kinetic energy, TKE; and δ_ij is the Kronecker delta tensor), is linearly proportional to the normalized strain rate tensor through the Boussinesq hypothesis (Boussinesq, 1897); see also Eq. (5). A modified model coined the “k–ε–f_P” model was developed and calibrated specifically for atmospheric wind farm flows by van der Laan (2014), and it showed much improvement over the standard k–ε model. It is essentially equivalent to the first-order model of the non-linear eddy viscosity model (NLEVM) of Apsley and Leschziner (1998).

Linear EVMs based on the Boussinesq hypothesis (e.g., mixing length, k–ε and k–ω models) are the de facto standard for RANS in many fields of research, including wind energy applications (e.g., Bleeg et al., 2018; Hornshøj-Møller et al., 2021; Heinz et al., 2021; Dicholkar et al., 2022; Letizia and Iungo, 2022). They are simple to implement and numerically robust but have a drawback in terms of physical correctness – the example of wake recovery has already been mentioned earlier, and one can also quickly derive that the Reynolds stresses become unrealizable for large normalized strain rates in general; see Fig. 2. Full Reynolds stress modeling (RSM), where a transport equation for each of the six unique Reynolds stress components is solved, was conducted by Launder et al. (1975), and this approach leads to more physical predictions because one avoids the use of the quite limiting Boussinesq hypothesis, and hence the production of turbulence needs no further modeling, which is a major improvement. The drawback of this method is that it is more computationally expensive and, perhaps more importantly, less numerically stable. Rodi (1976) deployed the weak equilibrium approximation (WEA) to the RSM equations, which transforms the set of differential equations into a set of algebraic equations while still retaining most of the physical behavior of the RSM. This approach is known as the algebraic Reynolds stress model (ARSM), and it can be formulated as a tensor equation for the anisotropy tensor, but unfortunately it is a non-linear implicit equation with multiple solutions, which is also prone to numerical stability issues. Pope (1975) simplified the ARSM by using the Cayley–Hamilton theorem and obtained the first explicit algebraic Reynolds stress model (EARSM) for the case of two-dimensional flow, which as the name suggests is an explicit algebraic relation between the Reynolds stress tensor (or equivalently the anisotropy tensor) and the normalized strain rate and rotation rate tensors. Different three-dimensional generalizations were subsequently made by Taulbee (1992), Gatski and Speziale (1993), and Wallin and Johansson (2000), which differ in the constants and more importantly in the way they treat the non-linearity of the ARSM. The Wallin and Johansson (2000) EARSM (hereafter WJ-EARSM) is chosen for this work because it is based on the concept of self-consistency (meaning that the explicit solution of a_ij satisfies the non-linear ARSM tensor equation exactly in two-dimensional mean flows and approximately in three-dimensional mean flows) and is only slightly more computationally expensive compared to a standard k–ε model. The self-consistent solution was independently formulated by Girimaji (1996), Ying and Canuto (1996), and Johansson and Wallin (1996) and facilitates a consistent solution in non-equilibrium conditions (where TKE production does not balance dissipation of TKE).

It appears that only Gómez-Elvira et al. (2005) and van der Laan (2014) have attempted to use EARSM for RANS simulations of wind turbine wakes. Many other research areas have however used the WJ-EARSM successfully (e.g., airfoil flow, Franke et al., 2005; vortex generators, Jirásek, 2005; stirred tanks, Feng et al., 2012; Kaplan turbines, Javadi and Nilsson, 2017; and high-speed trains, Munoz-Paniagua et al., 2017), so it seems that an untapped potential exists in turbulence modeling for wind energy applications. Gómez-Elvira et al. (2005) deployed the Taulbee (1992) model in a parabolic RANS setup, which is fast to execute but lacks physical features such as the upstream induction zone; van der Laan (2014) tested the Taulbee (1992), Gatski and Speziale (1993), and Apsley and Leschziner (1998) models in the elliptic RANS solver EllipSys3D (Sørensen, 1995) but found that they were more numerically unstable compared to the standard k–ε model.

The EARSM framework of Wallin and Johansson (2000) has been further exploited by including the strong coupling with density stratification present in the ABL (Lazeroms et al., 2013), capturing the effect of both stable and convective ABLs (Lazeroms et al., 2016; Želi et al., 2020, 2021). In this paper, we restrict the study to the neutral atmospheric surface layer (ASL; i.e., the lower part of the ABL), leaving stratified conditions for upcoming studies.

Section 2 describes the turbulence model formulations. In Sect. 3, we verify our implementation of the WJ-EARSM using several canonical flow cases, and finally in Sects. 4 and 5, it is applied to simulations of wind turbine wakes in the neutral ASL.

2 Turbulence modeling

The turbulence models utilized in this paper assume incompressible, non-stratified flow; no system rotation (no Coriolis or centrifugal contributions); isotropic dissipation of TKE; and high Reynolds number flow.

2.1 The standard k–ε model (Launder and Sharma, 1974)

The Boussinesq hypothesis is used to obtain the Reynolds stresses

$$\begin{array}{}\text{(1)}& \stackrel{\mathrm{‾}}{u_i^{\prime }{u}_j^{\prime }}=-{\mathit{\nu }}_t\left({\displaystyle \frac{\partial U_i}{\partial x_j}}+{\displaystyle \frac{\partial U_j}{\partial x_i}}\right)+{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}k{\mathit{\delta }}_{ij},\end{array}$$

which are needed to close the momentum equations (the equations for the mean velocity vector, U_i). The eddy viscosity in a linear EVM is defined as ${\mathit{\nu }}_t=C_{\mathit{\mu }}\frac{k^{\mathrm{2}}}{\mathit{\epsilon }}$, and transport equations are used to obtain TKE (k) and dissipation of TKE (ε):

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle \frac{\partial k}{\partial t}}+U_j{\displaystyle \frac{\partial k}{\partial x_j}}=\underset{\mathcal{P}}{\underbrace{-\stackrel{\mathrm{‾}}{u_i^{\prime }{u}_j^{\prime }}{\displaystyle \frac{\partial U_i}{\partial x_j}}}}-\mathit{\epsilon }+\underset{{\mathcal{D}}^{\left(k\right)}}{\underbrace{{\displaystyle \frac{\partial }{\partial x_j}}\left({\displaystyle \frac{{\mathit{\nu }}_t}{{\mathit{\sigma }}_k}}{\displaystyle \frac{\partial k}{\partial x_j}}\right)}},\text{(3)}& {\displaystyle }& {\displaystyle \frac{\partial \mathit{\epsilon }}{\partial t}}+U_j{\displaystyle \frac{\partial \mathit{\epsilon }}{\partial x_j}}=\left(C_{\mathit{\epsilon }\mathrm{1}}\mathcal{P}-C_{\mathit{\epsilon }\mathrm{2}}\mathit{\epsilon }\right){\displaystyle \frac{\mathit{\epsilon }}{k}}+\underset{{\mathcal{D}}^{\left(\mathit{\epsilon }\right)}}{\underbrace{{\displaystyle \frac{\partial }{\partial x_j}}\left({\displaystyle \frac{{\mathit{\nu }}_t}{{\mathit{\sigma }}_{\mathit{\epsilon }}}}{\displaystyle \frac{\partial \mathit{\epsilon }}{\partial x_j}}\right)}}.\end{array}$$

This paper focuses on the k–ε model because it is traditionally used for atmospheric flows (e.g., Crespo et al., 1985; Richards and Hoxey, 1993; Sørensen, 1995), but the widely used k–ω model could as well have been used. They can both be categorized as linear EVMs because they use the Boussinesq hypothesis to obtain the Reynolds stresses. Several empirical constants are also present for both models, though they are related to each other (Sogachev and Kelly, 2012); for the k–ε model there are C_μ, σ_k, C_ε1, C_ε2, and σ_ε (see review by Weaver and Mišković, 2021, for the most popular sets of constants used in the past). We shall use different sets of constants throughout the paper and remark the choice at each usage of the k–ε model.

To simplify the Boussinesq hypothesis, Eq. (1), we can use the anisotropy tensor also mentioned in the introduction:

$$\begin{array}{}\text{(4)}& \mathbf{a}=a_{ij}\equiv {\displaystyle \frac{\stackrel{\mathrm{‾}}{u_i^{\prime }{u}_j^{\prime }}}{k}}-{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}{\mathit{\delta }}_{ij}.\end{array}$$

It is dimensionless, symmetric, and traceless, and Eq. (1) can then be re-written as

$$\begin{array}{}\text{(5)}& a_{ij}=-\mathrm{2}{C}_{\mathit{\mu }}{S}_{ij}\left(\mathrm{linear}\mathrm{EVM}\right),\end{array}$$

where $\mathbf{S}=S_{ij}\equiv \frac{\mathrm{1}}{\mathrm{2}}\frac{k}{\mathit{\epsilon }}\left(\frac{\partial U_i}{\partial x_j}+\frac{\partial U_j}{\partial x_i}\right)$ is the normalized strain rate tensor. The k–ε model is independent of the normalized rotation rate tensor, $\mathbf{\Omega }={\mathrm{\Omega }}_{ij}\equiv \frac{\mathrm{1}}{\mathrm{2}}\frac{k}{\mathit{\epsilon }}\left(\frac{\partial U_i}{\partial x_j}-\frac{\partial U_j}{\partial x_i}\right)$; hence it will be unable to predict turbulence effects associated with rotation or curvature, e.g., damping/enhancement of turbulence in rotating homogeneous shear flow (Wallin and Johansson, 2002). The more advanced EARSMs also need to solve the k and ε transport equations because they depend on S_ij and Ω_ij and hence depend on the turbulence timescale $k/\mathit{\epsilon }$. The key difference between these more advanced models and the standard k–ε closure is that the Boussinesq hypothesis is replaced by a more general constitutive relation.

Finally, it can be noted that the time derivative is retained in the transport equations, Eqs. (2)–(3), to allow for unsteady RANS (URANS) simulations, e.g., homogeneous shear flow. URANS is generally only advisable when there is a clear separation between the timescale of the turbulence and the mean unsteadiness (see, e.g., Wallin, 2000) and should therefore be used with care.

2.2 The k–ε–f_P model

The k–ε–f_P model by van der Laan (2014) is equivalent to the first-order model of Apsley and Leschziner (1998), except for a re-tuning of the model coefficients. To summarize,

$$\begin{array}{}\text{(6)}& a_{ij}=-\mathrm{2}\underset{C_{\mathit{\mu }}^{\mathrm{eff}}}{\underbrace{C_{\mathit{\mu }}{f}_P}}{S}_{ij}\left(k\text{–}\mathit{\epsilon }\text{–}{f}_P\text{model}\right).\end{array}$$

Compared with the Boussinesq hypothesis, Eq. (5), we see that the only difference is that the k–ε–f_P model uses a variable or “effective” $C_{\mathit{\mu }}^{\mathrm{eff}}$, which is flow-dependent and thereby makes a_ij non-linear in the velocity gradient tensor; however, the k–ε–f_P model is still referred to as a “linear EVM” (e.g., van der Laan, 2014) because the direction of a_ij is aligned with S_ij and not any other higher-order tensors. The f_P function is

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& f_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}}}}},\\ & \mathit{\sigma }\equiv {\displaystyle \frac{k}{\mathit{\epsilon }}}\sqrt{{\left({\displaystyle \frac{\partial U_i}{\partial x_j}}\right)}^{\mathrm{2}}}=\sqrt{I{I}_S-I{I}_{\mathrm{\Omega }}},f_{\mathrm{0}}=\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{C_R-\mathrm{1}}}.\end{array}\end{array}$$

In the above equation σ is the “shear parameter”, which can be re-written using velocity gradient tensor invariants (see Eq. 14), while $\stackrel{\mathrm{̃}}{\mathit{\sigma }}$ is the shear parameter of the freestream calibration flow. Instead of using DNS channel flow data for the calibration, as was done by Apsley and Leschziner (1998), van der Laan (2014) chose to calibrate with a neutral ASL (logarithmic wind profile), which simply gives $\stackrel{\mathrm{̃}}{\mathit{\sigma }}=C_{\mathit{\mu }}^{-\mathrm{1}/\mathrm{2}}$. Finally, van der Laan (2014) took the Rotta constant as a free parameter and tuned it to C_R=4.5 using LES data of wake velocity and TI (pressure-strain data were not used) extracted from eight different wind turbine wake cases.

In the freestream of a neutral ASL, f_P=1, so it reduces to the standard k–ε model, while f_P<1 in regions with rapid strain compared to the turbulence timescale (i.e., large normalized velocity gradients), hence attenuating mixing in the wake shear layers and improving predictions of wind turbine wakes as found by van der Laan (2014). The f_P correction reduces mainly the turbulence length scale (but also the turbulence velocity scale) in the near wake and can therefore be interpreted as a local turbulence-length-scale limiter, as discussed by van der Laan and Andersen (2018). In fact, the damping introduced through f_P will preserve realizability for rapid shear (large normalized velocity gradient); i.e., in the rapid distortion theory (RDT) limit it will predict $a_{\mathrm{13}}>-\mathrm{1}$, whereas the standard k–ε model would predict $a_{\mathrm{13}}<-\mathrm{1}$ (unrealizable turbulence); see Fig. 2.

Table 1 summarizes the model constants of the k–ε–f_P model, where one can notice that C_μ=0.03 is used instead of the established value of C_μ=0.09 (e.g., Launder and Spalding, 1974). Several measurements of atmospheric flow in flat terrain (e.g., Panofsky and Dutton, 1984) point to a lower C_μ, and Bottema (1997) argues that this is due to “inactive” low-frequency atmospheric turbulence; see also discussion by Richards and Norris (2011). As a compromise $C_{\mathit{\mu }}=\sqrt{\mathrm{0.09}\cdot \mathrm{0.03}}\approx \mathrm{0.052}$ is used in WAsP-CFD (Bechmann, 2016).

Table 1 Model constants for the k–ε–f_P model as recommended by Sørensen (1995) and van der Laan (2014). Referred to later as “ABL coefs”.

2.3 Wallin and Johansson (2000) EARSM

The EARS model of Wallin and Johansson (2000) is derived from the ARSM of Rodi (1976) and therefore inherits the constants, c₁ and c₂, which are the Rotta coefficient and rapid pressure-strain coefficient of the Launder et al. (1975) model, respectively. The particular choice $c_{\mathrm{2}}=\mathrm{5}/\mathrm{9}$ reduces the model expressions significantly and is adopted in this study. This choice is also supported by the DNSs of Shabbir and Shih (1993). Moreover, we only consider the incompressible, high-Reynolds-number version (without near-wall corrections) due to the high Reynolds number and the use of rough-wall function boundary conditions at the ground in the considered flow cases of this paper. Additionally, we make comparisons between the 2D and 3D models. It is important to stress that the 2D model is fully general and invariant and can be used for simulation of three-dimensional mean flows as noted by Hellsten and Wallin (2009); this has for example also been commonly done with the EARSM of Gatski and Speziale (1993). Another important point to mention again is that all EARSMs are based on combinations of S_ij and Ω_ij (which depend on the turbulence timescale, $\mathit{\tau }=k/\mathit{\epsilon }$), so one still needs to solve either the k–ε model (used in this paper), k–ω model (used by Wallin and Johansson, 2000), or some other combination for obtaining the turbulence timescale.

2.3.1 2D WJ-EARSM

Using the complete two-dimensional tensor basis for the anisotropy tensor (Pope, 1975), an exact self-consistent 2D EARSM was found independently by Johansson and Wallin (1996), Girimaji (1996), and Ying and Canuto (1996). It was more thoroughly elaborated and tested by Wallin and Johansson (2000), where derivation details also can be found. Without loss of generality, and for the reason of numerical implementation, the anisotropy is split into linear and “extra” terms:

$$\begin{array}{}\text{(8)}& a_{ij}=-\mathrm{2}{C}_{\mathit{\mu }}^{\mathrm{eff}}{S}_{ij}+a_{ij}^{\left(\mathrm{ex}\right)},\end{array}$$

where

$$\begin{array}{}\text{(9)}& C_{\mathit{\mu }}^{\mathrm{eff}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\beta }}_{\mathrm{1}},a_{ij}^{\left(\mathrm{ex}\right)}={\mathit{\beta }}_{\mathrm{4}}\left(S_{ik}{\mathrm{\Omega }}_{kj}-{\mathrm{\Omega }}_{ik}{S}_{kj}\right).\end{array}$$

The tensor coefficients are

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }{\mathit{\beta }}_{\mathrm{1}}=-{\displaystyle \frac{\mathrm{6}}{\mathrm{5}}}{\displaystyle \frac{N}{N^{\mathrm{2}}-\mathrm{2}I{I}_{\mathrm{\Omega }}}},\text{(11)}& {\displaystyle }& {\displaystyle }{\mathit{\beta }}_{\mathrm{4}}=-{\displaystyle \frac{\mathrm{6}}{\mathrm{5}}}{\displaystyle \frac{\mathrm{1}}{N^{\mathrm{2}}-\mathrm{2}I{I}_{\mathrm{\Omega }}}},\end{array}$$

while N is the real and positive root of a cubic polynomial related to the non-linearity of the EARSM:

$$\begin{array}{}\text{(12)}& {\displaystyle }\begin{array}{rl}& N=\\ & \left\{\begin{array}{ll}\frac{c_{\mathrm{1}}^{\prime }}{\mathrm{3}}+{\left(P_{\mathrm{1}}+\sqrt{P_{\mathrm{2}}}\right)}^{\mathrm{1}/\mathrm{3}}+\mathrm{sign}\left(P_{\mathrm{1}}-\sqrt{P_{\mathrm{2}}}\right)|P_{\mathrm{1}}-\sqrt{P_{\mathrm{2}}}{|}^{\mathrm{1}/\mathrm{3}},& P_{\mathrm{2}}\ge \mathrm{0}\\ \frac{c_{\mathrm{1}}^{\prime }}{\mathrm{3}}+\mathrm{2}{\left(P_{\mathrm{1}}^{\mathrm{2}}-P_{\mathrm{2}}\right)}^{\mathrm{1}/\mathrm{6}}+\cos \left(\frac{\mathrm{1}}{\mathrm{3}}\left(\frac{P_{\mathrm{1}}}{\sqrt{P_{\mathrm{1}}^{\mathrm{2}}-P_{\mathrm{2}}}}\right)\right),& P_{\mathrm{2}}<\mathrm{0}\end{array}\right.\end{array}\text{(13)}& {\displaystyle }\begin{array}{rl}& P_{\mathrm{1}}=\left({\displaystyle \frac{\mathrm{1}}{\mathrm{27}}}c{{}^{\prime }}_{\mathrm{1}}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{9}}{\mathrm{20}}}I{I}_S-{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}I{I}_{\mathrm{\Omega }}\right)c_{\mathrm{1}}^{\prime },\\ & P_{\mathrm{2}}=P_{\mathrm{1}}^{\mathrm{2}}-{\left({\displaystyle \frac{\mathrm{1}}{\mathrm{9}}}c{{}^{\prime }}_{\mathrm{1}}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{9}}{\mathrm{10}}}I{I}_S+{\displaystyle \frac{\mathrm{2}}{\mathrm{3}}}I{I}_{\mathrm{\Omega }}\right)}^{\mathrm{3}},\\ & c_{\mathrm{1}}^{\prime }={\displaystyle \frac{\mathrm{9}}{\mathrm{4}}}\left(c_{\mathrm{1}}-\mathrm{1}\right).\end{array}\end{array}$$

The model only depends on the first two velocity gradient tensor invariants:

$$\begin{array}{}\text{(14)}& I{I}_S\equiv S_{ij}{S}_{ji},I{I}_{\mathrm{\Omega }}\equiv {\mathrm{\Omega }}_{ij}{\mathrm{\Omega }}_{ji}.\end{array}$$

The solution procedure is more thoroughly described in Sect. 2.4, but one can just notice that given S_ij and Ω_ij, there is a closed and explicit solution of a_ij. Since N, Eq. (12), is an exact solution of the underlying cubic polynomial problem, a_ij is an exact solution of the ARSM; hence it is a self-consistent EARSM, which ensures good model predictions in non-equilibrium conditions for $\mathcal{P}/\mathit{\epsilon }>\mathrm{1}$. We suspect that the treatment of the EARSM non-linearity by Taulbee (1992), Gatski and Speziale (1993), and Apsley and Leschziner (1998) is one of the main causes for the numerical instability when these EARSMs/NLEVMs are employed to model wind turbine wakes (van der Laan, 2014).

C_μ is not an input parameter but an output parameter as a result of applying the WJ-EARSM; see Eq. (9). However, one still needs it for the rough-wall boundary condition (BC) used in our code and for the diffusion terms in the k and ε transport equations in our implementation; here we shall use C_μ=0.087, which is the equilibrium value (the value obtained by fixing $\mathcal{P}/\mathit{\epsilon }=\mathrm{1}$ and considering the log-layer relations) with the current set of constants in Table 2; see Sect. 4.1. Also note that κ=0.38 is not an input parameter but is calculated from the other model constants in order to satisfy the log-layer balance; see Eq. (23).

Table 2 Model constants for the WJ-EARSM as recommended by Želi et al. (2020). Referred to later as “Zeli coefs”.

2.3.2 3D WJ-EARSM

The 3D model is derived in an analogous way to the 2D model (Wallin and Johansson, 2000) but with the complete three-dimensional tensor representation of a_ij, which gives

$$\begin{array}{}\text{(15)}& a_{ij}=-\mathrm{2}{C}_{\mathit{\mu }}^{\mathrm{eff}}{S}_{ij}+a_{ij}^{\left(\mathrm{ex}\right)},\end{array}$$

with

$$\begin{array}{}\text{(16)}& \begin{array}{rl}& C_{\mathit{\mu }}^{\mathrm{eff}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\mathit{\beta }}_{\mathrm{1}}+I{I}_{\mathrm{\Omega }}{\mathit{\beta }}_{\mathrm{6}}\right),\\ & a_{ij}^{\left(\mathrm{ex}\right)}={\mathit{\beta }}_{\mathrm{3}}{T}_{ij}^{\left(\mathrm{3}\right)}+{\mathit{\beta }}_{\mathrm{4}}{T}_{ij}^{\left(\mathrm{4}\right)}+{\mathit{\beta }}_{\mathrm{6}}\left(T_{ij}^{\left(\mathrm{6}\right)}-I{I}_{\mathrm{\Omega }}{S}_{ij}\right)+{\mathit{\beta }}_{\mathrm{9}}{T}_{ij}^{\left(\mathrm{9}\right)}.\end{array}\end{array}$$

The tensor coefficients and basis tensors are

${\mathit{\beta }}_{\mathrm{1}}=-\frac{N\left(\mathrm{2}{N}^{\mathrm{2}}-\mathrm{7}I{I}_{\mathrm{\Omega }}\right)}{Q}$			${\mathbf{T}}^{\left(\mathrm{3}\right)}={\mathbf{\Omega }}^{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{3}}I{I}_{\mathrm{\Omega }}\mathbf{I}$
${\mathit{\beta }}_{\mathrm{3}}=-\frac{\mathrm{12}{N}^{-\mathrm{1}}IV}{Q}$			${\mathbf{T}}^{\left(\mathrm{4}\right)}=\mathbf{S}\mathbf{\Omega }-\mathbf{\Omega }\mathbf{S}$
${\mathit{\beta }}_{\mathrm{4}}=-\frac{\mathrm{2}\left(N^{\mathrm{2}}-\mathrm{2}I{I}_{\mathrm{\Omega }}\right)}{Q}$			${\mathbf{T}}^{\left(\mathrm{6}\right)}=\mathbf{S}{\mathbf{\Omega }}^{\mathrm{2}}+{\mathbf{\Omega }}^{\mathrm{2}}\mathbf{S}-\frac{\mathrm{2}}{\mathrm{3}}IV\mathbf{I}$
${\mathit{\beta }}_{\mathrm{6}}=-\frac{\mathrm{6}N}{Q}$			${\mathbf{T}}^{\left(\mathrm{9}\right)}=\mathbf{\Omega }\mathbf{S}{\mathbf{\Omega }}^{\mathrm{2}}-{\mathbf{\Omega }}^{\mathrm{2}}\mathbf{S}\mathbf{\Omega }$,
${\mathit{\beta }}_{\mathrm{9}}=\frac{\mathrm{6}}{Q}$

where

$$\begin{array}{}\text{(17)}& {\displaystyle }& {\displaystyle }Q={\displaystyle \frac{\mathrm{5}}{\mathrm{6}}}\left(N^{\mathrm{2}}-\mathrm{2}I{I}_{\mathrm{\Omega }}\right)\left(\mathrm{2}{N}^{\mathrm{2}}-I{I}_{\mathrm{\Omega }}\right),\text{(18)}& {\displaystyle }& {\displaystyle }IV=S_{ij}{\mathrm{\Omega }}_{jk}{\mathrm{\Omega }}_{ki}.\end{array}$$

Unfortunately, there does not exist an analytical solution for N in the 3D model; therefore the cubic N solution from 2D, Eq. (12), is used; hence the 3D model is not exactly self-consistent, but the cubic N solution is however still a quite good approximation in many cases. The same model constants are also used; see Table 2.

One can show that the 3D model reduces to the 2D model in two-dimensional mean flows, where ${\mathbf{T}}^{\left(\mathrm{6}\right)}-I{I}_{\mathrm{\Omega }}\mathbf{S}=\mathrm{0}$, IV=0, and ${\mathbf{T}}^{\left(\mathrm{9}\right)}=-\frac{\mathrm{1}}{\mathrm{2}}I{I}_{\mathrm{\Omega }}{\mathbf{T}}^{\left(\mathrm{4}\right)}$ are valid.

2.4 Implementation details

The flow cases are simulated with EllipSys3D, which is a finite-volume CFD solver developed and described in detail by Michelsen (1992) and Sørensen (1995). The solver already has implementations of the standard k–ε model and the k–ε–f_P model (van der Laan, 2014), so the following focuses on the implementation of the WJ-EARSM. As emphasized in the previous sections, both the 2D and 3D models can be written in the same form; see Eqs. (8) and (15) (although with different expressions for $C_{\mathit{\mu }}^{\mathrm{eff}}$ and $a_{ij}^{\left(\mathrm{ex}\right)}$). The splitting of the anisotropy tensor into a linear and an extra part makes the implementation relatively straightforward in codes that already have a k–ε model implemented, as also noted in Appendix A of Wallin and Johansson (2000); for the momentum equations, simply use ${\mathit{\nu }}_t^{\mathrm{eff}}=C_{\mathit{\mu }}^{\mathrm{eff}}\frac{k^{\mathrm{2}}}{\mathit{\epsilon }}$ instead of ${\mathit{\nu }}_t=C_{\mathit{\mu }}\frac{k^{\mathrm{2}}}{\mathit{\epsilon }}$, and add $-\frac{\partial a_{ij}^{\left(\mathrm{ex}\right)}k}{\partial x_j}$ as a source term:

$$\begin{array}{}\text{(19)}& \begin{array}{rl}-{\displaystyle \frac{\partial \stackrel{\mathrm{‾}}{u_i^{\prime }{u}_j^{\prime }}}{\partial x_j}}& =-{\displaystyle \frac{\partial \left(a_{ij}k+\frac{\mathrm{2}}{\mathrm{3}}k{\mathit{\delta }}_{ij}\right)}{\partial x_j}}\\ & =\underset{\mathrm{Treat}\mathrm{implicit}}{\underbrace{{\displaystyle \frac{\partial \mathrm{2}{C}_{\mathit{\mu }}^{\mathrm{eff}}k{S}_{ij}}{\partial x_j}}}}-\underset{\mathrm{Treat}\mathrm{explicit}}{\underbrace{{\displaystyle \frac{\partial a_{ij}^{\left(\mathrm{ex}\right)}k}{\partial x_j}}}}-\underset{\mathrm{Absorb}\mathrm{into}\mathrm{pressure}}{\underbrace{{\displaystyle \frac{\partial \frac{\mathrm{2}}{\mathrm{3}}k{\mathit{\delta }}_{ij}}{\partial x_j}}}}.\end{array}\end{array}$$

For numerical stability, it is recommended to include the first term in the system matrix and the second term in the source vector (i.e., treat the terms implicitly and explicitly, respectively). The third term is isotropic and can be absorbed into a modified pressure.

In the k and ε transport equations, we use the standard ν_t in the diffusion terms (𝒟^(k) and 𝒟^(ε) in Eqs. 2–3) since these are calibrated using the standard model. This practice is also used by Apsley and Leschziner (1998), Myllerup (2000), and Menter et al. (2009). In the original paper of Wallin and Johansson (2000), it was proposed to use the Daly–Harlow diffusion model or an eddy-diffusion model with the effective ${\mathit{\nu }}_t^{\mathrm{eff}}$, which was later abandoned for the standard eddy-diffusivity model (see, e.g., Menter et al., 2012). The TKE production is calculated consistently as $\mathcal{P}=-\mathit{\epsilon }{a}_{ik}{S}_{ki}$ with the full anisotropy tensor; this expression is in fact valid for all turbulence models (see for example Wallin and Johansson, 2000).

A segregated solver (i.e., solving the U equation, then the V equation, then the W equation, etc.) is used in EllipSys, and the same shear production 𝒫 is used in both turbulence transport equations; an overview of the procedure is sketched in Fig. 1. The boundary conditions for each variable depend on the type of boundary (rough wall, inlet, outlet, cyclic, or symmetric), and the implementation of these in EllipSys3D is described by Sørensen (1995) and Sørensen et al. (2007).

Figure 1 Segregated solver procedure for WJ-EARSM in EllipSys3D. Note that P is the mean pressure, while 𝒫 is the TKE shear production.

The 1D version of EllipSys3D, EllipSys1D (van der Laan and Sørensen, 2017), is used for the 1D verification cases in Sect. 3. Its code structure is very similar to EllipSys3D's, but it has no W momentum equation (W=0) and no pressure correction equation ($\partial /\partial x=\partial /\partial y=\mathrm{0}$, so continuity is ensured already), and it is written using 1D simplifications (mean variations only in time and z direction), which makes it possible to run 1D simulations, or so-called “single-column models”, on a regular laptop in a few seconds.

For both EllipSys1D and EllipSys3D, the procedure of calling the WJ-EARSM is as such:

1. Use the most recent solution of momentum and turbulence transport equations to calculate the normalized strain rate and rotation rate tensors, S_ij and Ω_ij.

2. Calculate tensors and invariants.

3. Calculate N.

4. Calculate coefficients, β₁, …, β₁₀.

5. Calculate anisotropy tensor, a_ij.

6. Calculate TKE shear production, $\mathcal{P}=-\mathit{\epsilon }{a}_{ik}{S}_{ki}$.

7. Calculate $C_{\mathit{\mu }}^{\mathrm{eff}}$ and $a_{ij}^{\left(\mathrm{ex}\right)}$.

3 Verification cases

As is clear from the previous section, the expressions in the WJ-EARSM are considerably longer compared to the ones of the k–ε and k–ε–f_P models. Three canonical flows (homogeneous shear flow, half-channel flow, and square duct flow) are therefore used as verification cases to gain confidence in the numerical implementation (verification is the process of ensuring correct implementation, whereas validation is the process of assessing a model's accuracy; see Réthoré et al., 2014); full three-dimensional wind turbine wake simulations are first considered in Sect. 4. The first two cases can be simulated in a 1D setup, which make them ideal for initial testing, while the last case needs to be simulated in 3D due to the phenomenon of secondary motions. For the latter, a quasi-2D setup – homogeneous in the streamwise direction – would be possible, but a fully three-dimensional setup is chosen for verification of the 3D implementation. All cases are compared to either analytic expressions or DNS data to verify correct behavior of the implementation.

Even before running the verification cases, we consider the general class of “simple shear flows” (a.k.a. 1D parallel flows), where the normalized strain rate and rotation rate tensors are

$$\begin{array}{}\text{(20)}& \begin{array}{rl}& {\mathbf{S}}_{\mathrm{simple}}=\left(\begin{array}{ccc}\mathrm{0}& \mathrm{0}& S_{\mathrm{13}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ S_{\mathrm{13}}& \mathrm{0}& \mathrm{0}\end{array}\right),\\ & {\mathbf{\Omega }}_{\mathrm{simple}}=\left(\begin{array}{ccc}\mathrm{0}& \mathrm{0}& S_{\mathrm{13}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -S_{\mathrm{13}}& \mathrm{0}& \mathrm{0}\end{array}\right).\end{array}\end{array}$$

Prescribing S₁₃ enables us to evaluate the turbulence models analytically and obtain the anisotropy tensor, a_ij. The a₁₃ component is shown in Fig. 2, and from the Cauchy–Schwarz inequality one can show generally that $-\mathrm{1}\le a_{\mathrm{13}}\le \mathrm{1}$; if a turbulence model violates this, then it predicts “unrealizable” turbulence, meaning unphysical turbulence. For large S₁₃, the k–ε model leads to unrealizable turbulence, while both the k–ε–f_P model and WJ-EARSM (the 2D and 3D WJ-EARSMs are identical in 1D) are realizable even for large S₁₃. This is very desirable because large normalized velocity gradients are typically encountered in regions of non-equilibrium turbulence, e.g., in the vicinity of the rotor and in the wake shear layers. Indeed Réthoré (2009) noted large regions of unrealizable turbulence near the wind turbine when using the standard k–ε model. Also note that the k–ε–f_P model (using C_μ=0.03) predicts a similar a₁₃ as the WJ-EARSM solution for large strains, while the main difference between the models is found for small strains.

Figure 2 Analytical off-diagonal anisotropy in simple shear flow.

3.1 Homogeneous shear flow

Homogeneous shear flow (see review by Pope, 2000, pp. 154–157) has a simple setup but can be challenging to simulate and is conceptually a strange case. The momentum equations are not solved, and instead a constant velocity gradient is artificially fixed at all times; see Fig. 3. This is, indeed, rather unphysical since the velocity gradient should gradually decrease as turbulence is created until the velocity gradient eventually becomes zero, and turbulence dies out. Nevertheless, homogeneous shear flow constitutes an interesting test case because it only involves the solution of the turbulence transport equations and turbulence closure; hence it is a “pure” test of the turbulence model. Moreover, free-shear layers can locally be approximated with homogeneous shear. Only EllipSys1D is used for this case.

Figure 3 Mesh and the prescribed velocity profile in EllipSys1D for the homogeneous-shear-flow case.

In principle, homogeneous shear flow is an unsteady 0D case. The turbulence evolves identically at all positions in space; hence there is actually no need for a spatial discretization as shown in Fig. 3 (one could simply integrate the k and ε equations forward in time, given an initial turbulent state, as done by for example Taulbee, 1992), but this extra complexity is chosen for the present case since the goal is to verify the turbulence model implemented in a finite-volume CFD code. In practice this means that BCs need to be set at the top and bottom of the domain: a fixed velocity, U, is used as illustrated in Fig. 3, while symmetry BCs are used for k and ε. The implicit Euler scheme is used for time integration, while the second-order central scheme is used for diffusive terms.

The simulation parameters used are as follows.

• Height of domain, L_z=100 m. Uniform spacing and 10 finite-volume cells are used (if ghost cells are counted in, there are 12 cells); hence grid spacing is Δz=10 m.

• Fixed velocity profile, U(z)=Sz, where S=0.1 s⁻¹. This gives $U(z=\mathrm{0})=\mathrm{0}$ and U(L_z)=10 m s⁻¹.

• Normalized time step, $\mathrm{\Delta }{t}^{*}\equiv \mathrm{\Delta }tS=\mathrm{0.1}$ (in physical time Δt=1 s). Ten subiterations per time step are used, and the total normalized simulation time is $t^{*}=\mathrm{80}$ (in physical time t=800 s).

• Initial turbulent state, $S\frac{k}{\mathit{\epsilon }}=\mathrm{3.4}$ at $t^{*}=\mathrm{0}$. This non-dimensional quantity determines the evolution of the non-dimensional metrics, i.e., anisotropy components, production-to-dissipation ratio, etc.; it is the same initial turbulent state as used in the homogeneous shear flow simulations by Bardina et al. (1983), Gatski and Speziale (1993), Girimaji (1996), and Wallin and Johansson (2002).

The RANS model equations here are independent of ρ and ν because no momentum equation is solved, and the diffusion terms in the k and ε equations (𝒟^(k) and 𝒟^(ε)) in Eqs. (2)–(3) will be zero.

Figure 4 Simulation of homogeneous shear flow with WJ-EARSM (full line) and the model's analytical asymptotic values (dashed lines) with different sets of constants. The simulation data are extracted at $z/L_z=\mathrm{0.45}$, but the evolution is identical at all grid points.

A nice feature of homogeneous shear flow for verification purposes is that it evolves to an asymptotic state and that an analytic solution exists for this state; e.g., the asymptotic production-to-dissipation ratio can be derived to be (see for example Gatski and Speziale, 1993)

$$\begin{array}{}\text{(21)}& {\displaystyle \frac{\mathcal{P}}{\mathit{\epsilon }}}={\displaystyle \frac{\mathrm{1}-C_{\mathit{\epsilon }\mathrm{2}}}{\mathrm{1}-C_{\mathit{\epsilon }\mathrm{1}}}}(\mathrm{asymptotic}\mathrm{limit}t\to \mathrm{\infty }).\end{array}$$

Table 3 Analytic, asymptotic homogeneous shear flow formulas. In homogeneous shear flow, $S_{\mathrm{13}}=S_{\mathrm{31}}={\mathrm{\Omega }}_{\mathrm{13}}=-{\mathrm{\Omega }}_{\mathrm{31}}$, while all other components are 0. When $\frac{\mathcal{P}}{\mathit{\epsilon }}$ is known, one can use the direct definition of $N\equiv c_{\mathrm{1}}^{\prime }+\frac{\mathrm{9}}{\mathrm{4}}\frac{\mathcal{P}}{\mathit{\epsilon }}$ instead of Eq. (12); see Wallin and Johansson (2000).

In 1D parallel flows, $\mathcal{P}/\mathit{\epsilon }=-\mathrm{2}{a}_{\mathrm{13}}{S}_{\mathrm{13}}$, and this gives a second equation for $\mathcal{P}/\mathit{\epsilon }$. By inserting the expression for a₁₃ (the expression differs depending on the turbulence model), one can isolate and obtain a value for $S_{\mathrm{13}}=\frac{\mathrm{1}}{\mathrm{2}}\frac{k}{\mathit{\epsilon }}S$. This value can then be used to obtain S_ij and Ω_ij and hence the anisotropy tensor. The asymptotic formulas are summarized in Table 3, where one can note that the asymptotic WJ-EARSM formulas only depend on C_ε1 and C_ε2, while the k–ε formulas additionally depend on C_μ. Figure 4 shows the time evolution of the WJ-EARSM as well as the theoretical asymptotic values for three common sets of coefficients (C_ε1, C_ε2). All simulations indeed go to the analytical asymptote, which gives some confidence in the implementation of the WJ-EARSM. For validation purposes we also include some LES (Bardina et al., 1983) and experimental (Tavoularis and Corrsin, 1981) data in Fig. 4, which show that the Zeli constants perform better for all quantities except a₁₁. The choice of $c_{\mathrm{2}}=\mathrm{5}/\mathrm{9}$ makes a₂₂=0 for all simulations.

Figure 5 Half-channel EllipSys1D setup. The sketch to the right shows the convention of the coordinate system used for rough-wall simulations.

3.2 Half-channel flow

The second reference flow is the fully developed, steady-state half-channel flow, a.k.a. pressure-driven boundary layer (PDBL) flow, and can be solved in 1D using a grid as sketched in Fig. 5. The lower BC is a rough wall implemented as in Sørensen et al. (2007), while the upper BC uses symmetry; hence turbulence now depends on the z coordinate contrary to in the homogeneous-shear-flow case. The flow is driven by a constant streamwise pressure gradient force, which is an input parameter for the simulation (if one does not use such a forcing, then the rough wall will extract momentum from the flow, and the velocity eventually becomes zero throughout the domain).

The input parameters for the simulations are as follows:

• Height of domain, L_z=6000 m. The grid is stretched using the hyperbolic tangent function (Thompson et al., 1985), and 192 finite-volume cells are used (if ghost cells are counted in, there are 194 cells).

• Aerodynamic roughness height, z₀=0.03 m. This is used for the rough-wall BC in the first wall-adjacent cell, which has a height of Δz=0.10 m.

• Streamwise pressure gradient force per unit mass, $F_p=\mathrm{1.5}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻². One can show that the squared friction velocity then becomes $u_{*}^{\mathrm{2}}\equiv -{\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}}_s=F_p{L}_z=\mathrm{0.09}$ m s⁻¹.

The flow is independent of ν because a rough-wall BC is used and the flow is fully turbulent, hence ν≪ν_t.

The region in the lower part of the domain is known as the “log layer”, which is characterized by equilibrium turbulence, i.e., $\frac{\mathcal{P}}{\mathit{\epsilon }}\approx \mathrm{1}$. As in the homogeneous-shear-flow case, we can use this ratio to obtain analytical results for shear and anisotropy; see Table 4.

Table 4 Analytical values in the log layer of the half-channel flow.

Figure 6 Half-channel simulation results (full lines) and analytical log-layer solutions (dashed lines).

The flow profiles with the analytical log-layer values are plotted in Fig. 6. Both sets of model constants, Tables 1 and 2, are used for the k–ε model, which explains the different U profiles. A spike in the k profiles is seen near the wall, which is a well-known problem (see Blocken et al., 2007), but the value of k(z≈0) is nevertheless close to the equilibrium value $k_{\mathrm{eq}}=\frac{u_{*}^{\mathrm{2}}}{C_{\mathit{\mu }}^{\mathrm{1}/\mathrm{2}}}$ for both models (again reminding that C_μ differs between the two sets of model constants; compare Tables 1 and 2). The kinematic wall shear stress, $\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}(z=\mathrm{0})$, is close to $-u_{*}^{\mathrm{2}}$ for both models, confirming that the pressure gradient force is applied correctly, and the respective analytical log-layer solutions are also approximately obtained for both models at $z/L_z<\mathrm{0.3}$ in accordance with Pope (2000). One can notice that a feature of the WJ-EARSM is that the normal streamwise and vertical anisotropies are non-zero; this means $\stackrel{\mathrm{‾}}{u^{\prime }{u}^{\prime }}>\stackrel{\mathrm{‾}}{v^{\prime }{v}^{\prime }}=\frac{\mathrm{2}}{\mathrm{3}}k>\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }}$. This behavior is also seen in the neutral ASL (Panofsky and Dutton, 1984), which we return to in Sect. 4.

3.3 Square duct flow

The square duct geometry and boundary conditions are shown in Fig. 7. Similar to the half-channel flow, the flow is driven by a streamwise pressure gradient, but the difference is that the duct flow has walls on all four sides. Due to symmetry, we only simulate the lower left quadrant of the duct. As this paper focuses on high-Reynolds-number turbulence models, we choose to model the walls as rough walls instead of smooth walls, which are traditionally used in square duct flow simulations.

Figure 7 Square duct geometry and a cross-section showing the secondary corner flows. Sketch made with inspiration from Wu et al. (2016).

Although, the fully developed square duct flow might appear as a two-dimensional problem, it in fact features a full three-dimensional flow field, due to the secondary corner flows, also sketched in Fig. 7, which were first observed experimentally by Nikuradse (1930) and later with DNS by Gavrilakis (1992). The secondary motions are only on the order of $\frac{V}{U}\sim \frac{W}{U}\sim \mathcal{O}\left({\mathrm{10}}^{-\mathrm{2}}\right)$ but still have a notable effect on the bulk flow as they transport momentum from the center of the duct toward the corners. Perhaps the most interesting aspect of square duct flow from a turbulence modeling perspective is that linear EVMs are unable to predict the secondary corner flows because the secondary motions are caused by the normal anisotropy components, which are zero in linear EVMs for fully developed flow; see Eq. (5) and discussions by Menter et al. (2009) and Emory et al. (2013). However more sophisticated turbulence models such as EARSM and uncertainty quantification models (Emory et al., 2013) are able to predict this physical phenomenon.

Figure 8 Streamwise (a–c) and vertical (d–f) velocity contours in the lower left quadrant of the square duct. Normalized by the bulk velocity, $U_{\mathrm{bulk}}=\frac{\mathrm{1}}{A}\int U\mathrm{d}A$.

In this paper, we simulate a fully turbulent square duct flow (high Reynolds number) with rough-wall BCs; hence the flow is independent of ν similar to the half-channel flow. Currently, the DNS of Pirozzoli et al. (2018) is the most turbulent DNS available ($R{e}_{\mathit{\tau }}\equiv \frac{\langle u_{*}\rangle h}{\mathit{\nu }}=\mathrm{1055}$), so this is chosen for reference, although it uses smooth walls and might not correspond exactly to our high-Reynolds-number case. For this reason the RANS and DNS should not be compared directly, but the DNS will at least show typical characteristics and the order of magnitude to be expected. The parameters for the RANS simulation are given in Table 5, and a rectilinear $\mathrm{32}\times \mathrm{64}\times \mathrm{64}$ grid is used, which is stretched towards the walls to obtain a first wall-adjacent cell height on the order of the roughness length (no stretching used in streamwise direction).

Table 5 Parameters used for simulation of square duct flow in EllipSys3D. The brackets 〈〉 signify the average over the wall.

Table 6 Simulation parameters for the single-wake V80 case. Note that the LES uses neutral PDBL for the atmospheric model.

Figure 9 Streamwise velocity (a) and vertical velocity (b) on the corner bisector line. Results extracted at $x/L_x=\mathrm{0.5}$.

The streamwise and vertical velocities (U and V components) are shown in Fig. 8, which shows that the 2D WJ-EARSM is indeed capable of predicting secondary flows similar to the DNS of Pirozzoli et al. (2018). Both Figs. 8 and 9 clearly show that the prediction of the secondary flow is necessary to capture the correct shape of the U distribution. In contrast, the standard k–ε model predicts zero vertical velocity and therefore no secondary flow.

The 3D WJ-EARSM performed similarly to the 2D WJ-EARSM, although with a slightly weaker secondary flow, which can be seen in Fig. 9, where the velocity profiles are extracted on the corner bisector line (the diagonal line). Both WJ-EARSMs predict similar profiles to the DNS, although without the near-corner peak, which according to Pirozzoli et al. (2018) is caused by scale separation (this phenomena only occurs at higher Re_τ and is thus not visible in the earlier DNSs by Gavrilakis, 1992, and Huser and Biringen, 1993), and this effect is not captured by RANS.

Figure 10 Flow domain (a), finely resolved wake domain (b), x–y cut at hub height (c), and x–z cut at centerline (d). Every fourth cell is displayed in panels (c) and (d), and green arrows show velocity profile (not to scale).

This concludes the verification studies, where the WJ-EARSM has been seen to give expected results for three canonical flows. Furthermore, the last case of square duct flow clearly demonstrates that EARS models are able to predict physical phenomena that two-equation models based on the linear Boussinesq hypothesis cannot do.

4 Single-wind-turbine wake

This section concerns the application of the EARS model to a single-wind-turbine wake. The numerical CFD setup is similar to that used in many previous RANS studies (e.g., van der Laan, 2014; van der Laan et al., 2021; Baungaard et al., 2022): the RANS equations are solved using the SIMPLE method using a modified Rhie–Chow algorithm (Troldborg et al., 2015), while the convective terms are discretized with the QUICK scheme; see more details in Sørensen (1995). To model the wind turbine, an AD with uniform distribution of forces is used (the thrust forces are fixed, and no tangential forces are present), and the forces are transferred to the rectilinear flow domain with the intersectional method of Réthoré et al. (2014). The flow domain has a finely resolved “wake domain” in the center with uniform spacing of $D/\mathrm{8}$, and the grid is stretched outwards in all directions from this using the hyperbolic tangent method of Thompson et al. (1985); see Fig. 10. A grid convergence study was made to confirm that this grid resolution is also suitable for the EARS turbulence model; see Appendix A.

The case simulated is similar to the case used by Hornshøj-Møller et al. (2021), namely a single V80 turbine subject to neutral inflow; see Table 6. The authors of the aforementioned study have provided LES data to us, which are used as a reference in the following. It should be noted that their RANS simulations use ASL inflow (like we shall also use for our RANS simulations), while their LES is based on PDBL inflow. Although there will be differences between ASL and PDBL inflow, the latter is likely a good approximation of the former in the lower part of the domain, but its bias on wake simulations could be a subject for future studies.

4.1 Inflow

The neutral ASL inflow profile (Panofsky and Dutton, 1984) is prescribed at the inlet BC and top of the domain:

$$\begin{array}{}\text{(22)}& U\left(z\right)={\displaystyle \frac{u_{*}}{\mathit{\kappa }}}\ln \left({\displaystyle \frac{z}{z_{\mathrm{0}}}}\right),k\left(z\right)={\displaystyle \frac{u_{*}^{\mathrm{2}}}{\sqrt{C_{\mathit{\mu }}}}},\mathit{\epsilon }\left(z\right)={\displaystyle \frac{u_{*}^{\mathrm{3}}}{\mathit{\kappa }z}}.\end{array}$$

This type of inflow is routinely used in wind energy applications, and for a wake simulation it can be adopted to give a desired hub height velocity, U_ref, and hub height turbulence intensity (TI), I_ref, by adjusting u_* and z₀ (van der Laan et al., 2015b). One could alternatively adjust C_μ instead of z₀ to obtain a desired I_ref, but this was shown by van der Laan (2014) to give inconsistent results with the k–ε–f_P model and NLEVMs, so this practice is not recommended. Although the adjusted u_* and z₀ do not correspond to the physical values at the site, it is of higher priority to have the correct U_ref (and thereby correct thrust coefficient) and I_ref.

Figure 11 WJ-EARSM dependence on c₁ in the neutral ASL ($\mathcal{P}/\mathit{\epsilon }=\mathrm{1}$).

Figure 12 Inflow profiles for the single-wake V80 case.

In the freestream, the neutral ASL profiles, Eq. (22), should satisfy $Dk/Dt=\mathrm{0}$ and $D\mathit{\epsilon }/Dt=\mathrm{0}$ to be in balance and mitigate development of the inflow profiles in the streamwise direction. There will inevitably be a slight development (see Blocken et al., 2007), and for this reason a long domain is used to ensure fully developed profiles at the entrance of the wake domain. To satisfy the balance criteria, the turbulence constants should follow the relation (Richards and Hoxey, 1993)

$$\begin{array}{}\text{(23)}& C_{\mathit{\epsilon }\mathrm{1}}=C_{\mathit{\epsilon }\mathrm{2}}-{\displaystyle \frac{{\mathit{\kappa }}^{\mathrm{2}}}{\sqrt{C_{\mathit{\mu }}}{\mathit{\sigma }}_{\mathit{\epsilon }}}}.\end{array}$$

Indeed, both sets of constants in Tables 1 and 2 satisfy Eq. (23).

In the neutral ASL, we have $\mathcal{P}/\mathit{\epsilon }=\mathrm{1}$, similar to the log layer of the half-channel flow but at all heights. By utilizing the WJ-EARSM the equilibrium relations between the c₁ constant and various other variables, as shown in Fig. 11, can be derived analytically. For example, the standard value of c₁=1.8, which is used in this paper, gives C_μ=0.087 (the “equilibrium value” mentioned in Sect. 2.3.1). A closer agreement of the velocity standard deviation ratios, ${\mathit{\sigma }}_v/{\mathit{\sigma }}_u$ and ${\mathit{\sigma }}_w/{\mathit{\sigma }}_u$, is seen with the WJ-EARSM, when benchmarking against the ASL ratios of Panofsky and Dutton (1984); the k–ε and k–ε–f_P models have ${\mathit{\sigma }}_v/{\mathit{\sigma }}_u={\mathit{\sigma }}_w/{\mathit{\sigma }}_u=\mathrm{1}$, while the WJ-EARSM has ${\mathit{\sigma }}_v/{\mathit{\sigma }}_u=\mathrm{0.85}$ and ${\mathit{\sigma }}_w/{\mathit{\sigma }}_u=\mathrm{0.68}$ for c₁=1.8. In fact, the WJ-EARSM underpredicts the streamwise fluctuations due to the simplification introduced by setting $c_{\mathrm{2}}=\mathrm{5}/\mathrm{9}$ (this is also the case for half-channel flow; see Wallin and Johansson, 2000), resulting in the observed overprediction of the ${\mathit{\sigma }}_v/{\mathit{\sigma }}_u$ and ${\mathit{\sigma }}_w/{\mathit{\sigma }}_u$ ratios. Other c₁ values (three other choices than the standard value are marked with dots in Fig. 11) will enhance/decrease the ASL anisotropy (see details in Appendix B), but for the present simulations the standard model coefficients of Table 2 are used.

The inflow profiles of the LES and RANS simulations are shown in Fig. 12. Since the k–ε and k–ε–f_P models are identical in the freestream, which is also true for the 2D and 3D WJ-EARSMs, only one of the other is shown in the figure. The LES is driven with a streamwise pressure gradient and thus differs in the stress profiles compared to the pure ASL profiles used in our RANS, but U_ref and I_ref do match. As in the half-channel case, the velocity shear differs due to the different κ and C_μ used in the k–ε–f_P model and WJ-EARSM, and the stress profiles are also different between the two because of the anisotropic nature of WJ-EARSM. In terms of turbulence anisotropy, quantified as the distribution of TKE between the normal stress components in Fig. 12, the WJ-EARSM is clearly closer to the LES data as was also expected from the analytical results in Fig. 11.

Eigendecomposition of the Reynolds stress tensors can be used to describe the turbulence state through its three real eigenvalues, which describe the fluctuations in the three orthogonal, principal directions. Several techniques (e.g., eigenvalue map, invariant map, barycentric map, and Lumley triangle) combine the eigenvalues and visualize them with two-dimensional maps; in Fig. 13 the barycentric map is used with the RGB color scheme of Emory and Iaccarino (2014). Both RANS models use ASL inflow; hence the turbulence state is the same at all heights, whereas the LES state varies with height due to its PDBL inflow. The k–ε–f_P turbulence is mostly isotropic (the x_3c corner) as is expected for all Boussinesq-type models because normal stresses are $\stackrel{\mathrm{‾}}{u_{\mathit{\alpha }}^{\prime }{u}_{\mathit{\alpha }}^{\prime }}=\frac{\mathrm{2}}{\mathrm{3}}k$, whereas the WJ-EARSM turbulence is perturbed more towards 2D turbulence (the line connecting x_2c and x_1c) and thereby closer to LES turbulence.

Figure 13 RGB-colored barycentric triangle, with RANS/LES inflow data shaded by height; white lines at right mark (z_ref−R), z_ref, and (z_ref+R), respectively.

Figure 14 Streamwise velocity (a–c) and TI (d–f) contours at hub height for the single-wake V80 case. Full lines mark ADs, while dashed lines mark where wake profiles are extracted.

4.2 Velocity and turbulence intensity

Wake data in the form of velocity and TI contours at hub height and profiles at three downstream positions are shown in Figs. 14 and 15. In addition to the 2D WJ-EARSM, results using the 3D WJ-EARSM with/without tangential AD forces are also shown in Fig. 15, and the first conclusion to draw from this is that the 2D and 3D versions of the WJ-EARSM give similar wake profiles, and we therefore focus on the simpler 2D WJ-EARSM in the following. A simple diffusion correction to the WJ-EARSM was suggested by Wallin and Johansson (2000) to correct the model in regions with low normalized velocity gradients (e.g., at the top of the half-channel), but it only has a small effect on velocity deficit and TI, as seen in Fig. 15, so it is not considered in the following discussions. Finally, results of the standard k–ε model are also shown in Fig. 15 to show its overly diffusive behavior.

Figure 15 Streamwise velocity (a–c) and TI (d–f) profiles extracted at various downstream positions for the single-wake V80 case.

Overall the wake velocity contours in Fig. 14 appear similar, while the TI contours of the WJ-EARSM are improved over the k–ε–f_P model; however, similar to the k–ε–f_P model, it still fails to predict the TI delay in the near wake seen in the LES, which is also clearly visible in the disk-averaged TI recovery profiles in Fig. 16b. Also, there is a TI induction zone in both RANS simulations, which is not present in the LES; these two effects seem to be a general issue of k–ε-based RANS models as this was also observed with the standard k–ε, realizable k–ε, and re-normalization group (RNG) k–ε models by Hornshøj-Møller et al. (2021).

Figure 16 Disk-averaged streamwise velocity and TI for the single-wake V80 case.

Considering the wake profiles in Fig. 15, it is clear that WJ-EARSM produces more “top-hat-shaped” profiles similar to the classic Jensen model (Jensen, 1983). This has also been seen with other EARSMs (van der Laan, 2014) and with full differential RSMs (Cabezón et al., 2011; Tian et al., 2019). One can either apply the previously mentioned diffusion correction or increase the c₁ constant to obtain a more Gaussian-shaped profile for the WJ-EARSM, but the latter is not recommended as it will deteriorate the ASL anisotropy; see Appendix B. Another view is that the top-hat-shaped profile is a consequence of not taking various physical phenomena into account, e.g., large-scale atmospheric turbulence and wake instabilities, and applying a unidirectional wind direction in our idealized RANS setup (in the transient LES there will be a varying instantaneous wind direction throughout the simulation). These types of effects could be interpreted as a Gaussian filter on the wake profiles; see discussion in Appendix C. Note that the k–ε–f_P model already accounts for at least the wind direction effect as it has been calibrated with LES, which includes a wind direction distribution (van der Laan, 2014), and that it could be recalibrated to obtain a better match with the velocity deficit of the present LES case (since the latter LES was simulated with a different solver and AD implementation).

Figure 17 Normal stresses at hub height for the single-wake V80 case.

Figure 18 Lateral U momentum flux at hub height x–y plane (a–c) and vertical U momentum flux at the center x–z plane (d–f) for the single-wake V80 case.

Lastly, we want to emphasize that the turbulence model is not solely responsible for the wake results: the same turbulence model applied for the same case but with different codes/solvers can yield significantly different results, as can be seen in the comparison of the k–ε–f_P results in Fig. 16. This reminds us to be careful with general conclusions on which is the better turbulence model.

4.3 Stresses

For more insights on the wake mixing and turbulence, we now turn to the second-order statistics of turbulence, namely the individual Reynolds stress components. The normal components are shown in Fig. 17, which shows that the WJ-EARSM correctly dampens the lateral and vertical components, which are otherwise overestimated severely by the k–ε–f_P model due to its Boussinesq closure. This also explains the lower TI of the WJ-EARSM in Fig. 14–16 as the TI is composed of the sum of the normal components. On the other hand, the normal stress contours seem too elongated in the streamwise direction with the WJ-EARSM; for example the low $\stackrel{\mathrm{‾}}{u^{\prime }{u}^{\prime }}$ core in the center of the wake in Fig. 17b extends too far downstream compared to the LES data in Fig. 17c, which is possibly connected with the decreased turbulence mixing also causing the top-hat profiles of wake deficit shown in Fig. 15. Increasing the c₁ constant will alleviate this specific issue; see Appendix B.

Transport of U momentum by turbulence, $-\partial \stackrel{\mathrm{‾}}{u^{\prime }{u}_j}/\partial x_j$, from the ambient high-speed surroundings to the low-speed wake region is mainly determined by the off-diagonal components of the Reynolds stress tensor, $\stackrel{\mathrm{‾}}{u^{\prime }{v}^{\prime }}$ and $\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}$, i.e., the lateral and vertical turbulent fluxes of U, respectively, which are shown as contours in Fig. 18. Both are larger in absolute magnitude for the k–ε–f_P model compared to WJ-EARSM, leading to larger gradients of the stresses and explaining the increased wake recovery of the former, but when comparing to the LES data, it can be noticed that the magnitude of $\stackrel{\mathrm{‾}}{u^{\prime }{w}^{\prime }}$ is overestimated, hence the overestimated wake recovery of the k–ε–f_P model in Fig. 16. On the other hand, the WJ-EARSM underestimates the magnitude of the fluxes, especially notable in the $\stackrel{\mathrm{‾}}{u^{\prime }{v}^{\prime }}$ contours, hence the slight underestimation of wake recovery compared to LES in Fig. 16.

Table 7 Simulation parameters for the aligned-row case from TotalControl. The AD scaling and 1D mom'm control methods are described in detail by van der Laan et al. (2015a).

To conclude, we see some advantages but also disadvantages with using WJ-EARSM over the k–ε–f_P model for prediction of the Reynolds stresses, which also has direct consequences for the prediction of the velocity deficit and turbulence intensity.

4.4 Turbulence state

In the inflow section, Sect. 4.1, it is shown that the turbulence state of the freestream turbulence was mainly isotropic in the k–ε–f_P model, while both the WJ-EARSM and LES data were more perturbed towards 2D turbulence (the lower edge in the barycentric triangle). This is also seen in the RGB-colored x–y plane at hub height in Fig. 19, where one can see that the ambient flow of the k–ε–f_P model is predominantly colored blue and is hence isotropic. As also alluded to in Sect. 4.1, the ambient WJ-EARSM turbulence is more oblate (left edge of barycentric triangle), while the LES is more prolate (right edge of barycentric triangle), hence the green and purple coloring, respectively, of the ambient flow in Fig. 19. In the wake shear layers, both RANS models predict perturbations towards 2D turbulence, while LES is more perturbed towards 1D turbulence (the x_1c corner). These observations fit with the increase in $\stackrel{\mathrm{‾}}{u^{\prime }{u}^{\prime }}$ and strong damping of $\stackrel{\mathrm{‾}}{v^{\prime }{v}^{\prime }}$ and $\stackrel{\mathrm{‾}}{w^{\prime }{w}^{\prime }}$ in the LES shown in Fig. 17. To improve the prediction of the turbulence state, we therefore suspect that normal stress predictions are essential, especially the ratios of those – the WJ-EARSM definitely has some improvement from the overly isotropic k–ε–f_P model but does not capture the completely right ratio between normal stresses. We note that a resolution of 20 cells per diameter was used in the reference LES, which might be sufficient for first-order statistics and shear stresses but could bias the normal stresses and therefore also the turbulence state.

Figure 19 RGB-colored turbulence componentiality at hub height for the single-wake V80 case.

5 Aligned row of wind turbines

To test the WJ-EARSM in a wind farm scenario, we simulate the lower row in the TotalControl rot90 reference wind farm, which consists of eight aligned wind turbines with 5 D inter-spacing; see Fig. 20. The DTU 10 MW turbines of the wind farm are modeled with the scaling AD method with 1D momentum control (van der Laan et al., 2015a) and including tangential forces (hence there will be wake swirl), and the thrust, rotational speed, and power curves are taken from the DTU 10 MW report (Bak et al., 2013). We compare the results with LES conducted at KU Leuven (the “PDk90 case”; Sood and Meyers, 2020), where the turbines were modeled with an actuator surface (AS) model coupled to an aeroelastic code; hence the LES and RANS are not directly comparable but nevertheless give a reference to compare against. Due to the natural streaks appearing in LES, we choose to calculate U_ref and I_ref with planar averages of the LES data in the region upstream of the lower row rather than using the time-averaged precursor profiles. The overview of the simulation parameters is given in Table 7. The numerical setup of the RANS simulation is identical to the one used for the single-wake case, except that tangential forces and non-uniform thrust forces are applied on the AD, the domain is scaled with the new rotor diameter, and the wake region is extended to encompass all eight turbines.

Figure 20 A x–y cut of the domain with every eighth cell shown (a) and a 3D view of the wake domain (b) for the aligned-row case. The size of the domain is $L_x/D=\mathrm{104}$, $L_y/D=\mathrm{60}$, and $L_z/D=\mathrm{25}$. A rough-wall BC and an inlet BC are used for the lower and upper BCs, respectively.

Figure 21 Streamwise velocity contour at hub height for the aligned-row TotalControl case.

Figure 22 TI contour at hub height for the aligned-row TotalControl case.

The velocity contours in Fig. 21 show that both the k–ε–f_P model and WJ-EARSM qualitatively share some of the same flow features as the LES; e.g., a faster wake recovery is seen on the left side of the first wake (seen from upstream) because of the combined effect of wind shear and wake rotation (a.k.a. swirl), and they predict the largest wake deficit in the second wake. Neither of these effects are predicted by the standard k–ε model, but it can capture other phenomena as for example the induction zones and the expanding wake tube surrounding the whole row. Figure 23a shows the streamwise velocity along the axial line of the ADs and shows that the WJ-EARSM is closer to the LES in the near wake of each turbine, while the k–ε and k–ε–f_P models are better at the far wake of each turbine. The RANS models become more similar further down in the row, and the recovery behind the last turbine (a.k.a. the wind farm recovery) is very similar. In van der Laan et al. (2015c), a similar observation was made when comparing the standard k–ε and k–ε–f_P models for wind farm cases since the difference between both turbulence models reduces with increased levels of turbulence.

The turbulence intensity contours in Fig. 22 show more pronounced differences between the models; e.g., the LES has peaks of TI in the wake shear layers, whereas the RANS models have TI more evenly distributed over the wake. As is also seen and discussed in the V80 case, there is no induction zone of TI for the LES, and the development of TI is delayed. Again also the WJ-EARSM has lower TI compared to the k–ε and k–ε–f_P models and is in better agreement with the LES.

Figure 23 Streamwise velocity and turbulence intensity at the axial line going through the AD centers for the aligned-row TotalControl case.

In conclusion, the WJ-EARSM appears numerically stable and well behaved (e.g., no monotonic decreasing velocity deficit or other unphysical effects) for interacting wakes, and as in the single-wake case, there are both some improvements and some less desirable effects of the model over the k–ε–f_P model.

6 Conclusions

This paper documents and explains our implementation and application of an EARSM (Wallin and Johansson, 2000) as a turbulence model for RANS simulations of wind turbine wakes, in the neutral ASL. To our knowledge, EARSM is rarely – if ever – used in the wind energy community; but we show that it is actually relatively straightforward to implement in CFD codes that already employ two-equation turbulence models and importantly that the WJ-EARSM also appears to be numerically stable for wake simulations. Previous attempts by van der Laan (2014) of applying EARSMs (models of Taulbee, 1992; Gatski and Speziale, 1993; and Apsley and Leschziner, 1998) for wake simulations showed problems with numerical stability, even for single-wake cases, but this appears to not be an issue for the EARSM of Wallin and Johansson (2000). The reason for the better numerical behavior most likely lies in the self-consistent formulation of particular importance in regions with rapid shear, hence preserving physical realizability.

Three canonical flow cases – homogeneous shear flow, half-channel flow, and square duct flow – were used to verify the implementation of the model and also showcased some of the advantages with an EARSM over traditional linear EVMs, namely the prediction of freestream turbulence anisotropy and secondary flow phenomena. All three cases have either analytical asymptotes or DNS data to compare against and are easy to set up, which makes them ideal for verification purposes.

For neutral ASL inflow we show that there is a delicate relationship between the turbulence constants that needs to be fulfilled to ensure a non-developing freestream solution and that it also dictates the amount of freestream turbulence anisotropy. It should be noted that this balance of constants is also important for numerical robustness. Comparing the RANS inflow with reference LES data shows that the WJ-EARSM is capable of predicting similar freestream anisotropy, whereas the turbulence of the k–ε–f_P model is nearly isotropic (although not exactly) by definition of its Boussinesq closure. This is also clear from the eigendecomposition of the Reynolds stress tensor, which was visualized with the barycentric map technique.

A single-wake case was considered first, and it was observed that the 2D version of the EARSM yielded almost identical results to the 3D version, even when tangential forces were applied on the AD; thus we use the 2D model for the remainder of the paper. It should be noted that the 2D version of the WJ-EARSM is a complete and invariant model for general three-dimensional mean flows. Only the particular dependency of pure three-dimensional effects is simplified, which will have a minor impact on most three-dimensional mean flows of interest. The wake profiles of the EARSM were more top-hat-shaped than the profiles observed in the LES data, which might be related to the underlying weak-equilibrium assumption and limitations in the length-scale-determining ε model equation but could also be caused by wind direction variations (see Appendix C). The disk-averaged velocity deficit and turbulence intensity recovery profiles were however improved over the k–ε–f_P model for the specific case considered. It is also possible to obtain a more Gaussian-shaped wake profile by increasing the c₁ constant and re-tuning the turbulence model constants, but this will deteriorate the prediction of the underlying ASL anisotropy (Appendix B). A notable difference between the WJ-EARSM and the k–ε–f_P model is also that the latter predicts large peaks of lateral and vertical normal Reynolds stress components in the wake, which are not present in the LES data, because of the Boussinesq closure; this deficiency and the low value of C_μ used for the k–ε–f_P model are possibly the reasons why it tends to overpredict turbulence intensity in the wake.

Finally, we simulated a row of eight aligned turbines, where the trends from the single-wake case could also be seen, e.g., the top-hat-shaped profiles and better turbulence intensity prediction with EARSM. There were more uncertainties in the comparison with the LES data in this case because different turbine modeling techniques were used, but the case nevertheless shows that the EARSM also behaves sensibly in cases with wake–wake interaction in the sense that the code still converges and that no unphysical trends (such as monotonic increase in wake deficit throughout the row of turbines) are observed.

In conclusion, the EARSM of Wallin and Johansson (2000) can be used for wake simulations in a numerically robust way and only has a small computational overhead compared to standard two-equation models (on the order of 5 % for the V80 single-wake case); hence it is at least 3 orders of magnitude faster than LES. It provides an advantage over two-equation models in the sense that it has more realistic inflow with anisotropic turbulence and that the wake turbulence also becomes more anisotropic, which indeed is also observed in LES. The turbulence intensity prediction was improved for both test cases considered, while velocity deficit was only considerably improved for the single-wake case; more cases (both LES and experimental data) are needed to draw general conclusions about its performance in this regard.

Atmospheric conditions in thermally stable stratification and thermal convection strongly influence the turbulence states, anisotropies, and in particular the vertical mixing in the ASL. This will have a fundamental influence on the wake development and the performance of wind parks. The extension of the EARSM to non-neutral conditions has over recent years been developed by Lazeroms et al. (2013), and Želi et al. (2019) have subsequently demonstrated the model's capability of capturing these effects. This will be of interest for future wind turbine wake studies.

Appendix A: Grid study of WJ-EARSM for wind turbine wake simulations

Earlier studies (van der Laan et al., 2015b) have shown that a grid spacing of $D/\mathrm{8}$ in the wake region is sufficient for grid convergence of wake velocity deficits with the k–ε–f_P model, and the same conclusion can be drawn for the 2D WJ-EARSM; see Fig. A1. We also plot the TI profiles in Fig. A1, which are more sensitive to grid resolution, but we nevertheless decide to use $D/\mathrm{8}$ in this paper because it represents the typical resolution used in wind farm studies and saves a considerable amount of computational resources.

Figure A1 Grid study of streamwise velocity (upper row) and TI (lower row) profiles at hub height for the V80 case with the 2D WJ-EARSM using different mesh resolutions.

Appendix B: Tuning the turbulence model constants in the WJ-EARSM

The c₁ constant, a.k.a. Rotta coefficient, in the WJ-EARSM originates from the pressure-redistribution term and can in principle be re-tuned; e.g., the original LRR (Launder–Reece–Rodi; Launder et al., 1975) model has c₁=1.5, and van der Laan (2014) used c₁=4.5 for the f_P model. As shown in Sect. 4.1, for the neutral ASL the WJ-EARSM can be used to obtain a direct relationship between c₁ and many other variables in equilibrium; see Fig. 11. From this it is clear that the ASL anisotropy is enhanced for decreasing c₁ and vice versa. It is important to emphasize that the variables in Fig. 11 are only dependent on c₁ and not the other turbulence constants. In this section we test three different c₁ values in addition to the standard value of c₁=1.8; see Table B1.

Table B1 Tested sets of turbulence model constants and derived variables for the single-wake V80 case. For all sets, we use C_ε2=1.82, σ_k=1.0, σ_ε=1.3, κ=0.38, and $c_{\mathrm{2}}=\frac{\mathrm{5}}{\mathrm{9}}$.

Figure B1 Inflow profiles for the single-wake V80 case with different sets of model constants.

Equation (23) needs to be satisfied to have a balanced RANS solution; hence when c₁ is adjusted, and C_μ thereby changes, then either κ, σ_ε, C_ε1, C_ε2, or a combination of all needs to be adjusted. We choose to adjust C_ε1 and fix the others, and the resulting sets of constants are shown in Table B1. The first set (c₁=1.2) gives anisotropic freestream turbulence close to the Panofsky and Dutton (1984) neutral ASL values; the second set (c₁=1.8) is equivalent to the standard set (see Table 2); the third set (c₁=4.0) gives C_μ=0.054, which is close to the value used by WAsP CFD (Bechmann, 2016); and finally the last set (c₁=8.0) gives C_μ=0.03, which is often used for atmospheric applications (Sørensen, 1995; Richards and Hoxey, 1993).

Figure B2 Streamwise velocity (a–c) and TI (d–f) profiles extracted at various downstream positions for the single-wake V80 case with different sets of model constants.

Figure B3 Normal stresses at hub height for the single-wake V80 case with different sets of model constants.

Another consideration is that the roughness length and friction velocity also need to be modified to give the same hub height velocity and turbulence intensity according to Eq. (22), again because C_μ changes with changing c₁. The values of these are therefore also included in Table B1 and explain why the velocity inflow profiles differ slightly in Fig. B1. The figure also clearly demonstrates that freestream turbulence anisotropy decreases for increasing c₁ and vice versa, which is also evident from the combination of Eqs. (10)–(11):

$$\begin{array}{}\text{(B1)}& {\displaystyle \frac{{\mathit{\beta }}_{\mathrm{4}}}{{\mathit{\beta }}_{\mathrm{1}}}}=N^{-\mathrm{1}}={\left({\displaystyle \frac{\mathrm{9}}{\mathrm{4}}}{c}_{\mathrm{1}}\right)}^{-\mathrm{1}}.\end{array}$$

In the derivation of Eq. (B1) we use the definition of N and $c_{\mathrm{1}}^{\prime }$ as well as $\mathcal{P}/\mathit{\epsilon }=\mathrm{1}$ in the neutral ASL. From Eq. (B1), we see that for increasing c₁ there is a decreasing $a_{ij}^{\left(\mathrm{ex}\right)}$ (see Eq. 9), which is the part of the closure responsible for anisotropy and therefore explains why the turbulence becomes more isotropic with increasing c₁.

Figure B2 shows how the wake is effected by the new sets of constants. The velocity deficit shape is more Gaussian for larger c₁ and thus more similar to the LES shape, while the turbulence intensity increases. From these observations, one could argue that c₁=4.0 would perhaps be a better choice for modeling of velocity deficit, while c₁=1.8 is better for TI and anisotropy predictions.

Increasing c₁ also has a significant impact on the normal stress contours as shown in Fig. B3, where one especially can note how the length of the inner low $\stackrel{\mathrm{‾}}{u^{\prime }{u}^{\prime }}$ core decreases. In this regard a larger c₁ than the standard c₁=1.8 also seems desirable, although we again have to remind the reader that this will also result in less correct freestream anisotropy. However, it is notable that although a larger c₁ leads to more isotropic freestream turbulence, in contrast to the k–ε–f_P model (see Fig. 17), there is still significant wake turbulence anisotropy.

Appendix C: Effect of wind direction uncertainty

Figure C1 (a) Gaussian wind direction distribution and (b) the Gaussian filter width for the convolution with σ_φ=2.8^∘.

As discussed in Sect. 4.2 and Appendix B, the WJ-EARSM velocity deficit is rather top-hat-shaped, whereas the LES shape is more Gaussian. This could possibly be a consequence of the unidirectional inflow used in our steady RANS simulation and hence the lack of wake meandering (Larsen et al., 2008), which otherwise tends to smear out the time-averaged wake. For the k–ε–f_P model it was effectively included in the turbulence model through the calibration of the C_R model constant, but for the WJ-EARSM there is less room for calibration, and instead one should therefore ideally do a post-processing step to include the effect. It is however beyond the scope of this paper to carefully design such a post-processing step, but we can at least qualitatively show that the effect of wind direction uncertainty is to smoothen the top-hat-shaped wake profile into a more Gaussian shape.

Figure C2 Gaussian-filtered (GF) wake profiles.

The wind direction $\mathit{\phi }\equiv \arctan \left(\frac{V}{U}\right)\approx \frac{V}{U}$ (small angle approximation) variance is

$$\begin{array}{}\text{(C1)}& \begin{array}{rl}{\mathit{\sigma }}_{\mathit{\phi }}^{\mathrm{2}}& ={\left({\displaystyle \frac{\partial \mathit{\phi }}{\partial U}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_u^{\mathrm{2}}+{\left({\displaystyle \frac{\partial \mathit{\phi }}{\partial V}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_v^{\mathrm{2}}\\ & ={\left(-{\displaystyle \frac{V}{U^{\mathrm{2}}}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_u^{\mathrm{2}}+{\left({\displaystyle \frac{\mathrm{1}}{U}}\right)}^{\mathrm{2}}{\mathit{\sigma }}_v^{\mathrm{2}}\approx {\displaystyle \frac{{\mathit{\sigma }}_v^{\mathrm{2}}}{U^{\mathrm{2}}}}.\end{array}\end{array}$$

For example the hub height LES values for the V80 case gives ${\mathit{\sigma }}_{\mathit{\phi }}=\pm \mathrm{0.0484}$, rad = ±2.8^∘. However, according to Larsen et al. (2008) the wake meandering motion is only caused by the large eddies (the “slow-moving” part of the wind direction changes), so one should in principle use a smaller variance to model the meandering motion. However, as most of the variance is given by the large eddies, the usage of ${\mathit{\sigma }}_{\mathit{\phi }}^{\mathrm{2}}$ is only expected to give a small overprediction of meandering.

Assuming zero mean wind direction and a Gaussian distribution gives a simple model for the wind direction variability:

$$\begin{array}{}\text{(C2)}& f_{\mathit{\phi }}\left(\mathit{\phi }\right)={\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathit{\phi }}\sqrt{\mathrm{2}\mathit{\pi }}}}{e}^{-\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{\mathit{\phi }}{{\mathit{\sigma }}_{\mathit{\phi }}}\right)}^{\mathrm{2}}}.\end{array}$$

It should be intuitively clear that such a distribution of wind directions would act to smoothen out the wake profile compared to a simulation with φ=0^∘, but we can also illustrate this by applying a convolution to the φ=0^∘ wake profile:

$$\begin{array}{}\text{(C3)}& \stackrel{\mathrm{̃}}{U}(x,y)=\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}U(x,y-y^{\prime })f_y(x,y^{\prime })\mathrm{d}{y}^{\prime }.\end{array}$$

To obtain f_y from f_φ one could simply make a change in variables from φ to y and assume σ_y=xσ_φ, but the latter assumption is equivalent to assuming that the φ=0^∘ wake is moved “rigidly” from side to side, which is not realistic for at least two reasons. First, it assumes “frozen turbulence” in the sense that a wake released in the φ direction continues in that direction until infinity. Secondly, the yaw mechanism of a turbine is slower than the turbulence timescale; hence the wind direction uncertainty corresponds to yaw offsets, which will bend the wake inwards. It must therefore be presumed that σ_y increases more slowly than linearly in x and could for example be modeled as a quadratic polynomial instead:

$$\begin{array}{}\text{(C4)}& {\displaystyle }& {\displaystyle }{f}_y(x,y)={\displaystyle \frac{\mathrm{1}}{{\mathit{\sigma }}_y\sqrt{\mathrm{2}\mathit{\pi }}}}{e}^{-\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{y}{{\mathit{\sigma }}_y}\right)}^{\mathrm{2}}}\text{(C5)}& {\displaystyle }& {\displaystyle \frac{{\mathit{\sigma }}_y}{D}}=\left\{\begin{array}{ll}\left(\frac{x}{D}-\mathrm{0.05}{\left(\frac{x}{D}\right)}^{\mathrm{2}}\right){\mathit{\sigma }}_{\mathit{\phi }},& \mathrm{0}\le x/D\le \mathrm{10}\\ \mathrm{5}{\mathit{\sigma }}_{\mathit{\phi }},& x/D\ge \mathrm{10}\end{array}\right..\end{array}$$

The above model for σ_y was chosen such that it approaches the linear filter width model for $x/D\to \mathrm{0}$, has $\mathrm{d}{\mathit{\sigma }}_y/\mathrm{d}x\ge \mathrm{0}$, and is continuous; see Fig. C1. The constant −0.05 is rather arbitrary, and it must again be stressed that this is just a qualitative model to model the two aforementioned effects (one could just as well have used any other simple function to model σ_y).

Both the linear and quadratic filter width models are used in Gaussian convolutions of the 2D WJ-EARSM wake profiles, and the results are shown in Fig. C2. Although σ_φ=2.8^∘ at first might appear as a small number, a significant smoothing effect is seen, which effectively removes the top-hat-shaped velocity deficits and shows that missing meandering could actually explain the top-hat shape. The linear width model overestimates the smoothing as expected from its inherent frozen turbulence assumption, while the quadratic model gives results closer to the LES. The TI results are not improved by the Gaussian convolution correction.