RANS modeling of a single wind turbine wake in the unstable surface layer
- DTU Wind Energy, Technical University of Denmark, Risø Campus, Frederiksborgvej 399, 4000 Roskilde, Denmark
- DTU Wind Energy, Technical University of Denmark, Risø Campus, Frederiksborgvej 399, 4000 Roskilde, Denmark
Correspondence: Mads Baungaard (mchba@dtu.dk)
Hide author detailsCorrespondence: Mads Baungaard (mchba@dtu.dk)
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.
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.
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^{3} 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:
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^{−1} (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^{−1} (Kelly and Gryning, 2010), but offshore there tends to be a bias towards more unstable conditions, i.e., more negative L^{−1} compared to onshore (Sathe et al., 2013). Various parametrizations have been suggested for wind speed, k and ε in terms of ζ; in this paper we use the widely accepted forms of Dyer (1974) for U (namely the Ψ_{m} and Φ_{m} functions) and those found in Kaimal and Finnigan (1994) for ε and k (see van der Laan et al., 2017, for more details). Under unstable conditions these are
where
and
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^{−1}, which is not rare), we nevertheless still choose to apply the profiles and in fact use them all the way up to the top boundary. More realistic inflow profiles for RANS covering the whole atmospheric boundary layer (ABL) are indeed a current research topic but will not be discussed further in this paper. Maronga and Reuder (2017) reason that MOST is a pragmatic solution because the parameters needed for more realistic inflow profiles are often not available.
The roughness length z_{0} 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}):
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.
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^{−2} (±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.
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}.
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_{0} 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.
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):
where
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}.
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}.
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{\u203e}}{{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:
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^{−1} is used for all plots in this section.
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{\u203e}}{{\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{\u0303}}{\mathit{\sigma}}$, is however not equal for the two cases (see inlets in third row), which means that $\mathit{\sigma}/\stackrel{\mathrm{\u0303}}{\mathit{\sigma}}$ differs a lot between the two cases. This observation is important because $\mathit{\sigma}/\stackrel{\mathrm{\u0303}}{\mathit{\sigma}}$ is the main parameter used in the f_{P} correction to be discussed in the next section.
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
As shown in Fig. 7, $\mathit{\sigma}/\stackrel{\mathrm{\u0303}}{\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{\u0303}}{\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):
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_{0} parameter is also stability-dependent, i.e.,
This is actually the form suggested by Apsley and Leschziner (1998), but they were not considering stability effects, i.e., no variation of $\stackrel{\mathrm{\u0303}}{\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{\u0303}}{\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
where $\stackrel{\mathrm{\u0303}}{B}/\stackrel{\mathrm{\u0303}}{\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{\u0303}}{\mathit{\sigma}}$, modification 1 and modification 2 is shown in Fig. 8. When plotted with $\mathit{\sigma}/\stackrel{\mathrm{\u0303}}{\mathit{\sigma}}$ on the abscissa, the stability-dependent $\stackrel{\mathrm{\u0303}}{\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{\u0303}}{\mathit{\sigma}}>\mathrm{1}$, i.e., in the wake. This is necessary to compensate for the larger $\mathit{\sigma}/\stackrel{\mathrm{\u0303}}{\mathit{\sigma}}$ encountered with non-neutral inflow; see Fig. 7.
When both modifications are used, faster wake recovery for a given I_{ref} occurs in unstable conditions, as shown in Fig. 9, i.e., similar behavior as when f_{P} was “turned off” (i.e., f_{P}=1); see Fig. 6.
The cstB model with the f_{P} modifications described in the previous section is tested for the cases summarized in Table 2. Each case is simulated with a range of C_{B} parameters, ${C}_{\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.
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.
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^{−1} 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^{−1} 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.
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^{−1}, 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.
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.
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.
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^{−2}) on the AD are calculated in each iteration as
Here, f_{n,ij} and f_{θ,ij} are the normal and azimuthal surface force distributions at the (i,j)th AD element (i: radial direction, j: azimuthal direction), which are applied in the CFD domain using the methodology described by Réthoré et al. (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_{0} is the normalized circulation, g is the root correction and F is the tip correction. The latter two are obtained with Delery's root correction (with parameters a=2.335, b=4.0 and $\stackrel{\mathrm{\u203e}}{\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{\u203e}}{\mathit{\delta}}$'s, in Fig. A2. Clearly, the Joukowsky AD with $\stackrel{\mathrm{\u203e}}{\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.
EllipSys3D and PyWakeEllipSys are proprietary software of DTU. Information about the latter can, however, be freely accessed at https://topfarm.pages.windenergy.dtu.dk/cuttingedge/pywake/pywake_ellipsys/ (van der Laan, 2022).
The RANS results were generated with DTU's proprietary software, but the data presented can be made available by contacting the corresponding author. Interested parties are also welcome to hand-digitize the results and use them as reference in other publications.
MB performed the RANS simulations and proposed the modifications to the turbulence model. All authors (MB, MPvdL and MK) contributed to derivation of the new model and article writing.
The contact author has declared that neither they nor their co-authors have any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The authors are grateful to the creators of the LES and lidar results used in the validation section. We would also like to thank the two reviewers for their feedback and suggestions.
This paper was edited by Sandrine Aubrun and reviewed by two anonymous referees.
Abkar, M. and Porté-Agel, F.: The effect of atmospheric stability on wind-turbine wakes: A large-eddy simulation study, J. Phys.: Conf. Ser., 524, 1–9, https://doi.org/10.1088/1742-6596/524/1/012138, 2014. a
Abkar, M. and Porté-Agel, F.: Influence of atmospheric stability on wind-turbine wakes: A large-eddy simulation study, Phys. Fluids, 27, 035104, https://doi.org/10.1063/1.4913695, 2015. a, b, c, d
Alinot, C. and Masson, C.: Aerodynamic simulations of wind turbines operating in atmospheric boundary layer with various thermal stratifications, ASME 2002 Wind Energy Symposium, WIND2002, 206–215, https://doi.org/10.1115/wind2002-42, 2002. a, b
Apsley, D. D. and Castro, I. P.: A limited-length-scale k–ϵ model for the neutral and stably-stratified atmospheric boundary layer, Bound.-Lay. Meteorol., 83, 75–98, https://doi.org/10.1023/A:1000252210512, 1997. a
Apsley, D. D. and Leschziner, M. A.: A new low-Reynolds-number nonlinear two-equation turbulence model for complex flows, Int. J. Heat Fluid Flow, 19, 209–222, https://doi.org/10.1016/S0142-727X(97)10007-8, 1998. a
Chougule, A., Mann, J., Kelly, M., and Larsen, G. C.: Simplification and Validation of a Spectral-Tensor Model for Turbulence Including Atmospheric Stability, Bound.-Lay. Meteorol., 167, 371–397, https://doi.org/10.1007/s10546-018-0332-z, 2018. a
Churchfield, M. J., Lee, S., Michalakes, J., and Moriarty, P. J.: A numerical study of the effects of atmospheric and wake turbulence on wind turbine dynamics, J. Turbulence, 13, 1–32, https://doi.org/10.1080/14685248.2012.668191, 2012. a, b, c, d, e
Crespo, A., Manuel, F., Moreno, D., Fraga, E., and Hernandez, J.: Numerical analysis of wind turbine wakes, in: Proceedings of the Delphi Workshop on Wind Energy Applications, 15–25, https://www.academia.edu/15318618/Numerical_Analysis_of_Wind_Turbine_Wakes (last access: 30 March 2022), 1985. a
Doubrawa, P., Quon, E. W., Martinez‐Tossas, L. A., Shaler, K., Debnath, M., Hamilton, N., Herges, T. G., Maniaci, D., Kelley, C. L., Hsieh, A. S., Blaylock, M. L., Laan, P., Andersen, S. J., Krueger, S., Cathelain, M., Schlez, W., Jonkman, J., Branlard, E., Steinfeld, G., Schmidt, S., Blondel, F., Lukassen, L. J., and Moriarty, P.: Multimodel validation of single wakes in neutral and stratified atmospheric conditions, Wind Energy, 23, we.2543, https://doi.org/10.1002/we.2543, 2020. a, b, c, d, e, f
Dyer, A. J.: A review of flux-profile relationships, Bound.-Lay. Meteorol., 7, 363–372, 1974. a
El-Askary, W. A., Sakr, I. M., AbdelSalam, A. M., and Abuhegazy, M. R.: Modeling of wind turbine wakes under thermally-stratified atmospheric boundary layer, J. Wind Eng. Indust. Aerodynam., 160, 1–15, https://doi.org/10.1016/j.jweia.2016.11.001, 2017. a, b
Foken, T.: 50 years of the Monin-Obukhov similarity theory, Bound.-Lay. Meteorol., 119, 431–447, https://doi.org/10.1007/s10546-006-9048-6, 2006. a
Ghaisas, N. S., Archer, C. L., Xie, S., Wu, S., and Maguire, E.: Evaluation of layout and atmospheric stability effects in wind farms using large-eddy simulation, Wind Energy, 20, 1227–1240, https://doi.org/10.1002/we.2091, 2017. a
Gryning, S. E., Batchvarova, E., Brümmer, B., Jørgensen, H., and Larsen, S.: On the extension of the wind profile over homogeneous terrain beyond the surface boundary layer, Bound.-Lay. Meteorol., 124, 251–268, https://doi.org/10.1007/s10546-007-9166-9, 2007. a
Han, X. X., Liu, D. Y., Xu, C., Shen, W. Z., Li, L. M., and Xue, F. F.: A BEM-based actuator disk model for wind turbine wakes considering atmospheric stability, Preprint, 1–24, https://doi.org/10.20944/preprints201907.0029.v1, 2019. a
Hancock, P. E. and Zhang, S.: A Wind-Tunnel Simulation of the Wake of a Large Wind Turbine in a Weakly Unstable Boundary Layer, Bound.-Lay. Meteorol., 156, 395–413, https://doi.org/10.1007/s10546-015-0037-5, 2015. a
Hansen, K. S., Barthelmie, R. J., Jensen, L. E., and Sommer, A.: The impact of turbulence intensity and atmospheric stability on power deficits due to wind turbine wakes at Horns Rev wind farm, Wind Energy, 15, 183–196, https://doi.org/10.1002/we.512, 2012. a, b
Hansen, M. O. L.: Aerodynamics of Wind Turbines, 3rd Edn., Routledge, https://doi.org/10.4324/9781315769981, 2015. a
Jonkman, J., Butterfield, S., Musial, W., and Scott, G.: Definition of a 5-MW Reference Wind Turbine for Offshore System Development, NREL technical report, NREL, https://doi.org/10.2172/947422, 2009. a, b, c
Kaimal, J. and Finnigan, J. J.: Atmospheric Boundary Layer Flows, Oxford University Press, https://doi.org/10.1093/oso/9780195062397.001.0001, 1994. a
Kasmi, A. E. and Masson, C.: An extended k2 model for turbulent flow through horizontal-axis wind turbines, J. Wind Eng. Indust. Aerodynam., 96, 103–122, https://doi.org/10.1016/j.jweia.2007.03.007, 2008. a
Keck, R. E., de Maré, M., Churchfield, M. J., Lee, S., Larsen, G., and Aagaard Madsen, H.: On atmospheric stability in the dynamic wake meandering model, Wind Energy, 17, 1689–1710, https://doi.org/10.1002/we.1662, 2014. a, b, c, d, e, f, g
Kelly, M. and Gryning, S. E.: Long-Term Mean Wind Profiles Based on Similarity Theory, Bound.-Lay. Meteorol., 136, 377–390, https://doi.org/10.1007/s10546-010-9509-9, 2010. a, b
Li, N., Liu, Y., Li, L., Chang, S., Han, S., Zhao, H., and Meng, H.: Numerical simulation of wind turbine wake based on extended k–ϵ turbulence model coupling with actuator disc considering nacelle and tower, IET Renewable Power Generation, https://doi.org/10.1049/iet-rpg.2020.0416, 2020. a
Machefaux, E., Larsen, G. C., Koblitz, T., Troldborg, N., Kelly, M. C., Chougule, A., Hansen, K. S., and Rodrigo, J. S.: An experimental and numerical study of the atmospheric stability impact on wind turbine wakes, Wind Energy, 19, 1785–1805, https://doi.org/10.1002/we.1950, 2016. a, b, c, d, e, f, g
Magnusson, M. and Smedman, A. S.: Influence of Atmospheric Stability on Wind Turbine Wakes, Wind Eng., 18, 139–152, 1994. a
Maronga, B. and Reuder, J.: On the formulation and universality of Monin–Obukhov similarity functions for mean gradients and standard deviations in the unstable surface layer: Results from surface-layer-resolving large-eddy simulations, J. Atmos. Sci., 74, 989–1010, https://doi.org/10.1175/JAS-D-16-0186.1, 2017. a
Michelsen, J. A.: Basic3D : A platform for development of multiblock PDE solvers, Technical Report AFM 92-05, Technical University of Denmark, Lyngby, 1992. a, b
Monin, A. S. and Obukhov, A. M.: Basic laws of turbulent mixing in the surface layer of the atmosphere, Tr. Akad. Nauk SSSR Geophiz. Inst., 24, 163–187, 1954. a
Okulov, V. L. and Sørensen, J. N.: Maximum efficiency of wind turbine rotors using Joukowsky and Betz approaches, J. Fluid Mech., 649, 497–508, https://doi.org/10.1017/S0022112010000509, 2010. a
Panofsky, H. and Dutton, J.: Atmospheric Turbulence, John Wiley & Sons, Ltd, ISBN 0471057142, 1984. a
Pope, S. B.: Turbulent Flows, Cambridge University Press, https://doi.org/10.1017/CBO9780511840531, 2000. a
Porté-Agel, F., Wu, Y. T., Lu, H., and Conzemius, R. J.: Large-eddy simulation of atmospheric boundary layer flow through wind turbines and wind farms, J. Wind Eng. Indust. Aerodynam., 99, 154–168, https://doi.org/10.1016/j.jweia.2011.01.011, 2011. a, b
Porté-Agel, F., Bastankhah, M., and Shamsoddin, S.: Wind-Turbine and Wind-Farm Flows: A Review, vol. 174, Springer Netherlands, https://doi.org/10.1007/s10546-019-00473-0, 2020. a
Rados, K. G., Prospathopoulos, J. M., Stefanatos, N. C., Politis, E. S., Chaviaropoulos, P. K., and Zervos, A.: CFD modeling issues of wind turbine wakes under stable atmospheric conditions, in: European Wind Energy Conference and Exhibition 2009, EWEC 2009, 6, Marseille, 4085–4092, ISBN 9781615677467, https://www.researchgate.net/publication/237819589_CFD_modeling_issues_of_wind_turbine_wake_under_stable_atmospheric_conditions, (last access: 30 March 2022), 2009. a, b
Réthoré, P.-E.: Wind Turbine Wake in Atmospheric Turbulence, PhD thesis, Aalborg University, Aalborg, https://orbit.dtu.dk/en/publications/wind-turbine-wake-in-atmospheric-turbulence (last access: 30 March 2022), 2009. a
Réthoré, P.-E., van der Laan, P., Troldborg, N., Zahle, F., and Sørensen, N. N.: Verification and validation of an actuator disc model, Wind Energy, 17, 919–937, https://doi.org/10.1002/we.1607, 2014. a
Sathe, A., Mann, J., Barlas, T., Bierbooms, W., and van Bussel, G.: Influence of atmospheric stability on wind turbine loads, Wind Energy, 16, 1013–1032, https://doi.org/10.1002/we.1528, 2013. a
Sørensen, J. N. and Andersen, S. J.: Validation of analytical body force model for actuator disc computations, J. Phys.: Conf. Ser., 1618, 1–8, https://doi.org/10.1088/1742-6596/1618/5/052051, 2020. a
Sørensen, J. N. and Kock, C. W.: A model for unsteady rotor aerodynamics, J. Wind Eng. Indust. Aerodynam., 58, 259–275, https://doi.org/10.1016/0167-6105(95)00027-5, 1995. a
Sørensen, J. N., Nilsson, K., Ivanell, S., Asmuth, H., and Mikkelsen, R. F.: Analytical body forces in numerical actuator disc model of wind turbines, Renew. Energy, 147, 2259–2271, https://doi.org/10.1016/j.renene.2019.09.134, 2020. a, b
Sørensen, N. N.: General purpose flow solver applied to flow over hills, PhD thesis, Technical university of Denmark, Risø National Laboratory, https://orbit.dtu.dk/en/publications/general-purpose-flow-solver-applied-to-flow-over-hills (last access: 30 March 2022), 1995. a, b, c
Sørensen, N. N., Bechmann, A., Johansen, J., Myllerup, L., Botha, P., Vinther, S., and Nielsen, B. S.: Identification of severe wind conditions using a Reynolds Averaged Navier-Stokes solver, J. Phys.: Conf. Ser., 75, 012053, https://doi.org/10.1088/1742-6596/75/1/012053, 2007. a
Stull, R. B.: An introduction to boundary layer meteorology, An introduction to boundary layer meteorology, Kluwer Academic Publishers, https://doi.org/10.1007/978-94-009-3027-8, 1988. a
Troldborg, N., Sørensen, N. N., Réthoré, P.-E., and van der Laan, M. P.: A consistent method for finite volume discretization of body forces on collocated grids applied to flow through an actuator disk, Comput. Fluids, 119, 197–203, https://doi.org/10.1016/J.COMPFLUID.2015.06.028, 2015. a, b
van der Laan, M. P.: Efficient Turbulence Modeling for CFD Wake Simulations, PhD thesis, DTU Wind Energy, Technical University of Denmark, https://orbit.dtu.dk/en/publications/efficient-turbulence-modeling-for-cfd-wake-simulations (last access: 30 March 2022), 2014. a, b, c, d, e, f
van der Laan, M. P.: PyWakeEllipSys documentation, https://topfarm.pages.windenergy.dtu.dk/cuttingedge/pywake/pywake_ellipsys/, last access: 30 March 2022. a
van der Laan, M. P., Sørensen, N. N., Réthoré, P.-E., Mann, J., Kelly, M. C., and Troldborg, N.: The k-ϵ-fp model applied to double wind turbine wakes using different actuator disk force methods, Wind Energy, 18, 2223–2240, https://doi.org/10.1002/we.1816, 2015a. a, b
van der Laan, M. P., Sørensen, N. N., Réthoré, P. E., Mann, J., Kelly, M. C., Troldborg, N., Schepers, J. G., and Machefaux, E.: An improved k-ϵ model applied to a wind turbine wake in atmospheric turbulence, Wind Energy, 18, 889–907, https://doi.org/10.1002/we.1736, 2015b. a, b, c
van der Laan, P., Sørensen, N. N., Réthoré, P.-E., Mann, J., Kelly, M. C., Troldborg, N., Hansen, K. S., and Murcia, J. P.: The k-ϵ-fp model applied to wind farms, Wind Energy, 18, 2065–2084, https://doi.org/10.1002/we.1804, 2015c. a
van der Laan, M. P., Kelly, M. C., and Sørensen, N. N.: A new k-ϵ model consistent with Monin-Obukhov similarity theory, Wind Energy, 20, 479–489, https://doi.org/10.1002/we.2017, 2017. a, b, c, d, e, f, g, h, i, j
van der Laan, M. P., Andersen, S. J., Kelly, M., and Baungaard, M. C.: Fluid scaling laws of idealized wind farm simulations, J. Phys.: Conf. Ser., 1618, 1–10, https://doi.org/10.1088/1742-6596/1618/6/062018, 2020. a, b, c
van der Laan, M. P., Baungaard, M., and Kelly, M.: Inflow modeling for wind farm flows in RANS, J. Phys.: Conf. Seri., 1934, 012012, https://doi.org/10.1088/1742-6596/1934/1/012012, 2021. a, b, c, d, e, f
Wyngaard, J. C.: Turbulence in the atmosphere, Cambridge University Press, https://doi.org/10.1017/CBO9780511840524, 2010. a, b
Xie, S. and Archer, C. L.: A Numerical Study of Wind-Turbine Wakes for Three Atmospheric Stability Conditions, Bound.-Lay. Meteorol., 165, 87–112, https://doi.org/10.1007/s10546-017-0259-9, 2017. a
Zhang, W., Markfort, C. D., and Porté-Agel, F.: Wind-Turbine Wakes in a Convective Boundary Layer: A Wind-Tunnel Study, Bound.-Lay. Meteorol., 146, 161–179, https://doi.org/10.1007/s10546-012-9751-4, 2013. a
- Abstract
- Introduction
- Simulation setup
- Modification of the k–ε–f_{P} model in the unstable ASL
- Validation with experiments and LES
- Conclusions
- Appendix A: Simulation details
- Code availability
- Data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Review statement
- References
- Abstract
- Introduction
- Simulation setup
- Modification of the k–ε–f_{P} model in the unstable ASL
- Validation with experiments and LES
- Conclusions
- Appendix A: Simulation details
- Code availability
- Data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Review statement
- References