Inflow modeling for wind farm flows in RANS

Wind turbine interaction in wind farms can lead to energy loss and increased wind turbine loads, with the magnitude of these effects strongly depending on atmospheric conditions. One-dimensional Reynolds-averaged Navier-Stokes (RANS) models are able to represent the atmospheric boundary layer (ABL) over a range of model fidelity, and can be used for steady-state inflow profiles in three-dimensional RANS simulations of wind farms. In the present work, an overview of existing and recently developed atmospheric inflow models is presented. The inflow models are applied to separately demonstrate impacts on the velocity deficit of a row of ten wind turbines, due to: turbulence intensity and atmospheric stability in the atmospheric surface layer; ABL depth; and Coriolis-induced veer.


Introduction
Modeling wind turbine interaction in wind farms is important to obtain a better understanding of wake-induced energy losses and blade fatigue loads. In addition, wind turbine interaction can also lead to energy redistribution between wind turbines due to upstream effects known as global blockage. The magnitude of the wake and blockage effects in wind farms are strongly dependent on the atmospheric conditions as turbulence intensity and atmospheric stability. High-fidelity simulations methods as Computational Fluid Dynamics (CFD) can be employed to simulate atmospheric wind farm flows. Large-eddy simulation (LES) is a transient CFD method that resolves the large scale turbulence and models the small scales. LES is a popular method in academia but it is computationally too expensive for most industrial and practical applications because of the need for a fine grid resolution and time step. Reynolds-averaged Navier-Stokes (RANS) is a steady-state CFD method that simulates the mean statistics directly and it is roughly three to four orders of magnitude faster than LES. RANS models all turbulence scales and may not provide accurate results in all cases compared to turbulence resolving method as LES; however, it is a useful method to isolate relative effects, which can leverage the understanding of wind farm flows.
One-dimensional RANS models are able to represent the idealized atmospheric boundary layer (ABL) by a range of model fidelities, which can be used as steady-state inflow profiles for three-dimensional RANS simulations of wind farms. Here, idealized means neglecting meso-scale effects as horizontal pressure gradients (baroclinicity) and transient phenomena. Since RANS is a steady-state CFD method, it requires that the inflow is in balance with the momentum and turbulence model equations for a set of boundary conditions that corresponds to a flow domain with a homogeneous surface roughness and elevation (without wind turbines). If this requirement is not met, then the inflow will develop downstream and the wind farm will experience a different 2 inflow compared to the inflow set at the inlet boundary. One could place a wind farm relatively close to an inlet to reduce the downstream development of the inflow; however, this can influence the wind farm flow because the wind farm induction could be affected. Similar arguments can be made for atmospheric simulations in related scientific areas as flows in urban areas, complex terrain and forests. Therefore, atmospheric turbulence models have been developed that are meant to be in equilibrium with a chosen inflow as for example discussed by Richards and Hoxey (1993) [1], Alinot and Masson (2005) [2] and Parente et al. (2011) [3], and a type of atmospheric inflow will often correspond to a certain turbulence model. Atmospheric turbulence models in RANS are mostly based on either the k-ε or k-ω two-equation turbulence models. These models solve a transport equation for the turbulence kinetic energy k and a transport equation for the dissipation of k called ε or ω. The two turbulence quantities define the turbulence velocity scale and the turbulence length scale from which the eddy viscosity ν T is determined. In the present work, we choose the k-ε model, but it is also possible to use the k-ω model as discussed by Sogachev et al. (2012) [4]. It is important to note that the use of an active temperature equation may lead to a transient inflow model, especially for unstable conditions. An atmospheric inflow model including a temperature equation could be considered an unsteady RANS (URANS) model and this type of inflow model is not discussed in the present work. Furthermore, the effect of an active temperature equation on a single wake is minor, as shown by El-Askary et al. (2017) [5], RANS-based atmospheric inflow models can be subdivided into two types; an atmospheric surface layer (ASL) model and an ABL model. The ASL covers roughly the first 10% of the ABL. A widely used ASL inflow model is a logarithmic inflow, which neglects effects of atmospheric stability and Coriolis, and assumes a constant Reynolds-stress for all heights. A more general ASL model can be obtained following Monin-Obukhov Similarity Theory (MOST) [6], which includes effect of atmospheric stability expressed by analytic functions of the wind shear and potential temperature gradient. MOST assumes both a constant Reynolds-stress and vertical turbulent heat flux for all heights. The non-neutral conditions are modeled by source terms of buoyancy in the turbulence transport equations and a temperature equation is not necessary by relating the buoyancy to the wind shear. As a result, the MOST inflow model is independent of the Reynolds number for a fixed stability condition [7]. The analytic profiles following MOST are not a solution of the k-ε model and modifications are required when using MOST in combination with the k-ε model, as discussed by Alinot and Masson (2005) [2] and previous work [8]. The latter is used in the present work and it is referred to as the k-ε MO model. Note that the k-ε MO model is the same as the standard k-ε model for a neutral ASL when the stability parameter is set to resemble neutral conditions. Analytic inflow models beyond the surface layer do currently not exist unless an unrealistic constant or linearly varying eddy viscosity is assumed for all heights [9]. ABL models based on the k-ε model are numerical solutions of the one-dimensional RANS and turbulence transport equations. The ABL models include an ABL depth that can be modeled by the global turbulence length scale limiter of Apsley and Castro (1997) [10] and a temperature equation is not necessary. One of the ABL models discussed in this work is forced by a balance between a constant pressure gradient and the Coriolis force, and we refer to this model as the k-ε ABLc model, where c stands for Coriolis. The k-ε ABLc model includes effects of wind veer and can model effects of stable conditions by lowering the ABL depth. The k-ε ABLc model cannot model unstable conditions without additional modifications as discussed in previous work [9], where a source term of buoyancy based on MOST was added in the turbulence transport equations. The correction for unstable conditions increases the turbulence length scale near the surface beyond the neutral ASL value (as expected for unstable stratification) and compares well with measurements of the unstable ASL [9]. However, the obtained ABL depth is very large (O(10 km)) and the ABL model needs to be further developed to represent a more realistic convective ABL (beyond the 3 ASL). Therefore, we are not investigating the k-ε ABLc model including the unstable correction in the present work. While the k-ε ABLc model is dependent on four dimensional parameters the non-dimensional ABL profiles are only dependent on two Rossby numbers, one based on the roughness length and another one based on the maximum turbulence length scale that is used to force an ABL depth [9].
A pressure-driven ABL model has recently been developed [11], which is a simplification of the k-ε ABLc model by removing the interaction between the horizontal momentum equations, which results in zero wind veer. The model is labelled as the k-ε ABLp model and it can be used to follow both Reynolds (Re)-and Rossby number (Ro) similarity [11]. While k-ε ABLp model may not be realistic, it can be used to isolate the effect of wind veer when comparing with the k-ε ABLc model because both models can be employed to predict the same turbulence intensity and turbulence length scale. In previous work [12,13], the effect of wind veer has been investigated with a more traditional pressure-driven ABL model; however, such a model is more challenging for isolating the effect of wind veer because it does not predict a similar turbulence intensity and turbulence length scale as the k-ε ABLc model.

