Monin-Obukhov similarity theory for modeling of wind turbine wakes under atmospheric stable conditions: Breakdown and modifications

: Monin–Obukhov similarity theory (MOST) overestimates the mean vertical velocity gradient in some atmospheric stable conditions, i


Introduction
As a measure of turbulence exchanges in the atmospheric surface layer, atmospheric stability can significantly affect the wind turbine wake deficit and its recovery rate.In general, turbulence exchanges between the wake and atmosphere are depressed under stable conditions.High wake deficits and slow wake recovery, thus, are usually observed in the stable stratification boundary layer [1,2].As wakes play critical roles in wind farm energy production and the fatigue loads of wind turbines, there is an increasing interest in studying the effects of atmospheric stability on wakes and developing wake models for non-neutral conditions.
The impact of atmospheric stability on wakes is widely observed in wind tunnel measurements of small-scale models of wind turbines [3][4][5] and full-scale field experiments [1,2,[6][7][8].According to the wind tunnel measurements in Chamorror et al. [3], the stronger inlet velocity gradient in the stable case leads to a slightly stronger turbulence intensity and extends the region of enhanced turbulence intensity from a distance of approximately 4 to 5.5 rotor diameters to 3 and 6 rotor diameters downwind of the turbine location.In Zhang et al. [4], a 15% smaller velocity deficit at the wake center, a more rapid momentum recovery due to an enhanced radial momentum transport, and 20% higher peak turbulence intensity were observed in the unstable case, as compared to the wake of the same wind turbine under neutral conditions.According to the field experiments of wakes using LiDARs or masts, Magnusson et al. [6], Iungo et al. [7], and Han et al. [2] also found that the velocity deficit decreases more slowly under stable conditions and more quickly under unstable conditions.Those observations suggest that atmospheric stability should be considered for improved wake models and predictions of wind power harvesting.
Besides field measurements and wind tunnel experiments of wind turbine wakes under non-neutral conditions, numerical models based on the computational fluid dynamics (CFD) have also been used to investigate the effect of atmospheric stability on the wind turbine wakes [1,[9][10][11][12][13][14].High-fidelity large eddy simulations (LES) with the actuator line model [9,10] or the actuator disk model [1,11] are commonly used to study the structure and dynamics of wind turbine wakes in varying stability cases.This is mainly because the LES approaches allow capturing the near wake structure and resolving the interactions of tip vortices with large-scale eddies of the ambient flow and wake meandering.However, the high-fidelity approaches are computationally expensive for wind energy engineering applications.There are some efforts towards developing turbulence models for wind turbine wakes under non-neutral conditions in reynolds-averaged navier-Stokes (RANS) to reduce computational costs.
A widely used turbulence model for the thermal stratified boundary layer was developed by Alinot and Masson [15], where a coefficient of the buoyant terms in the transport equation of the turbulence kinetic energy dissipation is calibrated with atmospheric stability.However, van der Laan et al. [16] showed that this model cannot keep the flow homogeneity in a large domain under unstable conditions and thus developed a turbulence model consistent with monin-obukhov similarity theory (MOST).In [13], an additional buoyancy production based on the Richardson number, Ri, was added to the turbulent kinetic energy equation to model wind turbine wakes under stable conditions.El-Askary et al. [14] further applied this model with an additionally dissipation source term to investigate the wake behavior at different atmospheric stability conditions.Note that these turbulence models are based on MOST and could fail in modeling wind turbine wakes when MOST breaks down under stable conditions for the flux Richardson number R f > 0.25.According to the literature [17,18], the classical similarity functions that are widely used in MOST to determine the wind profile are only valid for R f < 0.25 and overestimate the velocity gradient for R f > 0.25.For a given aerodynamic roughness z 0 and a hub-height inflow velocity, this overestimated velocity gradient results in a smaller friction wind speed u * , and consequently a lower turbulence intensity, and could influence wake modeling of wind turbines.
To the best of our knowledge, the side effects of the breakdown of MOST on wake modeling under stable conditions have not been well investigated.The present paper aims to investigate these side effects through numerical simulations and to make a modification on MOST to eliminate these side effects.This modification proposes a new set of similarity functions based on field measurements to limit the velocity gradient in very stable conditions and introduce the new similarity functions into the turbulence model proposed by van der Laan et al. [16].The remainder of this paper is organized as follows.MOST is briefly described in Section 2. Models for wind turbine wakes under stable conditions are introduced in Section 3, which cover two actuator disk models based on the thrust coefficient and Blade Element Theory (BEM) calculations, and three turbulence models for stable conditions based on MOST.Section 3 also proposes a modified turbulence model that can be consistent with arbitrary similarity functions, e.g., the new similarity functions.The breakdown of MOST under stable conditions is experimentally investigated and new similarity functions are proposed and validated in Section 4. All test models are studied through numerical simulations: the simulation details are described in Section 5, and the results of the simulations are discussed in Section 6 and then concluded in Section 7.

Monin-Obukhov Similarity Theory
According to MOST [19], any dimensionless characteristics of the turbulence depends only on the dimensionless stability parameter ζ = z/L, where z is the height above the surface and the Obukhov length L that is defined by in which κ = 0.4, g is the gravity acceleration, Θ is the time-averaged potential temperature, θ is the fluctuation of the potential temperature, and u * = −u w is the friction speed, where u and w are the fluctuations of the longitudinal and vertical velocity components.
The flow similarity functions of the gradient of velocity, the gradient of potential temperature, the turbulence kinetic energy, and its dissipation can thus be defined as where U is the mean streamwise speed, θ * = −w θ /u * is the scaling temperature, C µ = 0.033, k is the turbulence kinetic energy (TKE), and ε is the TKE dissipation.According to the eddy viscosity hypothesis by Boussinesq [20], we have where ρ is the air density and the eddy viscosity is modeled in the k − ε model as [21] Combining the definitions of φ m , φ k , and φ ε with Equations ( 6) and ( 7) leads to Based on measurements of flows over flat terrain in the atmospheric surface layer, the classical similarity functions commonly used in literature are given as [22,23] where γ m = γ h = 16, β m = β h = 5 [23], and χ will be determined based on field measurements in Section 4. The classical similarity functions φ m and φ h are, however, only valid for the flux Richardson number R f below 0.25 [17].We thus proposed new functions, noted as φ m,exp and φ h,exp , in the full range of R f based on a field experiment in a wind farm.

Models for Wind Turbine Wakes under Stable Conditions
In the present work, the numerical simulation of wind turbine wakes under atmospheric stable conditions is based on RANS equations and the buoyancy effect due to atmospheric thermal stratifications is also taken into account.The following continuity, momentum, and energy equations are solved in wake simulations [15,24], ∂ ∂t where ρ is the air density, U i is the velocity component in the x i direction, p is the air pressure, µ is the laminar viscosity, µ t is the turbulent viscosity, S u,i is the momentum term source imposed by the wind turbine rotor in the x i direction, Θ is the potential temperature, Pr = 0.9 is the laminar Prandtl number, and α t ≡ µ t /Pr t is the turbulent heat conductivity, where Pr t is the turbulent Prandtl number.

Modeling of Wind Turbine
In this work, we introduce two kinds of actuator disk models: one is based on the thrust coefficient and another is based on BEM calculations to distribute forces through the disk.The second model, which provides more detailed information of the distributed forces caused by the rotor, is expected to capture near wake structure better than the first one.All simulations carried out in this work use the AD model based on BEM calculations in advance if the geometry information of the blades is available.

Actuator Disk Model Based on Thrust Coefficient (-CT)
In the origin actuator disk model, the momentum source term in the x i direction due to the thrust is uniformly distributed through the rotor: where T is the thrust, V disk is the disk volume, C T is the thrust coefficient, U ref is the upstream reference velocity at hub height, U ref, i is the component velocity U ref of in x i direction, ∆l is the disk depth, and the negative sign represents the drag effects of thrust on the flow.
For the upstream flow disturbed by a complex terrain or wind turbine wakes, U ref is unknown and difficult to evaluate from the local flow directly.According to the one-dimensional momentum theory and ignoring the velocity gradient, the reference velocity is a function of the rotor-averaged velocity U disk : U disk = (1 − a B )U ref (16) in which the induced factor a B [25] is related to the thrust coefficient C T by where a c = 1/3.In simulations, the disk-averaged velocity U disk is calculated by averaging the local velocity in the disk region, then is applied to estimate the upstream reference velocity U ref based on Equations ( 16) and (17) and the momentum source S u,i based on Equation (15), respectively.

Actuator Disk Model Based on Blade Element Method
In the BEM-based actuator disk model, the rotor plane consists of N actuator lines and each of the actuator lines is split into M element sections (Figure 1).The element section collects the local velocity and the rotor speed Ω to calculate the element force and applies this force in the neighbor cells of the element section.The reference velocity is firstly assessed from the disk-averaged velocity and is then applied to evaluate the rotor speed.By transforming the local velocity at the blade element into polar velocity components (u r , u θ , and u n ), the force of the blade element is where B is the number of the blades, c is the chord length, and ∆r is the length of the element section.The drag coefficient of the element section, C D , and its lift coefficient, C L , which are functions of the attack angle α, are estimated from XFOIL [26] and then corrected by three-dimensional rotational effects of the blades based on Du et al. [27].According to Figure 1, α = φ − (β + γ), where φ = arctan [u n /(Ωr + u θ )] is the flow angle, β is the blade installation angle, and γ is the pitch angle.The element force is distributed across neighbor cells.The force added to a cell is calculated by where s i the distance of the i-th element to the cell and s is the cut-off length scale that takes a value between 2 and 3 cell sizes [28].F tip and F hub are the Prandtl tip loss and hub loss functions [29]: where R is the rotor radius, R hub is the hub radius, and r is the radial distance of the element to the rotor center.

Turbulence Modeling
In this paper, we apply the k − ε turbulence model to close Equations ( 12)-( 14): ∂ ∂t where k is the turbulence kinetic energy (TKE), ε is the TKE dissipation, and C ε2 = 1.92.P and B are the TKE source production due to shear and buoyancy, respectively: B ≡ where u i is the fluctuation of the velocity component in x i direction.S ε,wake is a source term of TKE dissipation applied in the neighbor region of the turbine to correct the fast wake recovery of the standard k − ε model [30]: where C 4ε = 0.37.In this work, the methods of Alinot and Masson [15], the recently proposed model from M.P. van der Laan et al. [16], and the indirect model that cannot solve the energy equation applied in W.A. El-Askary et al. [14] are investigated, and the models are hereafter referred to as the AM, Laan, and El-Askary models, respectively.In the Laan model and its modified model (the proposed model), the TKE source term S ε,ABL and the coefficient C ε3 are functions of atmospheric stability and modeled to keep flow homogeneity under stable conditions.However, S ε,ABL is absent in the AM model and the El-Askary model, and C ε3 is set to −2.9 in the AM model [15] and set to 1 in the El-Askary model [14].
In the El-Askary model, the energy Equation ( 14) is not solved, and the TKE production due to buoyancy B is modeled as where The Laan model is derived from homogeneous steady flows over flat terrain where the kinematic viscosity µ can be ignored, and the k − ε model can be rewritten into normalized form by κz/u 3 * which yields C ε1 = 1.24.
Under stable conditions, we have where φ = ∂φ/∂ζ and φ = ∂ 2 φ/∂ζ 2 for φ ∈ {φ k , φ m , φ ε }.Equations ( 32) and ( 33) allow turbulence closure that is consistent with either the classical similarity functions in MOST (the Laan model) or any other similarity functions based on field measurements (the proposed model in this paper).One can prove that the above equations equal to those proposed by van der Laan et al. [16] if the classical similarity functions are applied.Note that the turbulent Prandtl number is set to 1 in the origin AM model and the Laan model, and the energy Equation ( 14) is also not solved in the Laan model.This work, however, applies Pr t ≡ φ h /φ m according to the definition of the turbulent Prandtl number [31] and solves the energy Equation ( 14) in the Laan model.There are negligible differences in wake modeling when considering these differences in implementations of the models.Details of the turbulence models investigated in this paper are summarized in Table 1.

Breakdown and Modifications of MOST under Stable Conditions
This section shows the breakdown of MOST in predicting the velocity gradient under stable conditions and proposes a set of new similarity functions to eliminate this side effect.All these works are based on a field experiment over complex terrain in a wind farm.The influence of complex terrain on estimating similarity functions is investigated and assessed for predicting wind profiles.

Experimental Data
The field experiment used to investigate the similarity functions was carried out in the Jingbian wind farm in northwest of China [2].In the experimental campaign, two masts, numbered M1 and M3, are installed near a wind turbine (No. 14) to capture the wake profiles (Figure 2).As the inflow from north is dominated by a small and constant terrain slope while the inflow from east and south are dominated by a significant terrain complexity, we applied measurements from the mast M1 to estimate similarity functions.On the mast M1, three Thies First Class Advanced cup anemometers were installed at 30 m, 50 m, and 70 m above the ground level (AGL) (Figure 3).Wind direction (WD) was measured at 30 m and 60 m AGL, and air temperature was measured at 30 m and 70 m AGL, respectively.Two Metek 3D sonic anemometers were installed on M1 to provide three-dimensional velocity components and sonic temperature at 30 m and 70 m AGL in a frequency of 35 Hz.A NRG logger was used to record the 10-min statistical values from the cup anemometers, wind vanes, and thermometers.The measured data are resampled to 10 min in ~190 days.All the anemometers are calibrated in wind tunnels to compensate flow distortion effects.The classification of cup anemometers after calibrations is class A 1.5-1.8[32], suggesting an uncertainty within ±2.3% for the wind speed above 3 m/s according to IEC 61400-12-1 [33].In [34], a total uncertainty of 3.6% for the wind speed can be assessed as the root sum of squared uncertainties, due to the sensor accuracy (2.3%), the sensor calibration (2%) [35], and the mounting of anemometers on the mast (2%) [36].For the sonic anemometers, the 10-min averaged wind speeds recorded by cup anemometers and sonic anemometers at the same heights agree well with slopes close to 1 and correlation coefficients above 0.993 (Figure 4).Note that the field experiment was carried out in a wind farm in complex terrain.Wind turbine wakes and complex terrain could influence estimating similarity functions.As the wake structure is nonlinear and aligned with wind direction, Han et al. selected the dataset in wind direction sectors with high wind speed correlation to eliminate wake effects for the same measurement campaign [2].Based on their work, we choose three wind direction intervals, (1) the north section: WD = 0 • ± 10 • ; (2) the east section: WD = 60 • ± 5 • ; and (3) the west section: WD = 250 • ± 40 • , to study similarity functions where WD stands for the wind direction, and the intervals are kept narrow, if possible, to ensure data quality with more than 400 points of data.According to Figure 2, the inflow wind in the east and west sections are dominated by complex terrain while the inflow from the north section is dominated by a small slope terrain (slope ≈ 2.4 • ).Therefore, the three wind direction intervals are also applied to study the influence of terrain on the estimation of similarity functions.

Data Processing
The heat flux w θ , the Obukhov length L, and the friction speed u * are evaluated based on measurements of the sonic anemometers [2].In reference to Högström [37], we enforce the criteria (1) u * ≥ 0.1 m/s for the φ m −analysis and (2) |H S | ≥ 10 W•m −2 and u * ≥ 0.1 m/s for the φ h −analysis, where H S = ρc p w θ is the flux of sensible heat and c p = 1004 J•(K•kg) −1 is the constant pressure specific heat of air.The local gradients could be derived using the profile-fitting methods of Nieuwstadt [38] or Högström [37] with more than four levels available to fit profiles of velocity and temperature.However, wind speeds and temperatures are only measured at three and two levels, respectively.To avoid overfitting profiles, the velocity and temperature gradients are approximated at the geometric mean height of levels z 1 = 30 m and z 2 = 70 m by the finite difference form [39,40] ∂Φ ∂z where Φ ∈ {U, Θ}.It is easy to prove that this expression is rigorous for logarithmic profiles and has an error within ±3% for log-linear profiles predicted by MOST.

Limitations and Breakdown of MOST
In this section, we show that the classical φ m consistent with MOST is not valid for R f > 0.25, based on the measured data of the sonic anemometer at 30 m a.g.l for WD = 0 • ± 10 • .The flux Richardson number R f is defined as ∂U ∂z (35) Combining the definitions of L and R f , we have which indicates that R f = 0.25 is a straight line in the plane ζ − φ m .Figure 5 shows that the line R f = 0.25 split the measured data into two groups: for R f < 0.25, the measured data have good agreements with the data predicted by MOST under unstable conditions and show a linear prediction of φ m , with β m = 0.35 under stable conditions; for R f > 0.25, the measured data are significantly lower than those predicted by MOST.For wind energy applications, measured data are usually mixed, and engineers commonly focus on wind turbine wakes under stable conditions without considering the limitation of R f .In this situation, MOST could overestimate the velocity gradient and potentially influence wake modeling of a wind turbine.

Assessment of Similarity Functions for Predicting Wind Profiles
Plots of the dimensionless velocity gradients versus the stability parameter in three wind direction intervals are presented in Figure 6.For the three wind direction intervals, the upslope inflows towards the mast are accelerated more at low levels than at high levels [41], and the dimensionless velocity gradient is decreased consequently.For the east and west sections in complex terrain, the fitted φ m based on field measurements is 30% lower than that of the north section, with a small slope terrain for 0 < ζ < 0.1.Similarly, a decrease of 30% in φ m (= 0.7 + 1.5ζ) was also observed in a stable boundary layer over the complex terrain of the Loess Plateau, China [42].The fitted φ m are also compared with the classical similarity function φ m,cls = 1 + β m ζ and the SHEBA profile function φ m,SHEBA = 1 + 6.5ζ(1 + ζ) 1/3 /(1.3 + ζ) [43].As shown in Figure 6, all the fitted φ m and the SHEBA profile function can limit the velocity gradient under strongly stable conditions: φ m ∝ ζ 1/4 (fitted) and The fitted similarity functions are evaluated for predicting wind profiles based on the Boulder Atmospheric Observatory (BAO) dataset [44].Among the seven runs presented in Sorbjan et al. [44], only runs 3, 5, and 7 are studied, due to the criteria |H S | ≥ 10 W•m −2 .The wind profiles predicted by varying similarity functions are presented in Figure 7 and estimated in numerical integration by where the aerodynamic rounghess z 0 is ~0.01 m for the BAO dataset [44].According to Figure 7, the fitted φ m of the north section performs better in predicting wind profiles for runs 3, 5, and 6, whereas φ m,cls and φ m,SHEBA overestimate the velocity gradient significantly.The wind speed is overestimated on average by 4.6 m/s, 3 m/s, and 0.5 m/s for φ m,cls , φ m,SHEBA and the fitted φ m of the north section, respectively.Influenced by complex terrain, the fitted φ m of the east and west sections underpredicts the velocity of 1.2 m/s averagely.As φ m (0) < 1 of the east and west sections disagrees with the logarithmic law of wind profiles in flat terrain under neutral conditions [45,46], we apply the similarity functions estimated from the dataset of the north section in the wake simulations.

The Proposed Similarity Functions
In the full range of R f , the corresponding proposed similarity functions based on the dataset of the north section under stable conditions are as follows, φ m,exp (ζ) = (1 + 40ζ) 1/4  (38) where χ = 0.9.The dimensionless temperature gradient is observed significantly lower than the classical value for ζ > 0.1 (Figure 8).Grachev et al. [31,43] also observed a similar phenomenon and then proposed a corresponding function:

Test Cases
In this work, under stable conditions, we investigate the wakes of a 500 kW NTK500/41 wind turbine at the Risø Campus test site of DTU in Denmark [1] and three 180 kW Danwin turbines on the island of Gotland in the Baltic Sea [6].The two test cases cover wake measured data 1D to 9D downstream of the turbines where D is the rotor diameter.Details of test cases are listed in Table 2.The first case is based on LiDAR measurements of the wakes of a NTK500/41 turbine.In this campaign, a pulsed LiDAR, which was mounted on a platform at the rear of the nacelle, pointed its laser downstream the turbine up to 5D downstream the turbine.A constant downhill slope of about 0.3% was observed downstream the turbine.Inlet meteorological properties, such as wind speed, wind direction, air temperature, and atmospheric pressure, were measured from a 57-m-tall meteorological mast located upstream the turbine.The NTK500/41 turbine is a stall-regulated 500 kW wind turbine equipped with LM 19.1 m blades and its rotor speed is fixed at 27.1 rpm.Its blades consist of two kinds of airfoils: FFA-W3-XX1 (from 16 to 50% span) and NACA 63-XXX (from 60 to 100% span) [47].Chord length and twist angle over the blades were presented in Johansen et al. [48].The power curve and the thrust coefficient curve based on BEM calculations are shown in Figure 9.The power curve from BEM calculations is shown to have good agreement with the RANS computations in the full turbulence model [48] below 15 m/s and measurements [48] below 10 m/s (Figure 9).The measured thrust coefficient from strain gauges measurements [49] is also shown to have good agreements with the BEM computations above 6 m/s.In this case, we set z 0 = 0.095 m, L = 29 m, and the hub-height reference wind speed U hub = 6.76 m/s, as presented in Machefaux et al. [1], and apply C T = 0.83 based on the BEM computations.The second case is based on wake measurements of three Danwin wind turbines using two 54 m meteorological masts, M1 and M2.In the campaign, the masts were equipped with wind sensors at 8 levels from 10 m to 53.3 m.Inlet conditions including temperature profiles were measured on mast M1, whereas wake profiles of the three turbines were measured on mast M2.The distances from the three turbines to mast M2 are 4.2D, 6.1D, and 9.6D, respectively.In this case, we set the roughness length, the Obukhov length, the reference wind speed and the thrust coefficient to be 0.0005 m, 35 m, 8 m/s and 0.82, respectively [6,50].The detail information of test cases is shown in Table 3 where TI= 2  3 k hub /U hub is the ambient turbulence intensity at hub height and k hub is the corresponding hub-height turbulence kinetic intensity [16].The measured ambient turbulence intensity in the first case is ~10%.

Computational Domain and Meshing
The computational domain has a length of 20D, a width of 10D, and a height of 10D (Figure 10).The background mesh, whose refinement level is 0, consists of 100 × 60 × 60 in length, width, and height.The vertical grids are clustered near the ground and the first cell height above the ground is set to 7.38 z 0 to produce appropriate kinetic energy [51].The mesh around the wind turbine and in the wake region is refined as shown in Figure 10.The mesh is refined around the disk region to ensure about 80 cells through the rotor diameter [52].The complete mesh is composed of ~1.6 million cells.

Boundary Conditions
The boundary conditions consistent with similarity functions are applied to modeling the atmospheric boundary stratification.Wind profile at inlet is estimated by Equation (37).The inlet profiles of potential temperature, TKE, and its dissipation are given by The vertical profiles of potential temperature are estimated in numerical integration.At the outlet, zero gradients of U, Θ, ε, k are applied.For the top boundary, the upstream flow properties are maintained constant.The left and right sides of the computational domain are set to be symmetry.
The near-wall treatment of Temel et al. [53] is implemented in the near ground region to calculate turbulent dissipation rate and the TKE production G k,p due to shear (P p ) and buoyancy (B p ): where the subscript p denotes the first cell center above the ground, the equivalent friction wind speed is u * k = C 0.25 µ φ −0.5 k,p k 0.5 p , and the wall shear stress is τ w = µ t,p (dU/dz).The eddy viscosity µ t,p and the turbulent heat conductivity are imposed at the first cell center above the ground as [54] where the laminar viscosity µ is not ignored to keep consistent with Equations ( 13) and ( 14).

Implementations of Wake Models
In this paper, wake models are implemented in a CFD framework Open-Source Field Operation And Manipulation (OpenFOAM) [55] for developing customized numerical solvers and pre-/postprocessing utilities for fluid flow problems using the finite-volume method.It allows developing user-defined boundary conditions, turbulence models, and solvers, or adding source terms in a convenient way by extending or modifying existing models.As shown in Figure 11, we implement all the test turbulence models (the proposed, Laan, El-Askary, and AM) and develop a large time-step transient solver using the PIMPLE (merged PISO-SIMPLE) algorithm to simulate wakes of wind turbines under stable conditions.At the beginning of simulations; this solver first initializes the flow field based on boundary conditions.During this period, similarity functions and their derivatives are also evaluated and then used in turbulence models to estimate S k,ABL and C ε3 if required.This solver then solves the continuity, momentum and energy equations based on the eddy viscosity µ provided by the turbulence models before in iterations the finish of the simulation.During the iteration, source terms S u and S ε,wake are modeled based on local flow information via user-specified finite volume options (fvOptions).In simulations, we apply a second-order backward difference scheme for temporal discretization and a second-order central difference scheme for spatial discretization of other derivative terms.

fvOptions
Estimating S u,i and S e,wake

Input: Wind Turbine Information
AD-CT: Thrust coefficient curve AD-BEM: radius, chord, twist, etc S e,wake

Input: Inflow conditions
Similarity functions f m,cls or f m,exp ?

Results of Wake Simulations
6.1.Case 1: Wakes of a NTK500/41 Wind Turbine In Machefaux et al. [1], the wakes of a NTK500/41 wind turbine was studied experimentally and numerically using two LES methods: a classical method based on the ELLIPSYS3D flow solver [56] and an extended approach that explicitly introduces thermal and Coriolis effects as external force terms into the solver.As the wind profile predicted by MOST was observed to be much more severely sheared as compared with measurement, the classical approach applied a power law for the inlet velocity profile.The extended LES approach additionally carried out a transient precursor computation to simulate the time-varying vertical structure of the whole ABL, and it used the results of the precursor simulations to impose the mean potential temperature and velocity profiles at the inlet for wake modeling.The precursor simulations are shown to have effects on limiting the overestimated velocity gradient predicted by MOST [1].Numerical simulations have shown that there are no major differences in wake deficit predictions between the classical approach and the extended model for stable and unstable cases.Therefore, we only compare the wake deficit predicted by the extended approach with the test models in RANS in this work.
Figure 12 shows contours of wake deficits, ∆U/U 0 , predicted by different models at hub height under stable conditions, where ∆U is the velocity deficit in the wake and U 0 is the inlet velocity.Simulations using the LES technology show a faster recovery wake than that in the proposed model, whereas the Laan model predicts a significant slower recovery wake with a larger deficit.One possible explanation for the stronger wake effects predicted by the Laan model is that they are due to the lower turbulence intensity predicted by the classical similarity functions.For a given reference wind speed U hub and a fixed roughness length z 0 , a higher velocity gradient in the Laan model results in a smaller friction speed u * and a lower turbulence intensity I hub than the proposed model, and thus it weakens the wake recovery.The ambient turbulence intensities for the LES approach of 9.7% and the proposed model of 10.2% are close to the measured value of 10% while the value for the Laan model is reduced by 40% as compared to the measured data.As compared with measurements, the proposed model shows better performance than other approaches based on MOST combined with RANS for the longitudinal distance of 1-5D and the LES approach for the longitudinal distance below 3D (Figures 13 and 14).The AM, El-Askary, and Laan models overestimate the wake deficits in both vertical and lateral directions for longitudinal distances above 2D.All the test models in RANS fail to predict wake deficits at 4-5D downstream of the rotor.This disagreement is probably due to a combination of terrain effects and experimental uncertainties using LiDAR.According to the Space Shuttle Topography Mission (STRM)-based terrain data, the test site terrain has a downhill slope characterized by a height difference of 5.5 m from the rotor location to the most downstream cross section.As the wake has been observed to follow terrain under stable conditions [1,8], the wind speed for the longitudinal distance above 5D at the hub height is actually the value of 5.5 m above the wake center, which results in the overestimation of wake by various turbulence models in RANS.∆U is the velocity deficit in the wake and U 0 is the inlet velocity.
For the vertical profiles of wake deficits, Figure 14 shows that the proposed model overestimates wake deficits below the hub height at 5D.A double-bell near-wake shape, due to a lower energy extraction around the blade root [1], is clearly observed under stable conditions at 1D and is captured by the actuator models based on BEM calculations.There are no major differences between the vertical wake profiles predicted by the proposed models based on either the thrust coefficient (Proposed-CT) or BEM calculations (Proposed) above 2D.As the LES model applies slip wall conditions at the bottom of the computational domain and does not consider the roughness effect of the ground, large negative wake deficits are observed near the ground [1].

Case 2: Wakes of Danwin 180 kW Wind Turbines
The numerical results of the turbulence models based on MOST and the proposed model based on the new similarity functions are compared with the experimental data reported by Magnusson et al. [6] in  A comparison is made at distances 4.2, 6.1, and 9.6D downstream of the turbine.As the detailed rotor geometry of the Danwin 180 kW wind turbine is not available, the actuator disk model based on the thrust coefficient is applied.the change of the wake deficit with the distance downstream of the turbine in both lateral and vertical directions.From this figure, it can be depicted that there is no discernible difference between AM, El-Askary, and Laan models at all downstream positions.As compared to measurements, these models that are based on MOST, predict larger wake deficits, and the wake deficit can be overestimated by 0.25 at 6.1D downstream of the rotor.Introducing the new similarity functions into the Laan model improves the wake prediction at 6.1D and 9.6D downstream, as compared with the models based on MOST.Note that the measured wake center or the position of the maximum wake deficit, shifts ~0.2D in both vertical and lateral directions at 6.1D and 9.6D downstream of the rotor.This shift of the wake center was also observed in full-scale turbine wakes with scanning-LiDAR measurements from the Crop Wind Energy eXperiment (CWEX) in Iowa [57].The stretching of the wake structures can be attributed to the strong wind veer associated with stable conditions [6,57].As wind veer is a direct result of Coriolis effects caused by the Earth's rotation and the Coriolis force is not modeled in this work, all the test models fail to capture this stretching of the wake structures.
Figure 17 shows the distribution of the normalized added turbulence σ 2 u /σ 2 u0 − 1 in the lateral direction at hub height, where σ u and σ u0 is the standard derivations of wind speed in wakes and in the atmosphere, respectively.A dual-peak pattern (lateral distance y − y c = ±0.5D) is detected experimentally at 4D downstream of the rotor, as a result of rotor tip vortices and high shear production of turbulent kinetic energy caused by strong velocity gradients at the wake boundary [3,4].All test models capture this dual-peak pattern with lower values: the measured peaks reduce by 3 in the simulations using models based on MOST and reduce by 5 in the simulations using the proposed model.This underestimation of the added turbulence may be due to representing the turbine rotor with actuator disk instead of the actual rotor geometry.Note that the dual-peak pattern predicted by the models based on MOST can be observed at distance downstream of the rotor up to 9.6D, which becomes indistinctly at 6.1D and 9.6D downstream of the rotor for both measurements and the simulations using the proposed model.As compared to measurements, the models based on MOST also overestimate the normalized added turbulence which is improved by the proposed model at 6.1D and 9.6D downstream of the rotor.

Model Assessment
In this section, we adopt the root-mean-square error (RSME) to analyze the accuracy of the CFD predictions using various turbulence models [58]: where n is the number of measurement points in the evaluation, y EXP is the experimental data, and y CFD is the simulated values from CFD.Note that measurement points out of the wake range are not applied to calculate the RSME, and the wake width is set 2D due to wake simulations.The RMSE analysis of the wake deficit is presented in Figures 18 and 19.Over a range of 0.04 to 0.2, the lower the RSME value, the better the behavior of the numerical model.According to Figure 18, it can be depicted that the proposed model, either based on BEM in Case 1 or based on the thrust coefficient in Case 2, has the best performance with the lowest RMSE values, except for 4-5D in Case 1.The AM, El-Askary, and Laan models have similar performance in the wake prediction, but the El-Askary model has lower RMSE values among the various turbulence models based on MOST.As compared to AM and Laan model, the El-Askary model could reduce the RSME of the lateral wake deficit by 0.03 in Case 1 and 0.01 in Case 2. However, this improvement is smaller than that of the proposed model: introducing the new similarity functions into the Laan model could reduce the RSME of the lateral wake deficit by 0.05 averagely (Figure 18).The LES approach presented in Machefaux et al. [1] predicts a decreased RSME of the lateral wake deficit with the distance downstream of the rotor.Note that all test models, except for the LES approach, have poor performance in predicting wake structure at 5D downstream of the rotor.This could be a result of the downslope terrain or uncertainties due to LiDAR measurements [1].
Using BEM to distribute force through the rotor significantly improves near wake prediction at 1-3D downstream of the rotor in Case 1: the RSME of vertical wake deficit predicted by the proposed model using BEM approach is lower, 0.02-0.05,than the proposed model using the C T approach.However, it should be noted that the proposed model using C T approach has a similar performance in the wake prediction with that using BEM approach at 4D-5D downstream of the rotor.This suggests that the distribution of the momentum source through the rotor only affects near wake structures up to 3D and that one can expect similar accuracy of wake prediction for far wake regardless of which method to distribute the force through the rotor.
The RMSE analysis of the added turbulence is presented in Figure 20.The RMSE of the added turbulence is observed to decrease with the distance downstream the rotor.Turbulence intensities are significantly underpredicted at 4.2D, where the RMSE of the added turbulence is ~4.More than 4.2D the RMSE of the added turbulence predicted by the proposed model is lower by 30% than those of models based on MOST.This indicates that introducing the new similarity functions into the Laan model can eliminate the side effects of the breakdown of MOST on wake modeling.

Conclusions
In this paper, the breakdown of MOST is investigated experimentally, and its side effects on wake modeling through numerical simulations are also investigated.New similarity functions based on a field measurement in a wind farm are proposed and applied in a turbulence model to eliminate the side effects of the breakdown of MOST on wake modeling.Wake simulations of two types of turbines under different stability conditions are carried out and compared with measurements from a LiDAR and cup anemometers.The main findings are as follows.(1) Field measurements show that MOST overestimates velocity gradient for ζ > 0.1 in the full range of R f .The proposed similarity functions can limit the velocity gradient as the stability increases.(2) Due to the breakdown of MOST for ζ > 0.1, the test models based on MOST overestimate wake effects in both wake deficits and the added turbulence under stable conditions.(3) The new similarity functions constrain velocity gradient for ζ > 0.1 as compared with MOST, and introducing the new similarity functions into the Laan model improves the wake prediction under stable conditions.(4) By distributing the blade force through the rotor, momentum effects of the rotor to the atmospheric boundary layer are simulated in a more detail way than the uniformly distributed blade force applied based on the thrust coefficient.This enables the proposed model to capture the double-bell near-wake shape.

Figure 2 .
Figure 2. Complex terrain around M1, where Wind direction (WD) stands for the wind direction and D = 93 m is the diameter of wind turbine 14.

Figure 4 .
Figure 4. Comparison of the wind speeds from sonic anemometers (U sonic ) and cup anemometers (U cup ) for the wind direction in 0 • ± 10 • : measurements in symbols and fit curves in solid lines.r xy is the Pearson correlation coefficient.

mFigure 6 .
Figure 6.Plots of φ m of different models as functions of ζ: dots (measured data) and circles (bin-averaged measured data).n is the number of the measured data.

Figure 7 .
Figure 7. Velocity profiles estimated based on different similarity functions φ m (ζ) (the upper part) and boxplots of errors of velocity predictions (the lower part).U exp is the measured velocity in runs 3, 5, and 6, whereas U predict is the value predicted at the corresponding height.

Figure 8 .
Figure 8. Same as Figure 6 but for φ h of the north section.

Figure 12 .
Figure 12.Normalized velocity distribution at hub height under the stable condition, where large eddy simulation (LES) data are modified from Machefaux et al.[1], where (x c , y c , z c ) is the center point of the disk, ∆U is the velocity deficit in the wake, and U 0 is the inlet velocity.

Figure 13 .
Figure13.Lateral wake deficit of a NTK500/41 wind turbine at hub height under the stable condition; ∆U is the velocity deficit in the wake and U 0 is the inlet velocity.

Figure 15 .
Figure 15.Same as Figure 13 but for Danwin 180 kW wind turbines.All turbulence models are combined with actuator disk model (AD) based on thrust coefficient.

Figure 16 .
Figure 16.Same as Figure 15 but for vertical wake deficit profiles.

Figures 15 and 16
Figures 15 and 16  demonstrate the change of the wake deficit with the distance downstream of the turbine in both lateral and vertical directions.From this figure, it can be depicted that there is no discernible difference between AM, El-Askary, and Laan models at all downstream positions.As compared to measurements, these models that are based on MOST, predict larger wake deficits, and the wake deficit can be overestimated by 0.25 at 6.1D downstream of the rotor.Introducing the new similarity functions into the Laan model improves the wake prediction at 6.1D and 9.6D downstream, as compared with the models based on MOST.Note that the measured wake center or the position Figures 15 and 16  demonstrate the change of the wake deficit with the distance downstream of the turbine in both lateral and vertical directions.From this figure, it can be depicted that there is no discernible difference between AM, El-Askary, and Laan models at all downstream positions.As compared to measurements, these models that are based on MOST, predict larger wake deficits, and the wake deficit can be overestimated by 0.25 at 6.1D downstream of the rotor.Introducing the new similarity functions into the Laan model improves the wake prediction at 6.1D and 9.6D downstream, as compared with the models based on MOST.Note that the measured wake center or the position

Figure 17 .
Figure 17.Lateral added turbulence in wakes of Danwin 180 kW wind turbines at hub height under the stable condition; σ u and σ u0 represent the standard derivations of wind speed in wakes and in the atmosphere, respectively.All turbulence models are combined with AD based on thrust coefficient.

Figure 18 .
Figure18.RSME of the wake deficit at hub height under the stable condition: "-CT" represents models using AD based on the thrust coefficient.

Figure 19 .
Figure 19.Same as Figure 18 but for boxplots of the root-mean-squared error (RSME) of wake deficits.

Figure 20 .
Figure 20.RSME of the lateral added turbulence in wakes of Danwin 180 kW wind turbines under the stable condition at hub height; σ u and σ u0 represent the standard derivations of wind speed in wakes and in the atmosphere, respectively.All turbulence models are combined with AD based on thrust coefficient.

Table 2 .
Details of measurements in the two test cases.

Table 3 .
Test cases and corresponding parameters under stable conditions.