A CFD framework for offshore and onshore wind farm simulation

We present a wind simulation framework for offshore and onshore wind farms. The simulation framework involves an automatic hybrid high-quality mesh generation process, a pre-processing to impose initial and boundary conditions, and a solver for the Reynolds Averaged Navier-Stokes (RANS) equations with two different turbulence models, a modified standard k-ϵ model and a realizable k-ϵ model in which we included the Coriolis effects. Wind turbines are modeled as actuator discs. The wind farm simulation framework has been implemented in Alya, an in-house High Performance Computing (HPC) multi-physics finite element parallel solver. An application example is shown for an onshore wind farm composed of 165 turbines.


Introduction
Simulation of wind farms with Computational Fluid Dynamics (CFD) models involves the resolution of the turbulent Atmospheric Boundary Layer (ABL) and the effects induced by wind turbines, including wind speed deficit, increase of turbulent kinetic energy, and interaction among wakes. Commonly, wind turbines are modeled as uniformly loaded actuator discs [1], an option which allows a compromise between a reasonable model accuracy and a low computational cost.
However, it is well known that the Reynolds-Averaged Navier-Stokes (RANS) standard k-ε model significantly overestimates Reynolds stresses [2] behind actuator discs, resulting on a significant underestimation of velocities and on an excessive wake damping. Several modifications of the standard k-ε model have been proposed for delaying the simulated wake flow recovery [3], but most of them require the adjustment of model parameters that depend on the characteristics of the wind turbine at stake. Shi et al. [4] proposed a realizable k-ε model that enhances the wake predictions of actuator discs [5] and needs no adjustment of parameters. On the other hand, it is also well-established that the k-ε model needs an additional mixing length limitation model [6] to accurately predict wind intensity profiles and wind veering with height in the ABL when considering Coriolis forces. However, no attempts have been made to modify the realizable model for the simulation of Coriolis forces in the ABL.
Several aspects of wind farm modeling pose requirements on the computational mesh. Firstly, the existence of a boundary layer imposes stretching requirements in order to capture nearsurface sharp gradients. Secondly and for onshore farms, the underlying terrain must be discretized and properly approximated by the mesh. Finally, the actuator discs have to be embedded in the mesh with a subsequent mesh refinement around and downstream the turbines in order to capture near wake flow effects properly. This is a constrain for conformal structured meshes, which do not allow to prescribe finer resolutions around turbines without extending the finer discretization across all directions of the computational domain, thereby increasing dramatically and unnecessarily the number of computational cells (nodes). Several codes circumvent this problem by loosing mesh conformity (e.g. use of hexahedral cells/elements with hanging nodes). In contrast, we propose an alternative based on hybrid element meshes, which allow local mesh refinements while keeping mesh conformity. The objective is to present a High Performance Computing (HPC) simulation framework for offshore and onshore wind farms with an automatic high-quality mesh generation process.
The rest of the manuscript is arranged as follows. Section 2 presents the two RANS physical models; a standard k-ε adapted to ABL flows [6] (i.e. with Coriolis effects, mixing length limitation and appropriate model constants), and the k-ε realizable model of [4] modified to account for Coriolis forces and with appropriate coefficients for the simulation of ABL flows. In both cases, wind turbines are simulated under the actuator disc theory [7]. Section 3 focus on numerical aspects including the hybrid mesh generation, the implementation of the actuator disc model, the numerical method and its implementation in the multi-physics parallel solver (Alya) and, finally, a succinct model validation. Section 4 shows an application example for an onshore wind farm consisting on 165 turbines deployed on complex terrain. The effect of Coriolis terms is analyzed for both (k-ε modified and k-ε realizable) RANS models. Finally, Section 5 wraps-up and briefly discusses model limitations and concluding remarks.
2. RANS models governing equations 2.1. k-ε modified model Considering the flow as incompressible and isothermal (neutral stability), the modified k-ε RANS model accounting for Coriolis effects and using the Apsley and Castro correction for the mixing length limitation [6] are written as: where the unknowns are the velocity field u, pressure p, turbulent kinetic energy k, dissipation rate of turbulent kinetic energy ε, and turbulent viscosity ν t (computed with the diagnostic equation (5)). The fourth term on the left hand side (LHS) of momentum equation (2) models the Coriolis force being ω the Earth's angular velocity. The sixth term on the LHS of equation (2) is the actuator disc force, which is active only inside the disc volume, where C t is the thrust coefficient, U ∞ is the free-stream velocity at hub height, n d the disc normal unit vector (pointing opposite to inflow), and ∆ is the thickness of the disc. The forces inside each disc volume are uniformly distributed. In the turbulence equations (3)-(4), the term P k = ν t S is the kinetic energy production due to shear stress, with S = ∇ s u : ∇ s u (∇ s denotes the symmetrical gradient operator). For the coefficients of the k-ε modified model we follow Panofsky and Dutton [8] and adopt: C µ = 0.0333; C 1 = 1.176; C 2 = 1.92; σ k = 1.0; σ ε = κ 2 where κ is the Von Karman constant. The coefficient C 1 in the RHS of the dissipation equation (4) is a modified coefficient, originally proposed by Apsley and Castro [6], to prevent the increase of mixing length l m = C 3/4 µ k 3/2 /ε above a maximum value l max when accounting for Coriolis effects: where l max is calculated as [9] : being u g the geostrophic wind velocity and λ the latitude. Note that, if no Coriolis forces are considered (i.e. |ω| = 0), then l max → ∞ and C 1 = C 1 .