Model
Type Solution Stability Coriolis Wind veer Similarity Table 1. Overview of idealized atmospheric inflow turbulence models investigated in the present article, ranked by model fidelity.
In the present work, the four different atmospheric inflow and turbulence models are further discussed in Section 2, an overview is provided in Tab. 1. Here, the atmospheric stability is defined as neutral (N), stable (S) and unstable (U). The inflow turbulence models are applied to separately demonstrate the impact of atmospheric turbulence intensity and stability in the atmospheric surface layer, ABL depth, and wind veer on a row of ten wind turbines in Section 4 using a methodology presented in Section 3.

Atmospheric inflow and turbulence models in RANS
The atmospheric inflow models of Tab. 1 are a solution of incompressible RANS equations including the linear stress-strain relation of Boussinesq (1890) [14]: Here, U , V and W are the velocity components defined in streamwise (x), lateral (y) and vertical (z) directions, respectively. The inflow turbulence models are based on a homogeneous horizontal flow and they can be obtained from a one-dimensional version of incompressible RANS equations, where only the height z is the remaining spatial dimension and the momentum equation in the z direction is redundant. The general three-dimensional form is also applicable to wind farm flows and it can be written as: with ρ as the density, P as the pressure, x i = {x, y, z} as the Cartesian coordinates and S i = {S U , S V , S W } as a momentum source term. We have neglected the molecular viscosity ν in eq. (1) since ν << ν T , but one could simply add ν to ν T if desired. The eddy viscosity is where C µ is the eddy viscosity coefficient, f P is scalar function that reduces the eddy viscosity or turbulence length scale locally where the normalized shear, k/ε U 2 i,j , is much larger than the normalized shear of the inflow [15]. The f P function improves the velocity deficit of the near wake but has little to no influence for the inflow models since f P ≈ 1 or f P = 1 in the freestream depending on which inflow model is used. The f P function is further discussed in Section 2.4. The transport equations of k and ε are defined as: where C * ε,1 and C ε,3 are variable turbulence model parameters, P is the turbulence production due to shear, B is a source or sink of turbulence due to buoyancy, and S k and S ε are additional source terms of k and ε, respectively. The definition of the source terms and variable turbulence model parameters are summarized in Tab. 2 and they are further discussed in Sections 2.1, 2.2 and 2.3. The following model constants are employed (C µ , σ k , σ ε , C ε,2 , κ) = (0.03, 1, 1.3, 1.92, 0.4) and C ε,1 = C ε,2 − κ 2 /(σ ε C µ ) ≈ 1.21 is used to make the k-ε model in balance with the logarithmic law. The value of C µ is based on field measurements of the ASL [16] following Richards and Hoxey (1993) [1].
The input parameters of the atmospheric turbulence models are summarized in Tab. 3, and are further discussed in Sections 2.1, 2.2 and 2.3.

