Data-driven RANS for simulations of large wind farms

In the wind energy industry there is a growing need for real-time predictions of wind turbine wake flows in order to optimize power plant control and inhibit detrimental wake interactions. To this aim, a data-driven RANS approach is proposed in order to achieve very low computational costs and adequate accuracy through the data assimilation procedure. The RANS simulations are implemented with a classical Boussinesq hypothesis and a mixing length turbulence closure model, which is calibrated through the available data. High-fidelity LES simulations of a utility-scale wind turbine operating with different tip speed ratios are used as database. It is shown that the mixing length model for the RANS simulations can be calibrated accurately through the Reynolds stress of the axial and radial velocity components, and the gradient of the axial velocity in the radial direction. It is found that the mixing length is roughly invariant in the very near wake, then it increases linearly with the downstream distance in the diffusive region. The variation rate of the mixing length in the downstream direction is proposed as a criterion to detect the transition between near wake and transition region of a wind turbine wake. Finally, RANS simulations were performed with the calibrated mixing length model, and a good agreement with the LES simulations is observed.


Introduction
As recently published on the Wall Street Journal (http://blogs.wsj.com/numbers/what-is-themost-efficient-source-of-electricity-1754/), wind energy is the most efficient source of electricity on a US national average. However, in 2013 the US Department of Energy estimated that typical power losses for a wind power plant are about 20% of its annual production (http://energy.gov/eere/wind/atmosphere-electrons), which are mainly due to wind turbine wake effects, such as chaotic wake interactions and shadowing due to upstream wind turbines. Moreover, wake-related phenomena affect not only power production, but also the overall lifecycle of wind turbines, thus their reliability. Indeed, wind turbines operating in the wake produced by upstream turbines experience highly heterogeneous, unsteady, and turbulent wind flows, which generate larger fatigue loads compared to the typical loads connected to the atmospheric wind.
Today there is an increasing need for tools for wind turbine wake prediction in real-time, which can be executed within a relatively short time frame and coupled with wind turbine control systems in order to maximize power production from a power plant and increase durability of the turbines.
Computational fluid dynamics (CFD) techniques with different degrees of complexity and accuracy have been applied to simulate wind turbine wake flows. More recently, large eddy simulation (LES) frameworks have been developed for wind energy applications. LES computes explicitly large energy-containing non-universal scales of the flow, while the small scales are modeled [1,2,3,4,5]. The wind induced loads, namely thrust and torque, are parameterized using different types of wind turbine models [6], such as the actuator disc model, which distributes uniformly the loads over the rotor disk, and the actuator line model, which distributes the forces along a lifting line that rotates mimicking the blade motion. Although accuracy on the evaluation of wind turbine wake flows has been improved over the last decade, LES still remains far too costly, which makes it inapplicable for operations.
In practice simpler and less accurate analytical models are extensively uses for wind farm siting and prediction of wind turbine wakes. Indeed, even though their fidelity is very limited, these built-to-purpose wind resource methodologies can easily be implemented thanks to their low computational cost [7], such as the commercial software WAsP [8], WindPRO [9], WindSim [10], WindFarmer [11] and OpenWind [12]. One of the pioneering analytical wake model, which is available since the late 1980s, is the one proposed by Jensen [13], which assumes a top-hat shape for the wake velocity deficit. This simple wake model is obtained by only applying the mass conservation. Similar models are derived from the mass and momentum conservation, such as for the Frandsen model [14]. More recently an improved analytical wake model has been proposed in [15].
Summarizing, LES simulations of wind turbine wakes can achieve very accurate results; however, they still imply high computational costs, which makes them not suitable for real wind turbine operations. On the other hand, analytical models have very low computational costs, but they can be very inaccurate due to their simplistic assumptions. In addition, both LES and analytical models may only apply in the diffusive region of the wake, i.e. downstream to the rollup and formation of the wake vorticity structures.
Other CFD tools have been developed for the wind energy field by adopting the Reynolds-Averaged Navier-Stokes (RANS) approach, in which the mean flow is computed and the effects of turbulence are parameterized using classical closure models [16,17,18,19,20]. The RANS approach present the undeniable advantage of requiring very low computational costs compared to LES simulations, thus to be suitable for real-time applications. In [21], the RANS approach is used with a semi-empirical eddy-viscosity model. In that work, the eddy-viscosity model is formulated as a function of the wake velocity deficit evaluated at the downstream location of two rotor diameters, d, and of several empirical parameters representing the effects of wind shear, turbulence intensity, and atmospheric thermal stability. Moreover, a filter function is applied in order to take the turbulence mechanical production occurring in the near wake into account. This filter function is empirically limited for downstream locations with d < 5.5. The main limitation of this approach clearly consists in the extensive use of empirical parameters, which are difficult to tune and, more importantly, to dynamically update for dynamic and operational conditions. The aim of this paper is to propose a CFD tool that requires low computational costs, as typical for a 2D RANS algorithm, in order to provide real-time predictions of wind turbine wakes. At this same time, this CFD tool ensures an adequate accuracy for real-time wind turbine operations, which is achieved the data-driven calibration. This numerical tool is mainly meant for a framework in which a measuring system monitors wind turbine operations, such as a wind LiDAR network [22,23,24,25], and the available data are embedded in a numerical algorithm in order to improve accuracy of the numerical model. Therefore, this approach can be classified as a data-driven algorithm.
The low computational cost is achieved by using a RANS approach , thus by predicting the time-averaged velocity field and by modeling the wind turbulence effects through turbulence closure models [26]. The accuracy of the wake predictions is mainly dependent on the turbulence model, which is dynamically tuned through data assimilation. Specifically, a mixing length model is calibrated by measuring the Reynolds stress between the axial and radial velocity component, and the gradient of the axial velocity in the radial direction. Since the turbulence closure model is calibrated, by measuring a velocity profile at a certain downstream location, it is possible to predict accurately the downstream evolution of the wake velocity field, thus the wind field impacting on downstream waked turbines. The paper is organized as follows. The LES database used for the present work is described in section 2, while the the RANS equations and the formulation of the mixing length model are briefly summarized in section 3. The statistical characterization of the turbulent wake flows obtained for a turbine operating at different tip speed ratios follows in section 4. In this section, the calibration of the mixing length model and the RANS simulations are also presented. Finally, the conclusions are given in section 5.

Dataset obtained from LES simulations
In this work, high fidelity LES simulations of a utility scale wind turbine are used for the tuning of the turbulence closure model adopted for the RANS simulations. The LES simulations were performed by the CFD Lab at UT Dallas with the in-house developed code UTD Wind Farm (UTDWF). The flow around the wind turbine blades is modeled using the actuator line model. The lift and drag coefficients are obtained from documented airfoil data using the calculated local angle of attack. The latter is corrected taking into account the tip and root vortex by means of a modified Prandtl correction factor. The applied aerodynamic forces are distributed smoothly among the points corresponding to the section of the airfoil in order to avoid a singular point force that may causes numerical instabilities. The forces are distributed in a two-dimensional manner by applying a regularization Gaussian distribution kernel. The tower and nacelle are treated by the efficient immersed boundary technique described in detail by [27]. This method allows the flexibility of imposing complex geometries inside the domain without the need of body fitted grids, thus reducing the computational time.
The code UTDWF has been validated with experiments performed at Norwegian University of Science and Technology using a wind turbine placed inside a wind tunnel. The computational box is 9d × 10d × 3d (d is the rotor diameter of 126 m) in the streamwise, wall normal and spanwise direction, respectively. The grid is 512 × 256 × 384 in streamwise, spanwise and wall normal direction, respectively. The grid is stretched in vertical direction, and it is uniform at the rotor heights. Along the blade span, 40 grid points are present.
A uniform velocity profile is used as inflow condition. Radiative boundary conditions are imposed at the outlet, a no slip condition is imposed on the lower boundary (the ground) and on the tower and nacelle, while free slip is assumed at the upper boundary. Periodicity is imposed in the spanwise direction to mimic infinite columns of turbines.

RANS formulation and mixing length turbulence model
Data-driven and low computational cost simulations of the wind turbine wake were performed with the Reynolds-Averaged Navier-Stokes (RANS) equations formulated in cylindrical coordinates. The equations are made non-dimensional by taking as reference length the rotor diameter d and the oncoming velocity at hub-height U hub . The pressure is made non-dimensional by the fluid density and U hub . Consequently, the Reynolds number is defined as Re = U hub d/ν. The RANS equations are the following: (1) where U are p are the time-averaged velocity vector and pressure, respectively. The timederivatives in the Navier-Stokes equations are omitted because only steady state is simulated. The last term in the right-hand-side of the momentum equation represents the Reynolds stresses, and u is the turbulent fluctuating velocity vector. The operator (·) denotes time averaging.
In order to solve the RANS equations, the Reynolds stresses term has to be modeled [26]. Through the Boussinesq hypothesis, the deviatoric part of the Reynolds stress tensor is expressed as a function of the strain rate tensor S = 1 2 [∇ + ∇ T ]U, while the isotropic part is absorbed in the the pressure term: where q is the turbulent kinetic energy (TKE) and I is the 3x3 identity matrix. Each component of the Reynolds stress tensor evaluated through the Boussinesq hypothesis is reported in the following: The turbulent eddy-viscosity, ν t , is modeled by using the generalized mixing-length algebraic model suitable for shear-flows [26]. Specifically ν t is expressed as a function of the strain rate tensor and the mixing-length l m , as where l m is only a function of the streamwise position. By considering an axisymmetric wake flow, the term (2S : S) 1/2 reduces to: Therefore, the full set of RANS equations for an axisymmetric flow results: where the advection operator Γ is At the inlet a specific velocity profile is given. Free-stresses conditions are imposed at the outlet and at the lateral boundary, while symmetry conditions are set on the axis. Eq. 6 together with the boundary conditions are discretized using a staggered pseudospectral Chebyshev-Chebyshev collocation method implemented in Matlab. The three velocity components are collocated at the Gauss-Lobatto-Chebyshev (GLC) nodes, whereas the pressure is defined on the Gauss-Chebyshev nodes (GC). Specifically, the momentum equation is collocated at the GLC nodes, and the pressure is interpolated from GC points to GLC points. Conversely the continuity equation is defined on the GC grid and the velocity components are interpolated from the GLC grid. Consequently, the two grids are mapped in the physical domain through a linear mapping.
Being the time-averaged Navier-Stokes equations steady, the solution is obtained by Newton-Raphson method. The iterative procedure is arrested when the L 2 -norm of the difference between two consecutive solutions is less than 10 −8 .

Turbulence statistics of the wind turbine wake and data-driven RANS
The turbulent wake velocity field produced by the wind turbine is investigated in order to perform an efficient calibration of the RANS turblence closure model from the available data. To this aim, LES simulations are considered, which were performed for a turbine operating with different tip speed ratios (T SR), i.e. the ratio between the rotational velocity at the tip of the turbine blade and the velocity of the incoming wind at hub height. In Fig. 1, the time-average of the axial velocity evaluated over the horizontal plane at hub height is shown for the wind turbine operating with T SR equal to 5.8, 6.5, 7 and 7.5. The cross-width of the wake seems to be barely affected by the the T SR variation; however, an enhanced velocity deficit is observed by increasing T SR, which is related to a larger power harvesting. Indeed, T SR=7.5 is to a good approximation the optimal set-point for the tested turbine.
The different terms in Eq. 5 are now investigated in order to evaluate their contribution to the eddy-viscosity formulated in Eq. 4. In Fig. 2, the results related to the simulation performed with T RS=5.8 are presented. It is evident that the main contribution to (2S : S) 1/2 ( Fig. 2(a)) is represented by the term (∂U x /∂r + ∂U r /∂x) reported in Fig. 2(f). Moreover, within this term the most significant contribution is provided only by the term (∂U x /∂r) in Fig. 2(g), which is the gradient of the axial velocity in the radial direction. By considering the Boussinesq hypothesis used for modeling the Reynolds stress tensor, in Eq. 3 the only component affected by (∂U x /∂r) is u x u r . The latter represents the turbulent interaction between the wake flow and the surrounding unperturbed wind. Therefore, the higher u x u r , more flow entrainment is promoted while the wake moves downstream. In Fig.  3, the Reynolds stress component u x u r is evaluated over the horizontal plane at hub height for the turbine operating with different TSR values. Generally, it increases by proceeding downstream with the maximum intensity located at a radial distance roughly equal to the rotor radius. Moreover, the region characterized by a significant value of u x u r becomes wider by moving downstream, which represents the enhanced diffusive phenomena in the far wake. With increasing TSR, the intensity of u x u r becomes more significant, which is connected with the larger mechanically produced turbulence, thus to the higher gradient of the axial velocity in the radial direction, (∂U x /∂r). Therefore, it is concluded that, up to the optimal T SR the larger the T SR, more power is generated, thus producing a larger velocity deficit. Consequently, the larger gradient of the axial velocity in the radial direction leads to an enhanced mechanically produced turbulence, which promotes flow entrainment and wake recovery.
Accordingly to the discussion above, it derives that the most meaningful procedure to evaluate the mixing length through Eqs. 3 and 4, is to fit the Reynolds stress component u x u r through  the time-averaged velocity gradient (∂U x /∂r) as follows: The results of the best fit procedure are reported in Fig. 4 for the case with T SR=5.8. In Fig.  4(a), the absolute value of the Reynolds stress component u x u r is compared to l 2 m (∂U x /∂r) 2 in 4(b), where the mixing length is obtained at each downstream location through a best fit procedure.
The mixing length evaluated for the different T SR values over the horizontal plane at hub height and over the vertical symmetry plane of the wind turbine wake is reported in Fig. 5. For all the analyzed cases, the mixing length is roughly invariant for different downstream locations in the very near wake, i.e. for x/d < 1, then it increases for larger downstream distances. This wake feature can be very useful for a sharp detection of the near wake from the transition region, denoted also as diffusive region [28]. Indeed, in the near wake the rollup of the wake vorticity structures and the wake formation occur, whereas further downstream diffusive phenomena dominate leading to the wake recovery.
It is also interesting to observe that at the beginning of the diffusive region, i.e. at x/d ≈ 1, the mixing length has roughly the same value for the different T SR values, which is expectable due to the same incoming wind condition simulated for the different cases. Another general feature is the roughly linear increase of the mixing length with the downstream distance, which is in agreement with previous wind tunnel experiments presented in [29,30]. Moreover, the slope of the mixing length as a function of the downstream location increases with increasing T SR   value. Therefore, wind turbines operating with larger T SR capture more power and produce a larger velocity deficit, but their wakes are characterized by a higher recovery rate. Indeed, the enhanced diffusive phenomena are related to the enhanced mechanically produced turbulence that, in turn, are related to the high axial velocity shear in the radial direction.
Since the mixing length model is calibrated through the available data as explained in the previous section, RANS simulations for the evaluation of the downstream evolution of the wind turbine wake flow can be easily executed within a very limited amount of time. For this work, the velocity profile at a downstream location of x/d = 2 is provided as boundary condition at the inlet of the computational domain, and the downstream evolution is calculated till a downstream distance of x/d = 5. The velocity field over the horizontal plane at hub height is shown in Fig. 6 for different T SR values. A general good agreement is observed with respect to the time-averaged velocity field obtained through the LES simulations reported in Fig. 1.

Final discussion and conclusions
In this paper a data-driven RANS framework for low computational cost and accurate simulations of wind turbine wake flows is proposed. The limited computational requirements are obtained through the RANS formulation, and an adequate accuracy is achieved through a data assimilation procedure for the calibration of the turbulence closure model. Specifically, a mixing length model for swirling flow is adopted for which the eddy-viscosity is only function of the mixing length and the gradient of the axial velocity along the radial direction. Moreover, the mixing length is only a function of the downstream location. High fidelity LES simulations of a wind turbine operating with different tip speed ratios are used as test case for the present work. It has been shown that the main contribution of the time-averaged velocity field to the evaluation of the eddy-viscosity through the Boussinesq hypothesis is represented by the gradient of the axial velocity component along the radial direction. Consequently, the most meaningful turbulence parameter for the evaluation of the mixing length results to be the Reynolds stress component between the axial velocity and the radial velocity. Indeed, the mixing length is evaluated at each downstream location as best fit between the above-mentioned Reynolds stress component and the gradient of the axial velocity along the radial direction.
It has been found that the mixing length is roughly invariant in the very near wake, and its value is not affected by the operational conditions of the wind turbine. This is an important feature that allows in the near wake to evaluate the mixing length only through a characterization of the incoming wind field, and by neglecting the wind turbine characteristics. At a fixed location, which is one rotor diameter downstream for this work, the mixing length starts increasing linearly with the downstream location. This feature, i.e. the slope change of the mixing length as a function of the downstream location, is considered an effective criterion to sharply detect the transition from the near wake region and the transition region. The former is the wake region where the rollup of the wake vorticity structure occurs, whereas the latter is the area where the diffusive phenomena dominate and the wake gradually recovers to the incoming wind conditions. In the transition region, the slope of the mixing length in the downstream direction raises with increasing tip speed ratio. This suggests that a turbine operating with a higher tip speed ratio is characterized by a faster wake recovery. This feature is a consequence of the enhanced flow entrainment and mechanically produced turbulence generated by the higher shear of the axial velocity in the radial direction. Since the mixing length is calibrated through available wake measurements, the downstream evolution of a wind turbine wake can be quickly evaluated through RANS simulations by providing the velocity profile at a certain downstream location.
This data-driven RANS approach will be tested to simulate multiple turbine wakes, where the near wake velocity profile can be predicted through different model, such as actuator line or actuator disc. The main objective of the proposed tool is to perform quick predictions of wind turbine wakes by performing a limited number of wake measurements.