k-ε realizable model with Coriolis effects
The k-ε realizable model proposed by Shi et al. [4] presents the advantage of satisfying realizability conditions on the Reynolds stresses. This model is known to improve the accuracy of flows involving detachment and re-circulations, and specifically, to enhance the prediction of wakes when using actuator disc models [5]. The realizable model shares the same turbulent kinetic energy equation (3) with the standard k-ε model. However, differences exist for the dissipation rate (4) and turbulent viscosity (5) equations. The turbulent viscosity ν t is calculated from (5) with a variable C µ that depends on the local values of k, ε and velocity gradients [4]. In turn, the dissipation equation is derived from the vorticity fluctuation transport equation, which is modified in the present work to use the model coefficients for ABL flows (6) and corrected with the Apsley and Castro limitation model to account for Coriolis effects in the ABL: where ν is the laminar viscosity and C µ 0 is the value taken by coefficient C µ under homogeneous shear, boundary layer logarithmic profile.

Boundary conditions
Proper boundary conditions need to be added to the Navier Stokes (1)-(2) and turbulence k-ε (3)-(4) ( (3) and (9) for the realizable model) equations. The boundaries of the computational domain are split into inflow, outflow, bottom and top.
• On the inflow boundary a vertical profile is imposed for inflow velocity u and turbulence unknowns k and ε. The profiles are generated from a single-column (1D) precursor simulation (i.e. flat terrain and uniform roughness) . • On the outflow boundary geostrophic pressure and no shear stress are imposed for the momentum equation and symmetric boundary conditions (no gradient) are imposed for the turbulence unknowns. • On the top boundary symmetry boundary conditions are imposed for tangential velocity and turbulence unknowns. The normal velocity component is fixed to zero (i.e. u · n = 0) and pressure is set to geostrophic. • On the bottom boundary a wall law satisfying the Monin-Obukhov equilibrium is imposed to the momentum and turbulence equations removing a boundary layer of thickness δ w . The imposed shear stress τ w tangent to the wall is expressed in terms of two velocity scales, namely u * v and u * k , based on the tangent velocity and the turbulent kinetic energy respectively: where u is the component of the velocity tangent to the wall and |u| denotes its norm. The friction velocity u * v is obtained from the neutral atmospheric velocity profile at a distance δ w from the wall, being κ the Von Karman constant and z 0 the roughness length of the terrain. Finally, zero diffusion through the wall is imposed for the turbulent kinetic energy (∇k · n = 0) and ε is imposed as:

Hybrid mesh generation
The objective is to generate high-quality conformal hybrid meshes specifically designed for the simulation of ABL flows over complex terrains using actuator discs. The required mesh inputs include topography and terrain roughness, turbine characteristics (location, diameter and hub height), and wind inflow direction to determine the orientation of the actuator discs.
To compute the wind inflow direction for each turbine we perform a precursor simulation using the background ABL mesh without turbines. Given these data and parameters, the hybrid mesh generation procedure is fully automatic and consists on 3 main steps (see [10] for further details): i) Surface mesh generation. First, a 2D flat surface quadrilateral mesh is generated with three differentiated zones, namely farm, transition and buffer zones. The farm zone contains the area of interest and has a smaller element size. The transition zone surrounds the farm and has elements of increasing size outwards. Finally, the external buffer zone has coarser elements and is designed to accommodate the inflow and outflow conditions. This 2D flat mesh is then projected to the topography to obtain a 3D surface mesh. The underlying topography and terrain roughness are assimilated from a data file having any of the standard formats (point cloud, cartesian grid, contour levels, etc.) and filtered to remove data noise using the signal processing smoothing method presented in [11]. In certain zones with high topographic gradients, this simple projection can result on low quality elements (see Fig. 1(a)). For this reason, the 3D surface mesh is optimized maximizing the elemental quality presented in [12,13] and imposing the nodes position over the exact topography to avoid lost of real geometry representation, using the process presented in [14,15]. Fig. 1(b) shows the surface mesh after optimization, where a significant quality improvement is clearly observed. Only the transition and farm zones contain assimilated topography and roughness data whereas, in contrast, the buffer zone is left flat to guarantee consistence with the inflow profiles.
ii) Semi-structured volume mesh generation. Second, a semi-structured volume mesh of hexahedral elements is generated extruding the surface mesh into layers using a desired geometrical growing ratio to have higher vertical resolution near the ground. The extruding direction for a given node is computed as the pseudo-normal direction [16], using an average normal of the adjacent elements, that maximizes the orthogonality of the new generated layer. Each extruding step is combined with a non-linear optimization [17,18] of the element quality to improve the mesh configuration before generating a new layer of hexahedra (see Fig. 1(d)).
iii) Actuator disc insertion and hybrid mesh. Finally, the third step consists on inserting the discs and refining the mesh around and downstream each disc conserving conformity. Firstly, the area surrounding each turbine is emptied. The disc is discretized using hexahedra, and the resulting disc mesh is extruded upstream, downstream and radially. Different element growing factors are used across each direction in order to smoothly tend to the element size of the background mesh in the different regions. At this stage, the background mesh and the mesh surrounding each disc are disconnected and pyramids and tetrahedra are used to perform a conformal union, thus leading to a final hybrid mesh. Pyramid elements are generated at the faces of the hexahedra facing other hexahedra, while tetrahedra are generated at the faces facing the gap between meshes. In this way, the gap between the different meshes is now bounded by triangles, which it is exploded to conform using the tetrahedral mesher TetGen [19]. The process is depicted in Figure 2 for a single turbine case.

Actuator disc model
The force exerted by the wind turbine over the flow is modeled as a uniformly loaded disc by means of the force term Ct 2∆ U 2 ∞ n d , where the thrust coefficient C t is supplied by the manufacturers as a thrust coefficient curve depending on the undisturbed wind velocity C tm (U ∞ ). However, in the case of wind farms, there is no obvious approach to estimate the free stream velocity U ∞ (and therefore C t ) because the wind turbine power and thrust curves are usually provided for single-machine operation rather than operation in the wake of another turbine. For this reason, we relate the free stream velocity U ∞ to the velocity at hub height U hub in terms of the thrust coefficient using one dimensional momentum theory. For high thrust coefficients momentum theory is no longer valid and we use the empirical relationship developed by Glauert [20] to obtain the theoretical thrust coefficient C ta : where a is the axial induction factor. The velocity at hub height U hub is calculated as the wind velocity component perpendicular to the disc surface averaged over the entire disc volume. To compute the proper value of U ∞ (and C t ) it is standardly posed an iterative procedure [3]. Given an initial guess for U ∞ it is first computed C tm (U ∞ ), following, the induction factor a is updated in terms of C ta = C tm (U ∞ ) (inverse of Eq. (14)) and, finally, a new U ∞ is computed in terms of a (from Eq. (15)) until U ∞ converges to a fixed value (verifying C tm = C ta ). However, this iterative scheme diverged in some of the tested complex cases, even using relaxation methods. Herein, introducing (15) in (14) we rewrite the iterative problem in terms of U ∞ and we translate it into solving the non-linear equation To solve (16) we use the bisection method, exploiting that the equation is one-dimensional and avoiding to compute the derivatives of the target function f . We found this nonlinear method to always converge to a proper value of U ∞ in all the tested cases.

Numerical method: linearization, discretization scheme and linear solvers
The RANS models have been implemented in Alya, an in-house high performance computing (HPC) code able to run large-scale applications. The code was recently tested on 100, 000 processors with a parallel efficiency above 90% [21]. The Navier-Stokes and turbulence equations are discretized using a stabilized finite element method using equal interpolation for all the unknowns. As stabilization scheme we used the Algebraical Subgrid Scale method (ASGS) [22] extended for nonlinear equations [23], which gives stability to convection and Coriolis dominating terms in the momentum equation and to convection and reactive terms in the turbulence equations, removing spurious oscillations. The ASGS stabilization method gives also stability to pressure, allowing equal interpolation spaces for pressure and velocity. The velocity-pressure problem is decoupled using an Orthomin solver [24] that converges to the monolithic scheme.
A robust finite element scheme written in block-triangular form [25] is obtained for the k-ε equations (3)-(4)/(9). A priori, the k-ε equations are "well behaved" since the diffusion and reaction coefficients are positive. However, the numerical scheme cannot always guarantee positiveness (numerical oscillations occur) and sign variations in the reactive terms leading to loss of stability. In order to avoid instabilities and numerical convergence issues, ε and k are not allowed to drop below a predefined limit by applying a clipping. In addition, the innermost iterative loops of the k and ε equations (3)-(4)/(9) are linearized using a Newton-Raphson scheme for the quadratic terms, considering ν t and P k constants within the innermost loops. The stabilization of the reactive terms of the dissipation equation (4)/(9) has been found essential to achieve convergence of the unknowns. Once the algebraical system of equations are obtained, a Deflated Conjugate Gradient [26] solver with a linelet pre-conditioner [27] is used to solve the pressure, and a Generalized Minimizing Residual (GMRES) solver is used for the velocity and turbulence unknowns, leading to un-symmetric problems. Figure 3 shows a basic model validation for the well-known Sexbierum benchmark case [5], in which wind velocity deficits and added turbulence intensities are analyzed at different transects downstream of a single isolated turbine. Four different simulations were performed using the standard modified and realizable RANS models accounting and not accounting for Coriolis forces. The obtained results are in very good agreement with those presented in [5]. As expected, the realizable model gives more accurate results than the standard modified model downstream of the actuator disc. Note also how the effect of the Coriolis force is to increase the wind velocity deficit along the wake.