k-ε and k-ε MO models
The inflow profiles of k-ε MO model are based on MOST [6], which represents the effects of surface layer atmospheric stability as empirically derived analytic functions. These functions represent a normalized wind shear, φ m = (κz/u * )dU/dz, a normalized dissipation, φ ε = (κz/u 3 * )ε, and a normalized potential temperature gradient, φ h = (κz/θ * )dθ/dz, with θ as the potential temperature and θ * as the surface layer temperature scale. MOST assumes both

Model
Non-dimensional input Set input to define the shape Derived input Table 3. Input parameters of atmospheric inflow turbulence models.
a constant Reynolds-stress and vertical turbulent heat flux for all heights, and neglects effects of Coriolis. The profiles can be written as where u * is the friction velocity, z 0 is the roughness length and Ψ m is the integral of φ m . The functions φ m , φ ε , φ h and Ψ m are only dependent on a stability parameter ζ ≡ (z + z 0 )/L, with L as the Obukhov length, a full description is provided in previous work [8]. Note that we choose to use (z + z 0 ) instead of z because we employ a numerical implementation of the rough wall boundary condition where the roughness length is placed on top of the wall boundary [17]. For ζ = 0, we get Ψ m = 0, φ m = φ ε = 1, and the neutral ASL is obtained, where U is a logarithmic profile and k is a constant.
The k-ε MO model represents the effect of stability by a buoyancy source term B in both transport equations. The definition of the buoyancy source term is dependent on a temperature gradient but it is possible to rewrite it as function of the vertical wind shear, and stability functions φ h and φ m following MOST, the result is listed in Tab. 2. The latter is a convenient model choice because a temperature variable and temperature equation are not necessary and the k-ε MO model becomes independent of the Reynolds number for a fixed stability parameter [7]. In other words, the shape of inflow profiles is independent of U ref and the magnitude can be scaled by U ref , as indicated in Tab. 3.
The inflow profiles of eq. (5) are not a solution of the standard k-ε model plus buoyancy. This can be solved by an additional source term S k in the transport equation of k and using a C ε,3 parameter that is function of ζ, as shown in previous work [8]. The source term S k counteracts the non-zero diffusion in the k-equation, which is mainly active for unstable conditions.
The MOST inflow of eq. (5) is defined by the dimensional parameters z 0 , u * and L, where z 0 and L define the shape and u * the magnitude.
Subsequently, one needs to choose z ref and U ref to determine z 0 and u * .

k-ε ABLp model
The recently developed k-ε ABLp model [11] is a simple pressure-driven ABL model that includes an ABL depth but neglects Coriolis-induced wind veer. The latter is accomplished by removing the interaction between the U -and V -momentum equations. The model is forced by a constant geostrophic wind speed G = {U G , V G } representing a constant horizontal pressure gradient, and it is balanced by a remainder of the Coriolis forces in absence of turbulence: , with f pg as a forcing parameter that can be related to the Coriolis parameter f c . The model is not meant to represent a realistic idealized ABL, where wind veer is present; however, it can be employed to isolate the effect of ABL depth or wind veer when comparing its results with results of the k-ε or k-ε ABLc models, respectively.
The ABL depth is indirectly set by the global length scale limiter of Apsley and Castro (1997) [10], C * ε,1 , see Tab. 2. Here, is the local turbulence length scale, ≡ C 3/4 µ k 3/2 /ε, and max is its maximum set value, which is a proxy for ABL depth [9]. For low values of max (e.g. 5-10 m), a small ABL depth is obtained and the turbulence intensity and turbulence length scale are reduced, which could represent stable atmospheric conditions (without wind veer). For large values of max (e.g. 30-100 m), a taller ABL depth is obtained that could resemble neutral conditions. It is not possible to model unstable conditions without additional modeling [9].
Above the ABL depth, the turbulence variables k and ε decay to zero for increasing z, and this can lead to numerical problems in RANS. One can avoid this by the addition of small ambient sources terms of k and ε, using S k and S ε as listed in Tab. 2. The ambient source terms are dependent on k amb and ε amb , which we parametrize as [9]: with amb as an ambient turbulence length scale, C amb as a model constant that defines the ratio amb / max and I amb as an ambient turbulence intensity above the ABL. We choose C amb = 10 −6 and I amb = 10 −6 based on previous work [11]. The k-ε ABLp model has four dimensional input parameters: namely, the geostrophic wind speed G, the roughness length z 0 , the forcing parameter f pg and the maximum turbulence length scale max . However, it can be shown that the normalized results of numerical solutions of the k-ε ABLp model are only dependent on two Rossby numbers based on the forcing parameter f pg , each with a different length scale [11]: Here, Ro 0 is a surface layer Rossby number and Ro is a Rossby number related to the maximum turbulence length that defines the ABL depth. The Rossby similarity can be used to create a library of numerical solutions of the k-ε ABLp model, which can then be utilized to find an ABL profile that has the desired characteristics as I ref and U ref at z ref [11]. We employ the k-ε ABLp model in two different ways, which is indicated by the roman number in Tab. 3 (I or II). If one chooses to set a constant Ro 0 (by varying the forcing parameter f pg for different values of G such that G/f pg is constant), then the inflow model becomes independent of the Reynolds number for a fixed I ref (set by z 0 and max ), and the shape of the ABL profiles are independent of the wind speed similarly to the ASL models. The Reynolds similarity can be used to perform parametric studies of wind farms more quickly, as discussed in previous work [18]. In addition, one can investigate the effect of ABL depth on a wind farm by changing the maximum turbulence length scale max for a constant roughness length and Ro 0 . The second usage of the model, k-ε ABLp II in Tab. 3, is to isolate the effect of wind veer by comparing its results to results of the k-ε ABLc model. This is done by using the same values of max and z 0 , while using different values of G and f pg for the k-ε ABLp model in order to get the same I ref and U ref at z ref as predicted by the k-ε ABLc model.

k-ε ABLc model
The k-ε ABLc is the same as the k-ε ABLp model; however, the Coriolis-induced wind veer is included by using realistic momentum sources: where f c is the Coriolis parameter. The k-ε ABLc model also follows a Rossby similarity [9]: where the Rossby numbers are based on f c .