Validation
4. Application to an onshore wind farm As an illustrative application example, this section shows results for an onshore wind farm located in Spain and consisting on 165 turbines having a diameter of 77 m and a hub height of 80 m ( Fig. 4(a)). The computational domain is 17 × 14 km 2 with the top at 2 km above the highest terrain elevation. The final horizontal mesh resolution was determined after a preliminary mesh convergence analysis without considering the presence of turbines. We found that an horizontal mesh resolution of 25 m guarantees a numeric accuracy of L 2 -norm of the error ( e 2 := ( e 2 ) 1/2 ) below the 0.05% (discretization errors were computed against the solution from a finer mesh). Note from Fig. 4(b) how the mesh convergence of the solver is quadratic even for complex terrain. After introducing the actuator discs, the final hybrid mesh (25 m resolution at surface) is composed of 21.7 M elements (16M hexahedra, 4M tetrahedra, 1M pyramids). The actuator discs were meshed with elements of 11.6 m (15% of the turbine diameter) on the disc plane and assuming a disc thickness of 4 m (6% of its diameter).  The CPU time required to complete the entire meshing process was of 284 seconds using a single thread on a MacBookPro with core i7-2.7GHz. Figures 4(c) and 4(d) show the velocity speed up and the Turbulent Kinetic Energy (TKE) at hub height, respectively. Note that, for this particular wind direction, many wind turbines are affected by the upstream turbines. Figure 5 compares hub wind velocities and TKE from 4 different simulations (2 RANS models with and without Coriolis terms). Note how, when the Coriolis effect is not accounted for, the standard modified model predicts higher hub velocities than the realizable model. This is due to the over-prediction of the wake velocity in the standard model [5]. However, when accounting for Coriolis force, both models predict similar hub velocities and TKE at hub height. The mixing length limitation model counteracts the over-prediction of Reynolds stresses in the standard model behind the discs, leading to larger wakes. When accounting for Coriolis force the use of the mixing length limitation model decreases the obtained TKE at hub height for each RANS model in almost all wind turbines. This effect is stronger for the standard modified model, which predicts the lowers TKE values at hub height. Note also how, when using the realizable model, the effect of accounting for Coriolis force over the obtained wind velocity and TKE at hub height is much smaller than in the standard k-ε model.

Conclusions
A wind simulation framework for offshore and onshore wind farms has been implemented in the Alya HPC solver considering two different RANS models and simulating wind turbines as actuator discs. This includes an automatic hybrid mesh generator tailored to wind farms,    Figure 5. Simulated hub height wind velocities (a) and turbulent kinetic energy (b) at 165 turbines for the standard modified (STD) and realizable (REA) k − ε RANS models with and without Coriolis terms.
i.e. able to reproduce the ABL, capture the underlying terrain, and refine areas upstream and downstream turbines in a conformal way. The presented mesher is fully automatic and features quadratic convergence to the topography.
Regarding the actuator disc model, a new robust nonlinear method has been presented to calculate the thrust coefficient and free stream velocity for each wind turbine. This new approach has proved to converge in all the tested cases, in contrast with the standard iterative process, which was not convergent for some tested complex onshore wind farm configurations. The k-ε realizable model has been modified to account for Coriolis effects in the atmospheric boundary layer. The realizable and the standard modified k-ε models obtain very similar results over flat and homogeneous terrains without wind turbines. However, on complex terrains, the realizable model predicts lower velocities in the wakes of wind turbines, specially when not accounting for Coriolis effects. The obtained results of the standard and realizable models are in closer agreement when accounting for Coriolis force.