Coupling of inflow turbulence models to a wake turbulence model via the f P function
The standard k-ε model cannot represent the turbulence length scale of the near wind turbine wake and several researchers have proposed extended k-ε models to overcome this [19][20][21]. In previous work, we have developed such a model, called the k-ε-f P model [15], where a scalar f P is multiplied by the eddy viscosity [eq. (2)], which limits the turbulence length scale in the near wake [22]. The k-ε-f P model has been developed and calibrated for a neutral ASL using LES. In more recent works, the k-ε-f P model has been coupled with the k-ε MO model [23] and k-ε ABLc model [12], although these couplings have not yet been fully validated and may need to be revised. The f P function is defined as is the shear parameter of inflow assuming MOST. For an empty flow domain (without wind turbines) over a homogeneous rough flat surface, f P = 1 . The model constant C R is set to 4.5 based on calibration using LES for eight different single wake cases with neutral ASL inflow. When the f P function is coupled with one of the k-ε ABL turbulence models (k-ε ABLp or k-ε ABLc), the global / max term within C * ε,1 (defined in Table 2) [10,12] can be multiplied by a blending function f 1 (f P ) to switch off the global length scale limiter in regions where the local turbulence limiter f P is active. We use a hyperbolic tangent form (e.g. App. B of [24]): where w = 0.02 is the characteristic width of the 'switching' and f P,sw = 0.9 is the f P value where the switching happens.

Methodology of numerical simulations
The four inflow models of Section 2 are applied to RANS simulations of ten in-line NREL-5MW wind turbines [25] with a spacing of 7D, where the rotor diameter D is set to 126 m. The reference values of the inflow models are based on z ref equal to the hub height of 90 m, and a rowaligned wind direction is used. The RANS simulations are carried out with EllipSys3D, the inhouse flow solver of DTU Wind Energy (Michelsen, 1992 [26]; Sørensen, 1995 [27]). EllipSys3D is an incompressible finite volume code employing block-structured curvilinear grids. A onedimensional version of EllipSys3D [28] is used to generate inflow profiles for the k-ε ABLp and k-ε ABLc models. The wind farm simulation setup is nearly identical to that used in previous work [29]; a brief summary is provided here. The wind turbines are modeled as actuator disks [30] [32]. This method is known to overestimate the wind turbine forces by a few percent [33]; however, here we are only interested in the trends of the wake deficits subjected to different inflow models. The numerical domain is a Cartesian grid of about 100 km in both horizontal directions. The wind turbine row is placed in the center of the domain, where a refined horizontal cell size of D/10 is applied; this extends 3D upstream, 20D downstream and 3D in the lateral directions. The grid spacing is based on previously performed grid study [15]. At the ground, a rough wall boundary condition is applied [17]; above that the cells grow from a lowest cell height of 0.5 m, up to a spacing of D/10 at z = 3D. For z > 3D the cells continue to grow, with an expansion ratio of 1.2, up to a height of 25D. In total 14.7 million cells are used. An inflow boundary condition is set at the inlet and top of the domain, and all gradients across the outflow boundary are prescribed to be zero. Periodic conditions are set at the lateral boundaries.
We define eight inflow cases, with the input parameters listed in Tab. 4. Cases 1-3 are used to investigate the effect of ambient turbulence intensity under neutral conditions by employing the standard k-ε model. The effect of stable (Case 4) and unstable (Case 5) surface layer stability is studied using the k-ε MO model; the turbulence intensity is based on Cases 1 and 3, respectively, which allows a comparison with neutral conditions. Cases 6 and 7 are meant to isolate the effect of ABL depth, which is set by changing the turbulence intensity from 6 to 4%, for a constant roughness length of 5 × 10 −3 m. Finally, Case 8 is used to investigate the effect of wind veer when comparing with Case 7. A low ABL depth is chosen for Case 8, to obtain a strong wind veer over the rotor area. In addition, in Case 7 the surface-layer Rossby number Ro 0 is chosen to match the ABL profiles from Case 8; we use the same Ro 0 for Case 6 to isolate the effect of ABL depth when comparing Cases 6 and 7. The resulting inflow profiles are further discussed in the Section 4.

Results and Discussion
The eight inflow cases listed in Tab. 4 are employed to study the effect of turbulence intensity, surface layer atmospheric stability, ABL depth and wind veer on the wake deficit of a row of ten wind turbines using the wind farm simulation setup as discussed in Section 3. The inflow profiles are depicted in Fig. 1. Figures 1d and 1h show that the standard k-ε model can only predict a single turbulence length scale, namely = κz, while the k-ε MO model can represent larger and smaller values for unstable and stable conditions, respectively. The k-ε ABL models always predict a smaller turbulence scale compared to the k-ε model (as shown by comparing Cases 6, 7 and 8 with Case 2), and therefore the ABL profiles represent stable conditions. Cases 4 and 6 actually have a similar atmospheric stability since one can relate the maximum turbulence length scale to a stable Obukhov length: max ≈ 0.08L [10], which results in L = 181 m for Case 6. The wind shear (Figure 1e) of the inflow models is related to the behavior of the turbulence length scale, with the exception of the standard k-ε model. Figures 1b and 1f show that all models have a zero wind veer, with the exception of the k-ε ABLc model, which predicts a 13.3 • change of wind direction over the rotor area. This large amount of wind veer should be considered as an extreme case that can occur for very stable conditions. Cases 7 and 8 show nearly identical profiles of turbulence intensity and turbulence length scale as intended, which makes them ideal to isolate the effect of wind veer. Small differences between Cases 7 and 8 are obtained for the wind speed in the top part of the rotor area due to the lack of the supergeostrophic jet in k-ε ABLp model. The latter is a result from neglecting the Coriolis-induced interaction of the U and V momentum equations [11].  Figure 2 depicts the streamwise velocity integrated over a fictitious rotor area,Ū , normalized by its value without wind turbine forcesŪ 0 , as function of downstream distance. The inflow turbulence models are coupled with the f P function, as discussed in Section 2.4. The effect of the ambient turbulence intensity under neutral conditions is shown in Figure 2a. When increasing the ambient turbulence intensity the wakes recover faster and the peak of the velocity deficit of the first wind turbine moves upstream. The effect of the ambient turbulence intensity is less important further downstream because of the wake generated turbulence. Figure 2b and 2c depicts the effect of stable and unstable conditions, respectively, compared to neutral conditions with the corresponding ambient turbulence intensity. As expected, the  Fig. 2c also shows a similar conclusion for unstable conditions, which is not physical. The k-ε MO model is made in balance with MOST using an additional source term in the kequation, S k , as discussed in Section 2. S k is mainly active for unstable conditions [8] and acts as a sink, but is typically an order of magnitude smaller than the other terms in eq. (3), so it cannot explain the unphysical behavior. The unphysical behavior of the k-ε MO model for unstable conditions is instead related to an over-production of ε connected with our formulation of B or the lack of turbulence transport, which need to be revised in future work.
The effect of ABL depth on the velocity deficits is depicted in Fig. 2d. A shallow ABL is charecterized by a lower turbulence intensity and turbulence length scale compared to a taller ABL. Hence, the wake recovery is slower with decreasing ABL depth, which is a similar trend obtained for the wake recovery for neutral and stable conditions, as shown in Fig. 2b. Figure 2e depicts the effect of wind veer on the wake deficits. When a strong wind veer is present (Case 8), the wakes deflect clockwise as seen from above for the Northern hemisphere [13], and the downstream wind turbines are not fully aligned with the local wind direction, which reduces velocity deficit seen by the wind turbines compared to an inflow without wind veer.

Conclusions
An overview of four existing atmospheric inflow turbulence models in RANS is presented; the standard k-ε model for a neutral ASL, the k-ε MO model for a non-neutral ASL, the k-ε ABLp model for an ABL without wind veer and the k-ε ABLc model for an ABL with Coriolis-induced wind veer. The inflow turbulence models are coupled with a wake turbulence model and they are applied to RANS simulations of a row of ten wind turbines with an aligned wind direction to separately demonstrate the effect of ambient turbulence intensity, atmospheric stability, ABL depth and wind veer on the velocity deficit. Reducing the ambient turbulence intensity increases the velocity deficit and moves the peak velocity deficit of the first wind turbine wake downstream. Stable atmospheric conditions also increase the velocity deficit as shown by the k-ε MO model and the k-ε ABLp model by reducing the ABL depth. While the effect of turbulence intensity is less pronounced further downstream in the wind turbine row due the wake generated turbulence, the effect of stable conditions persist for entire wind turbine row. The k-ε MO model is also used to investigate the effect of unstable conditions and unphysical trends are obtained due to our buoyant production formulation. The unphysical model behavior is related to an inconsistent over-production of ε and the lack of turbulence transport modeling, which needs to be investigated in future work. The effect of strong wind veer is shown to reduce the velocity deficit by comparing results of the k-ε ABLp and k-ε ABLc models, due to fact that the wind veer deflects the wind turbine wakes. The accuracy of the velocity deficit predicted by the RANS models will be investigated in future work, employing both LES and measurements.