Development of a Computational System to Improve Wind Farm Layout , Part II : Wind Turbine Wakes Interaction

The second part of this work describes a wind turbine Computational Fluid Dynamics (CFD) simulation capable of modeling wake effects. The work is intended to establish a computational framework from which to investigate wind farm layout. Following the first part of this work that described the near wake flow field, the physical domain of the validated model in the near wake was adapted and extended to include the far wake. Additionally, the numerical approach implemented allowed to efficiently model the effects of the wake interaction between rows in a wind farm with reduced computational costs. The influence of some wind farm design parameters on the wake development was assessed: Tip Speed Ratio (TSR), free-stream velocity, and pitch angle. The results showed that the velocity and turbulence intensity profiles in the far wake are dependent on the TSR. The wake profile did not present significant sensitivity to the pitch angle for values kept close to the designed condition. The capability of the proposed CFD model showed to be consistent when compared with field data and kinematical models results, presenting similar ranges of wake deficit. In conclusion, the computational models proposed in this work can be used to improve wind farm layout considering wake effects.


Introduction
The necessity for improving wake models has become more apparent over the last decade with the continuous growth of the wind energy market.Literature shows several analytical wake models: Infinite wind farm boundary layer model, Jensen wake model, Larsen model, dynamic wake meandering model, FUGA (Linearized RANS Model), and EllipSys3D.All these models are excellent tools to estimate wake effects, but there is still room for improvement.Usually, analytical models do not consider wake characteristics according to variable operating conditions.However, Computational Fluid Dynamics (CFD) models have the capability to model wake velocity deficit and Turbulence Intensity (TI) according to variable operating conditions.Although computationally expensive, CFD models are powerful tools that can be applied to solve some of the most complex problems in engineering.This work describes how operational parameters affect the aerodynamic behavior of the near wake of a wind turbine up to 5 diameters downstream of the rotor.Moreover, this study proposes a CFD modeling technique to characterize three-dimensional far wake effects, and numerically quantify the influence of some important wind farm design parameters on the far wake aerodynamic behavior.The literature shows that there is a gap in attempting to solve the Wind Farm Layout Optimization Problem (WFLOP) while still considering a rigorous evaluation of the wake effects.The objective of this work is to develop a CFD model with such capabilities, applicable for future applications related to the WFLOP.
1.1.Review: Wind Farm Aerodynamics 1.1.1.Wake Aerodynamics Wake models are usually divided in literature [1][2][3][4] in two categories: (1) Analytical/empirical/ explicit wake models; and (2) computational/implicit wake models.The analytical models solve a set of equations based on the conservation of mass and empirical relations of wake decay, characterizing the energy content in the flow field, and ignoring the details of the exact nature of the flow field.Kinematic models such as Jensen, Larsen, and Frandsen's model assume self-similar velocity deficit profiles, not solving the turbulence field but only the momentum equation [2].The velocity deficit is derived from global momentum conversation, using thrust coefficient of the turbine as an input [1].The computational models solve the fluid flow equations for the wake velocity and turbulence field, whether simplified or not [2].

Wind Energy Computational Fluid Dynamics (CFD) Review
Although there are many CFD studies in the literature approaching wind energy, this is a field of study still in development.CFD modeling techniques applicable for wind turbines significantly vary in literature, showing that there is no well-stablished standard approach.This section presents a comprehensive literature review in CFD models applicable to wind energy, providing an overview on what has been done prior to this work.In regards to CFD techniques for modeling wind turbine flow field, the goal is to investigate what possibilities have not been explored yet, seeking to develop a novel wind turbine CFD model capable of evaluating far wake aerodynamics characteristics.As previously mentioned, a correct evaluation of such characteristics can help to achieve better solutions for the WFLOP.

NREL (National Renewable Energy Laboratory) Phase VI
Several studies utilized the NREL (National Renewable Energy Laboratory)/NASA (The National Aeronautics and Space Administration) Ames Phase VI experimental data campaign to validate their computational models, all of them using pressure coefficient on the blades and aerodynamic torque data for comparison.However, it is difficult to validate wake flow field since no wake measurements were performed in these experiments.Zhou et al. [5] performed Large-Eddy Simulation (LES) of the NREL phase VI, evaluating the effect of different inflow conditions (using user-defined functions) on aerodynamic loading and near wake characteristics.A structured multi-block mesh (with sliding mesh zone) was implemented with refinement on leading and trailing edges.They found that the wind shear and turbulence effects destroyed the uniform and symmetric wake profile in the far wake.Hsu et al. [6] validated a finite-element (Lagrangian-Eulerian) model of the NREL Phase VI using a non-structured rotating mesh.Wake characterization was not the focus of the study, which explains the wake made out of coarse non-structured cells with no refinement.Gundling et al. [7] evaluated low and high fidelity models using the NREL Phase VI for predicting wind turbine performance, aeroelastic behavior, and wakes: (1) The Blade Element Method (BEM) with a free-vortex wake; (2) the Actuator Disc Model (ADM); and (3) the Full Rotor Method (FRM).No specific information or sketch of the wake was provided or described.The FRM showed the largest wind deficits and the slowest dissipation rate for the far wake.Mo et al. [8] developed a study in more depth to understand wake aerodynamics performing a LES of the NREL Phase VI using the dynamic Smagorinsky model.Additionally, verification of the average TI was performed against an analytical model.They found that the downstream distance where instability and vortex breakdowns occur is dependent on wind free-stream inlet conditions: 7 m•s −1 happens at four rotor diameters, while 15.1 m•s −1 between 11 and 13 diameters.A decrease of the TI happened after instability and vortex breakdowns.The strategy for meshing the physical domain consisted of a virtual wind tunnel with the same dimensions of the NASA Ames with the rotor located at 2 diameters downstream of the inlet with a downstream domain of 20 rotor diameters in length.Choudhry et al. [9] performed a very similar CFD study of the NREL Phase VI using the same computational methods of the study conducted by Mo et al. [8], finding that regions of velocity deficit and high TI are within the high vorticity region.Choudry's study did not specify if the mesh is structured or unstructured.

NREL 5 MW
Many studies have developed CFD models considering the NREL 5 MW wind turbine.Among these studies, Troldborg et al. [10] developed a wake CFD (EllipSys3D) study for the NREL 5MW considering three different models: (1) A fully resolved rotor geometry; (2) the Actuator Line Model (ALM); and (3) the ADM.A comparison for wake properties in uniform and turbulent inflows was performed.All the models correctly predict mean axial velocity within 4 radii downstream of the turbine for laminar inflow.The agreement between ADM and ALM methods is acceptable for the wake deficit.They found that the ADM/ALM model is sufficient to simulate turbines under Atmospheric Boundary Layer (ABL) conditions.Storey et al. [11] implemented a CFD model using a modified actuator technique to develop transient simulations, considering the NREL 5MW turbine.They achieved reduction in the computational time for the simulation while still keeping flow solution fidelity compared to the standard ADM.Seydel et al. [12] performed a Reynolds Averaged Navier-Stokes (RANS) k-ω simulation of the NREL 5 MW to study wake effects between two wind turbines.Réthoré et al. [13] investigated CFD techniques based on permeable body forces including: ADM, ALM, and the Actuator Surface Model (ASM).These approaches can potentially reduce the necessity for mesh refinement next to the rotor.Verification for the ADM in comparison with analytical solution for heavily-loaded turbines demonstrated that the ADM can be a cost-effective way to model wind turbine wake.The verification of the ADM showed that 10 cells per diameter are adequate to describe the near wake flow characteristics, and the cell size becomes less critical in the far wake.The computational domain extends 10 diameters laterally and 25 diameters horizontally, and the wake computational grid is uniformly spaced with cells of the same size.Heinz et al. [14] developed a fluid-structure interaction simulation using EllipSys3D and aero-elastic HAWC2 for the NREL 5 MW considering yaw and standard conditions.Miao et al. [15] developed an unsteady CFD (STAR-CCM+) model for the NREL 5 MW rotor considering yawed flow to investigate wake deviation.The full rotor geometry was modeled considering the NREL 5MW wind turbine, under neutral ABL conditions.Wilson et al. [16] developed a CFD model based on the RANS (OpenFoam and ANSYS Fluent) equations, considering k-ε and k-ω SST (Shear Stress Transport) turbulence model to investigate interactions between wind turbines in neutral ABL conditions.The ADM, the ALM, and the FRM were compared considering the NREL 5 MW.Weipao et al. [17] considered the tilt and cone angle to maximize the power generation of a wind farm for the NREL 5 MW.
Other Topics CFD modeling techniques have been applied for designing and the analysis of floating offshore wind farms.Wu et al. [18] developed a CFD for an offshore floating wind turbine.The near-wake domain is defined as 3D downstream, whereas a 0.5 D distance upstream of the rotor is maintained with constant size mesh cells.Two different approaches for blade meshing were implemented: unstructured tetrahedral and unstructured hexahedral.Theunissen et al. [19] developed a computational and experimental study to optimize the layout of an offshore wind farm array with 80 turbines.Tran et al. [20] developed an unsteady CFD model for a floating offshore, using the software FAST (Fatigue, Aerodynamics, Structure and Turbulence) and Unsteady BEM equations for the analysis.
RANS techniques have been widely implemented in the literature.Zhale et al. [21] performed an unsteady yaw description for a 500 kW rotor modeling the RANS equations using EllipSys3D.A pressure-based incompressible flow was setup, considering an iterative SIMPLE (Semi-Implicit Method for Pressure-Linked Equations) and PISO (Pressure-Implicit with Splitting of Operator) second-order accurate scheme, the turbulence k-ω SST model (good performance for wall-bounded adverse pressure gradient flows).The computational mesh was generated using the software Gridgen, with structured elements.Prospathopoulos et al. [22] developed a RANS k-ω model modified for atmospheric flows, finding that CFD models underestimate near wake deficit even for single-wind turbine wake predictions especially under neutral atmospheric conditions.The accuracy was better for the far wake, and this study also considered the multi-wake interaction considering the case of five turbines in a row.AbdelSalam et al. [23] performed experimental procedure and numerical simulation considering a FRM, RANS k-ε modified for atmospheric flows, 2 MW wind turbine SODAR upstream measurements, and wake LIDAR (Light Detection and Ranging) measurements at downstream distances from 2 to 7 diameters.Boudreau et al. [24] studied the axial-flow and cross-flow configurations operating at respective optimal efficiency, with Reynolds' number around 10 7 , 3D DES (Detached-Eddy Simulation), and Unsteady RANS.Ammara et al. [25] developed a RANS steady CVFEM (Control Volume Finite-Element Method) model, considering a two-row periodic wind farm in a neutral ABL.Frau et al. [26] developed an unsteady CFD (ANSYS CFX) k-ω SST model to compare downwind and upwind configurations for offshore applications, using 9 million to 25 million cells.They concluded that the downwind turbine configuration is better suited for multimegawatt offshore wind turbines.Lann et al. [27] developed a new k-ε model consistent with Monin-Obukhov similarity theory (MOST), comparing it to other k-ε models.Lann et al. [28] developed a k-ε-fP viscosity model applied to one on-shore and two off-shore wind farms, and the results were compared with power measurements.The k-ε model underpredicts the power deficit of the first downstream wind turbine, while the k-ε-fP eddy viscosity shows good agreement with the measurements.The difference becomes smaller for wind turbines further downstream.
Computational models based on ADM and ALM have also been widely implemented in the literature.More recently, the ADS has been developed for some researchers.Models based on ADM and ALM are relatively less computationally demanding than computational models for the FRM, such as RANS and LES.Sarmast et al. [29] developed an ALM using a new vortex code on the Biot-Savart law, and by considering two different wind turbines: Constant and variable circulation along the blades.They concluded that a simplex vortex code has similar results to the ALM and a lower computational cost.Ivanell et al. [30] developed a CFD (EllipSys3D) ALM using 5 million mesh points to evaluate downstream wake flow field characteristics and the tip vortices positioning.Masson et al. [31] developed a RANS k-ε ADM to assess impacts of the variation of operational parameters influencing the turbulent flow around a wind turbine nacelle.Troldborg et al. [32] developed an unsteady RANS ALM to analyze wake interaction between two wind turbines under different degrees of ambient TI: Laminar, offshore, and onshore conditions.The results show the influence of the upstream turbine wakes on external blade loading of the downstream turbines.Makridis et al. [33] developed a CFD model in ANSYS Fluent solving the RANS equations, assuming ADM (based on BEM) and considering complex terrain and neutral atmospheric wind flow.A validation was performed against wake data over flat terrain.Neutral atmospheric flow conditions over a hill were tested and validated.
LES and DES models have been studied and implemented for wind energy applications over the last years.Although computationally more expensive, these models are capable of modeling the transient behavior of wind turbines.Schulz et al. [34] developed a CFD (FLOWer) study of the yawed flow (−50 • to +50 • ) on a generic 2.4 MW using DES.Ivanell et al. [35] studied stability properties of wind turbine wakes using a CFD model based on the LES ALM on the tip vortices of the Tjaereborg wind turbine.Bromm et al. [36] investigated the impact of directionally sheared inflow in the wake development, and analysis of the impact of wakes on energy production and loading on a downstream turbine.A LES was performed using the ALM representation.Storey et al. [37] developed a technique coupling transient wind simulation with an aero-elastic simulation to dynamically model turbine operation and wake structures.A LES with an ADM was performed for that study.Troldborg et al. [38] developed a LES with an ALM technique using 8.4 million grid points to study the near and far wake of a wind turbine at various Tip Speed Ratios (TSR).Lann et al. [39] achieved an improvement for the k-ε model, comparing this model with the original k-ε eddy viscosity model, the LES, and a total of eight field test case measurements.The results showed a better agreement with measurements and LES in comparison to the original k-ε.Transient unsteady models such as LES can account for velocity fluctuations by setting perturbation components using Reynolds stress components.This is important for modeling the fluctuations inherently present at the atmospheric wind.Examples of LES models that simulate wind turbines operating in the ABL can be found in the literature ( [40][41][42][43][44][45]). A full review of LES simulations of wind farm aerodynamics can be found in the literature [46].According to Rodrigo et al. [47], challenges for ABL modeling include relation between enhanced mixing in operational models, role of land surface heterogeneity, development of LES models with interactive land-surface, and climatology of boundary-layer parameters such as stability.
Moreover, on the topic of wind energy CFD techniques, some researchers have incorporated one-dimensional codes based on BEM to their models, developing a combined hybrid approach CFD-BEM.For instance, Choi et al. [48] developed a CFD model using ANSYS CFX for 2 MW wind turbines, using BEM theory for the blade design.The distance from upstream and downstream wind turbines changed from three to seven times the diameter, and obviously power output was affected.Esfahanian et al. [49] developed a CFD model of the NREL Phase II using ANSYS Fluent and BEM improved methodology.Furthermore, in CFD techniques, Gopalana et al. [50] developed a coupled mesoscale-microscale model (WINDWYO) coupled with WRF (weather research and forecasting) model and CFD codes of different complexity in order to assess the power predictions and wake visualization at the Lillgrund wind farm.Rosenberg et al. [51] extended efforts of the Vortex Lattice Method (VLM) to analyze aerodynamics of dual-rotor wind turbines.Sreenivas et al. [52] studied the interaction between two wind turbines (NREL S826 airfoils) operating in tandem for TSR of 2.5, 4, and 7 in a wind tunnel speed at 10 m•s −1 .Larsen et al. [53] reviewed several studies in wake aerodynamics.Mittal et al. [54] developed a CFD model (Tenasi: Finite Volume unstructured flow solver) of a wind turbine at various tip-speed ratios, evaluating the effect of temporal convergence on the predicted thrust and power coefficient.Three turbulence models were evaluated: Spalart-Allmaras, Menter SST two equations, and the DES version of the Menter SST.The results pointed that the DES model is significantly better for predicting velocity components in the wake.AbdelSalam et al. [55] modeled the near and far wake using the RANS rotating reference frame, k-ε turbulence model.A FRM and an ADM were compared, and two additional k-ε previously studied in the literature.Wake results were validated against the 180 kW Danwin (three-bladed), showing good agreement.

Gaps in the Literature
Basically the gap existent in the literature is related to CFD models capable of simulating a whole wind farm.The vast majority of the methods simulate single turbines, and only a few of them simulate more than one rotor.The computational resources may be a limiting factor for that, however the gap related to lack of CFD models to simulate whole wind farms can be overcome in other ways.Section 2.3 shows a novel approach of this work as an attempt to overcome the main gap identified in the literature.In regards to other aspects, there is no well-established approach to computationally model wind farms.The choice for boundary conditions and turbulence models widely vary in research and any pattern was identified.Moreover, lack of experimental data in controlled environments for the far wake do not allow researchers to validate their data and improve wake aerodynamics knowledge.Consequently, it is not possible to accurately evaluate wake CFD models found in the literature.The majority of the experimental data for far wake characterization comes from field experimental data, which are difficult to replicate in computational models.

Wake Effects
The wake of a wind turbine is characterized by decreased velocity and increased TI.There are many analytical methods to estimate the velocity-deficit in the wake, but models based on CFD are robust and reliable.In this work, a CFD model was developed to determine the wake velocity deficit and consequently its influence on the wind farm output power.The TI profile in the wake is also characterized using a CFD solver.A very important design parameter for wind farms is the TSR, which is defined as the ratio between the blade tip speed velocity and the free-stream velocity (Equation ( 1)).The TSR and other parameters such as free-stream velocity are critical to determine wake behavior: where ω is the rotor rotational speed, R is the blade radius, and U is the free stream velocity.
Another important design parameter is the TI.This parameter can be calculated using Equation ( 2):

CFD Model
The wind turbine modeled in this work was adapted from the previously validated wind turbine CFD model from part I [56], the MEXICO (Model Experiments in Controlled Conditions) rotor (4.5 m diameter) [57] tested in wind tunnel.The wind turbine blade geometry (MEXICO rotor) including twist angle was built using SolidWorks, and then imported to the ANSYS Design Modeler to build the other turbine components (tower, hub) and the physical domain (Figure 1).The geometry of the MEXICO rotor blades is shown in Figure 1c, the three-bladed model has three types of airfoil: DU91-W2-250 (20% to 45%), Riso-A1-21 (54% to 65%), and NACA 64-418 (75% until the blade tip).The blade is also twisted, and a pitch angle of −2.3 • was applied for the measurements.Since some of the airfoil data are not publicly available, a reverse engineering process was performed to find the airfoil coordinates.A rectangular physical domain was built, and it was broken into smaller pieces, allowing local wake mesh sizing.The largest rectangle in Figure 1a is an exterior part, and the first rectangle corresponds to the near wake until 2 diameters downstream of the rotor.The wake was simulated with a domain extending 13 diameters downstream of the rotor.The CFD model of this study was adapted from part I of this research [56], which is a validation and near wake analysis of the MEXICO rotor.In part I [56], the wind tunnel inlet is located 7 m upstream of the rotor.In the current study, the same distance was adopted as the length upstream of the wind turbine.The solution of the continuity equation generates fewer amounts of residuals for shorter upstream distances, resulting in a better convergence for the CFD solution.
Energies 2019, 11, x FOR PEER REVIEW 6 of 28 characterized using a CFD solver.A very important design parameter for wind farms is the TSR, which is defined as the ratio between the blade tip speed velocity and the free-stream velocity (Equation ( 1)).The TSR and other parameters such as free-stream velocity are critical to determine wake behavior: Where ω is the rotor rotational speed, R is the blade radius, and U is the free stream velocity.
Another important design parameter is the TI.This parameter can be calculated using Equation (2):

CFD Model
The wind turbine modeled in this work was adapted from the previously validated wind turbine CFD model from part I [56], the MEXICO (Model Experiments in Controlled Conditions) rotor (4.5 m diameter) [57] tested in wind tunnel.The wind turbine blade geometry (MEXICO rotor) including twist angle was built using SolidWorks, and then imported to the ANSYS Design Modeler to build the other turbine components (tower, hub) and the physical domain (Figure 1).The geometry of the MEXICO rotor blades is shown in Figure 1c, the three-bladed model has three types of airfoil: DU91-W2-250 (20% to 45%), Riso-A1-21 (54% to 65%), and NACA 64-418 (75% until the blade tip).The blade is also twisted, and a pitch angle of −2.3° was applied for the measurements.Since some of the airfoil data are not publicly available, a reverse engineering process was performed to find the airfoil coordinates.A rectangular physical domain was built, and it was broken into smaller pieces, allowing local wake mesh sizing.The largest rectangle in Figure 1a is an exterior part, and the first rectangle corresponds to the near wake until 2 diameters downstream of the rotor.The wake was simulated with a domain extending 13 diameters downstream of the rotor.The CFD model of this study was adapted from part I of this research [56], which is a validation and near wake analysis of the MEXICO rotor.In part I [56], the wind tunnel inlet is located 7 m upstream of the rotor.In the current study, the same distance was adopted as the length upstream of the wind turbine.The solution of the continuity equation generates fewer amounts of residuals for shorter upstream distances, resulting in a better convergence for the CFD solution.The strategy for meshing (Figure 2) the physical domain is to build a sphere of influence surrounding each rotor, and break the physical domain into smaller rectangles defining them as the same part in the ANSYS Design Modeler.The sphere of influence option allows for a better convergence of the flow field solution.The smaller rectangles allow the mesh element sizing of the near and far wake to be controlled locally, avoiding gradients in the mesh sizing in the interface of each sub-domain.The flow field solution is determined using the CFD solver ANSYS Fluent 17 (ANSYS, Canonsburg, PA, USA), housed in two computers with 64 GB RAM and 8 processes for each machine.ANSYS Fluent solves the equations of fluid flow and heat transfer by default using a stationary (or inertial) reference frame.However, a Moving (or non-inertial) Reference Frame (MRF) can bring advantages in solving the equations for some problems involving moving parts, such as rotating blades.In those problems, the flow around the moving parts is the variable of interest to be determined.In the case of this work, the region behind the wind turbine corresponding to the wake flow field is the region of interest.The MRF technique models the flow around the moving part as a steady-stead problem with respect to the moving frame, allowing to activate reference frames in selected cell zones.The ANSYS Fluent MRF modeling modify the equations of motion to incorporate additional acceleration terms that occur due to the transformation from the stationary to the moving reference.The main reason for employing a MRF is to solve a problem that is unsteady in the stationary (inertial) frame but steady with respect to the moving frame.In this work, the simulation was performed using a steady state MRF approach, and setting the rotational speed to match experimental conditions.Unlike the ADM, ALM, and ASM approaches, the CFD model of this work is a FRM approach which considers the exact 3D blade geometry, including variable chord length, local twist angle, and blade pitch angle.The boundary layer was solved using 10 inflation layers with a ratio of 1.1 to ensure y+ < 1 next to the blade surface.Even though the full blade geometry was resolved using the CFD model, the solid blade geometry was suppressed from the rotating disc, centrally located at the physical domain.The MRF approach essentially consists in building a central disc (a fluid zone) surrounding the solid three-dimensional blades (solid zone) inside the disc.At this point, there are two physical domains: 1) A solid zone representing the blades; and 2) a fluid zone (central disc) surrounding the blades, which is the central disc.The next step is to subtract the blade domain (solid zone) from the fluid zone corresponding to the central disc.After the subtraction operation, there is no more solid body (blade) inside the central disc, but only the external surfaces (walls) of the full three-dimensional blade geometry, meaning that the interior of the blade is now an empty space.The exterior blade surfaces (three-dimensional blade surface including chord, twist, and pitch) remains in the central disc (fluid zone), behaving exactly in the same way as if the blades The strategy for meshing (Figure 2) the physical domain is to build a sphere of influence surrounding each rotor, and break the physical domain into smaller rectangles defining them as the same part in the ANSYS Design Modeler.The sphere of influence option allows for a better convergence of the flow field solution.The smaller rectangles allow the mesh element sizing of the near and far wake to be controlled locally, avoiding gradients in the mesh sizing in the interface of each sub-domain.The flow field solution is determined using the CFD solver ANSYS Fluent 17 (ANSYS, Canonsburg, PA, USA), housed in two computers with 64 GB RAM and 8 processes for each machine.ANSYS Fluent solves the equations of fluid flow and heat transfer by default using a stationary (or inertial) reference frame.However, a Moving (or non-inertial) Reference Frame (MRF) can bring advantages in solving the equations for some problems involving moving parts, such as rotating blades.In those problems, the flow around the moving parts is the variable of interest to be determined.In the case of this work, the region behind the wind turbine corresponding to the wake flow field is the region of interest.The MRF technique models the flow around the moving part as a steady-stead problem with respect to the moving frame, allowing to activate reference frames in selected cell zones.The ANSYS Fluent MRF modeling modify the equations of motion to incorporate additional acceleration terms that occur due to the transformation from the stationary to the moving reference.The main reason for employing a MRF is to solve a problem that is unsteady in the stationary (inertial) frame but steady with respect to the moving frame.In this work, the simulation was performed using a steady state MRF approach, and setting the rotational speed to match experimental conditions.Unlike the ADM, ALM, and ASM approaches, the CFD model of this work is a FRM approach which considers the exact 3D blade geometry, including variable chord length, local twist angle, and blade pitch angle.The boundary layer was solved using 10 inflation layers with a ratio of 1.1 to ensure y+ < 1 next to the blade surface.Even though the full blade geometry was resolved using the CFD model, the solid blade geometry was suppressed from the rotating disc, centrally located at the physical domain.The MRF approach essentially consists in building a central disc (a fluid zone) surrounding the solid three-dimensional blades (solid zone) inside the disc.At this point, there are two physical domains: (1) A solid zone representing the blades; and (2) a fluid zone (central disc) surrounding the blades, which is the central disc.The next step is to subtract the blade domain (solid zone) from the fluid zone corresponding to the central disc.After the subtraction operation, there is no more solid body (blade) inside the central disc, but only the external surfaces (walls) of the full three-dimensional blade geometry, meaning that the interior of the blade is now an empty space.The exterior blade surfaces (three-dimensional blade surface including chord, twist, and pitch) remains in the central disc (fluid zone), behaving exactly in the same way as if the blades had not been suppressed: External walls.This Energies 2019, 12, 1328 8 of 27 procedure is performed because the remaining central disc is the rotating frame in the MRF approach.The rotating speed is set up to the frame and not the blade itself.The disc evolving the three-blade wind turbine shown in Figure 1b is the reference frame, which is set to rotate at the desired operating condition.The loading on the blades is represented using the rotating central disc from Figure 1b, but careful work has been taken to correctly represent these forces.The model validation process can be found in the first part of this research [56], including blade loading, pressure coefficient on the blades, and near wake velocity flow field.Additionally, more details about the numerical modeling process can be found in part I [56].The process to adapt the geometry from part I [56] to the extended geometry in this work included the use of Ansys Design Modeler functions.The operations utilized to build the physical domain shown in Figure 1 include extrusion, skin, Boolean, pattern (to duplicate the turbines and wake domain), construction of primitives (cylinder), slicing, rotation, and translation.

Second and Third Rows Simulation
In this work, we developed a new method to evaluate the second and third rows of turbines where the outlet of the first row becomes the inlet of the second row.This results in a significant reduction in the computational expenses, since there is no need to simulate multiple turbines at once.Multiple turbines would require a mesh with a significant higher number of elements.For instance, the three first rows would require three times more elements in comparison with our approach.The Moreover, the turbulence model selected was the k-ω SST, which is suitable for swirl flow and was used in the literature studies as their main turbulence modeling technique.Since there is no public information from the reports of the MEXICO experiment regarding the inlet inflow conditions, default values of 5% and 10% were assumed for the inflow TI and the viscosity ratio (VR) at the inlet, respectively.The Reynolds number based on the average chord length is approximately 1.5 × 10 5 .Pressure-far-field boundaries, which require the larger exterior rectangle to achieve convergence, were applied for the lateral and superior boundaries.The turbines were both rotating in a clockwise direction, which prevented the implementation of symmetry boundary conditions to simulate the two-turbine case.Pressure-far-field boundary conditions are suitable to model the lateral boundaries of the physical domain.If the lateral boundaries are placed far away from the region perturbed by the wind turbine fluid flow, the streamlines have a straight direction.Pressure-far-field boundaries have a good numerical convergence and stability for the case of straight streamlines.We also apply a pressure-outlet for the exit boundary, and a special type of wall with no shear for the inferior boundary.The use of a no shear wall intends to reduce the complexity of the problem by eliminating the need for modeling the surface roughness, which would require a much more refined mesh at the bottom of the physical domain.Essentially, the goal of the study was to develop a preliminary model capable of simulating a wake interaction effect in a wind farm.Originally, a no shear wall condition for the bottom was implemented for the wake validation presented in part I [56] of this research because the validation of the wake velocity field was not affected by the roughness of the bottom.In part II, the validated model from part I was adapted keeping the same type of boundary conditions.The implementation of a shear wall through the definition of ground surface roughness is an intended further improvement of the model.The mesh sensitivity study can be found at Appendix A, showing the need for using 10 million cells.Additionally, the dimensions of the cell elements in each of the fluid cell zones from Figure 2 can be found in the Appendix A.

Second and Third Rows Simulation
In this work, we developed a new method to evaluate the second and third rows of turbines where the outlet of the first row becomes the inlet of the second row.This results in a significant reduction in the computational expenses, since there is no need to simulate multiple turbines at once.Multiple turbines would require a mesh with a significant higher number of elements.For instance, the three first rows would require three times more elements in comparison with our approach.The goal of this approach was to propose a method to overcome the challenges pointed out in the section: The vast majority of the methods simulate single turbines.This method has never been applied to solve wind farm before in the literature.It is worth mentioning that there is not necessarily an improvement in terms of computational time, since three sequential simulations to simulate three rows take the same amount of time of the conventional simulation with three times more elements.On average, each simulation for case 3 of Table A1 (Appendix A) takes approximately 10 h.The referred reduction in computational expenses comes from the fact that less expensive computational resources are required to perform such simulations.One of the biggest challenges on wind farm computational modeling is the expensive computational resources required to simulate several rows in a wind farm.The use of the technique introduced in this work allows researchers to simulate the wake interaction effect without the need for expensive computational resources.In other words, there is a reduction in the capabilities (processors) required to develop wake interaction simulations.

Wake Similarity
The wind turbine modeled in this work (the MEXICO rotor) has an extensive wake flow field dataset, which allowed the validation performed in part I [56] of this research.The Reynolds number of a utility-scale turbine is higher than a small wind turbine prototype (such as the MEXICO rotor) mainly because of the differences in the chord length.Matching the Reynolds number of a utility-scale turbine and the MEXICO rotor was not achievable because of the extremely high velocities required to counter balance the difference in chord length.Wake characteristics are highly dependent on the TSR though, as shown in part I [56] of this work.This work relies on the assumption that similarity between large and small-scale turbines in regards to wake characteristics can be ensured by matching TSR operating conditions.Examples of small-scale experiments to study large scale wake aerodynamics can be found in the literature [58][59][60][61], supporting the assumptions made in this work.

Cases of Study and Motivation
In the study case of this work, two turbines were located side by side in the first row (Figure 3).Then, a second row of turbines was placed in a position which was totally aligned with the location of the first row.The downstream distance between the first and second rows was 10 rotor diameters.The operating conditions simulated in this work correspond to wind velocities of 10 m•s −1 and 15 m•s −1 , pitch angle (θ) within −1 • and −3 • , and a TSR within 4 and 10.The main motivation for positioning the second row aligned with the first row was to study the effect of wake interaction on the velocity and turbulence flow field (Figures 4-6).The motivation for positioning the turbines side by side in the first row was to study the effect of designing staggered rows.For instance, an increase in output power could be achieved in the second row of turbines by staggering the second row rotors out of the region affected by the wake of upstream rows.However, there could be consequences regarding increased turbulence levels at these locations because of wake expansion effects for the turbulence flow field.A further discussion in Figure 6 explains the importance of studying side-by-side distances in upstream rows and its effect on turbulence flow field.Moreover, the motivation for selecting the operating conditions in this work had to do with the MEXICO experiment and typical wind farm conditions:

Velocity and Turbulence Intensity (TI) Contours
The intensity of the velocity-deficit decayed along the axial distance downstream of the rotor, however the velocity in the wake did not fully recover its free-stream value even after more than 10 diameters downstream of the rotor.Figure 3 shows time-averaged velocity contours for the two-turbine case when considering the designed aerodynamic condition for this specific wind turbine (U = 15 m•s −1 , λ = 6.6, ω = 424.5 rpm, θ = −2.3• ).The region in red (15 m•s −1 ) represents the area where the velocity was not affected by wake effects.On the other hand, the velocity-deficit in the wake of the wind turbine is represented by green and yellow contours.
however the velocity in the wake did not fully recover its free-stream value even after more than 10 diameters downstream of the rotor.Figure 3 shows time-averaged velocity contours for the two-turbine case when considering the designed aerodynamic condition for this specific wind turbine (U = 15 m•s −1 , λ = 6.6, ω = 424.5 rpm, θ = −2.3°).The region in red (15 m•s −1 ) represents the area where the velocity was not affected by wake effects.On the other hand, the velocity-deficit in the wake of the wind turbine is represented by green and yellow contours.Wind farms experience effects from the interaction between wakes from the different rows, which changes the velocity and turbulence flow field.The region free of wake effects became smaller after each row of turbines.Figure 4 shows the velocity contours for a hypothetical second row of wind turbines, while Figure 5 shows TI contours.These simulations considered the designed operational conditions (U = 15 m•s −1 , θ = −2.3°,and TSR = 6.6) for both the first and second row of turbines.Instead of simulating 4 turbines, the methodology applied used data from the previous simulation (Figure 3) for the velocity inlet.Basically, the pressure-outlet of Figure 3 became the velocity-inlet profile for the simulation from Figure 4.This procedure significantly improved the Wind farms experience effects from the interaction between wakes from the different rows, which changes the velocity and turbulence flow field.The region free of wake effects became smaller after each row of turbines.Figure 4 shows the velocity contours for a hypothetical second row of wind turbines, while Figure 5 shows TI contours.These simulations considered the designed operational conditions (U = 15 m•s −1 , θ = −2.3• , and TSR = 6.6) for both the first and second row of turbines.Instead of simulating 4 turbines, the methodology applied used data from the previous simulation (Figure 3) for the velocity inlet.Basically, the pressure-outlet of Figure 3 became the velocity-inlet profile for the simulation from Figure 4.This procedure significantly improved the computational efficiency of the simulation with regards to computational time and convergence, since two turbines were simulated instead of four.The second row of wind turbines were not staggered from the first row of turbines, this way occupying a region affected by wake effects from the upstream first row.The wake velocity contours in Figure 4 show a smaller region of unaffected velocity in comparison with Figure 3, meaning that the region free of wake effects becomes smaller after each row of turbines.
since two turbines were simulated instead of four.The second row of wind turbines were not staggered from the first row of turbines, this way occupying a region affected by wake effects from the upstream first row.The wake velocity contours in Figure 4 show a smaller region of unaffected velocity in comparison with Figure 3, meaning that the region free of wake effects becomes smaller after each row of turbines.

Velocity and TI Plots
Wake data plots for the wake of the first and second turbine rows are shown in Figure 6, comparing the behavior of the velocity deficit (Figure 6a) and the TI (Figure 6b) at a radial traverse at 10 D (diameters) in the wake.The velocity deficit existing in the wake of the second row was slightly higher than the velocity deficit found in the wake of the first row (Figure 6a).The TI of the second row of turbines was considerably higher compared to the same downstream position (10 D) of the first row (Figure 6a).Moreover, interestingly there was an increase of the TI (Figure 6b) in the region between the two turbines (at r = 0), which can be attributed to wake expansion of the turbulence flow field.This can be extremely relevant in the context of wind farm layout optimization, since there is need for improving the turbine packing factor in a wind farm to take the highest benefit/output out of the windiest sites.For instance, the region between the two turbines would not be locally affected with reduced velocities, which could lead to a misleading decision of installing turbines at this position.However, the TI would have increased levels which could have an impact on the components (blades, tower, and turbine) fatigue lifetime.The increase in TI caused by wake interaction effects becomes much more significant as the lateral distance between turbines in the same row decreases.Figure 6b shows that there was a severe increase in TI for the first and second rows when the lateral spacing between turbines was too small (2 D), and Figure 6a shows that even the wake velocity profile was affected.Such effects tend to dissipate for larger lateral distances, as shown by Figure 6e,f which considered 4 D of lateral spacing.Further studies on staggered wind farms should address the influence of the spacing between upstream turbines on turbulence flow field characteristics at the wake, aiming to determine the areas where the level of TI is reduced.If chosen correctly, the side by side distance (from upstream rows) could result in a wake region in which TI levels are reduced, consequently these spots would be more suitable to place downstream turbines.

Velocity and TI Plots
Wake data plots for the wake of the first and second turbine rows are shown in Figure 6, comparing the behavior of the velocity deficit (Figure 6a) and the TI (Figure 6b) at a radial traverse at 10 D (diameters) in the wake.The velocity deficit existing in the wake of the second row was slightly higher than the velocity deficit found in the wake of the first row (Figure 6a).The TI of the second row of turbines was considerably higher compared to the same downstream position (10 D) of the first row (Figure 6a).Moreover, interestingly there was an increase of the TI (Figure 6b) in the region between the two turbines (at r = 0), which can be attributed to wake expansion of the turbulence flow field.This can be extremely relevant in the context of wind farm layout optimization, since there is need for improving the turbine packing factor in a wind farm to take the highest benefit/output out of the windiest sites.For instance, the region between the two turbines would not wake velocity profile was affected.effects tend to dissipate for larger lateral distances, as shown by Figure 6e and 6f which considered 4 D of lateral spacing.Further studies on staggered wind farms should address the influence of the spacing between upstream turbines on turbulence flow field characteristics at the wake, aiming to determine the areas where the level of TI is reduced.If chosen correctly, the side by side distance (from upstream rows) could result in a wake region in which TI levels are reduced, consequently these spots would be more suitable to place downstream turbines.The evolution of the near wake (up to 6 D) of a single turbine is shown in Figures 7 and 8 for two different free-stream and TSR values (U = 10 m•s −1 , U = 15 m•s −1 , TSR = 4 and 6.6).The velocity-deficit increased as the TSR increased from 4 to 6.6 for all the positions considered in the wake (Figure 7).In regards to a TSR = 4 and considering U = 10 m•s −1 , the wake velocity deficit had a peak of approximately 15% at x/D = 3 in the near wake, and the velocity deficit decreased at x/D = 6 to approximately 11%.The case of TSR = 6.6 and U = 10 m•s −1 presents a velocity deficit peak of 25% at The evolution of the near wake (up to 6 D) of a single turbine is shown in Figures 7 and 8 for two different free-stream and TSR values (U = 10 m•s −1 , U = 15 m•s −1 , TSR = 4 and 6.6).The velocity-deficit increased as the TSR increased from 4 to 6.6 for all the positions considered in the wake (Figure 7).In regards to a TSR = 4 and considering U = 10 m•s −1 , the wake velocity deficit had a peak of approximately 15% at x/D = 3 in the near wake, and the velocity deficit decreased at x/D = 6 to approximately 11%.The case of TSR = 6.The evolution of the near wake (up to 6 D) of a single turbine is shown in Figures 7 and 8 for two different free-stream and TSR values (U = 10 m•s − 1 , U = 15 m•s − 1 , TSR = 4 and 6.6).The velocity-deficit increased as the TSR increased from 4 to 6.6 for all the positions considered in the wake (Figure 7).In regards to a TSR = 4 and considering U = 10 m•s − 1 , the wake velocity deficit had a peak of approximately 15% at x/D = 3 in the near wake, and the velocity deficit decreased at x/D = 6 to approximately 11%.The case of TSR = 6.6 and U = 10 m•s − 1 presents a velocity deficit peak of 25% at x/D = 3 and 17.25% at x/D = 6, which was 9% and 6.25% smaller than the values for U = 10 m•s −

Far Wake Aerodynamics: Influence of Operating Conditions
The problem of optimizing a wind farm layout is very complex, therefore assumptions for the operational conditions are important to allow finding a solution to this type of problem.This explains

Far Wake Aerodynamics: Influence of Operating Conditions
The problem of optimizing a wind farm layout is very complex, therefore assumptions for the operational conditions are important to allow finding a solution to this type of problem.This explains the importance of this section; it is very important to verify the range of validity of the solution from the optimization routine.In this section, the influence of some important operating design parameters on the velocity deficit and the TI profile in the far-wake development was analyzed including: TSR, pitch angle (θ), and free-stream velocity (U).

Influence of the Tip Speed Ratio (TSR)
The TSR (or λ) critically influenced the far wake behavior.The velocity deficit increased as the TSR increased from 4 to 10, according to the plots from Figure 9 for axial velocity for a radial traverse in the wake at 10 D (diameters) axial location downstream the rotor.Comparing the two values of TSR from Figure 9a, the highest TSR value (λ = 10) presented the highest velocity-deficit in the far wake behavior for the downstream position considered.Consequently, the TSR was a critical design parameter affecting the three-dimensional extension of the wake.This parameter must be considered to determine the minimal distances between rotors, since a wind farm experiences several different operational conditions with regards to TSR.The TSR (λ) also critically influences the TI in the far wake (Figure 9b), increasing the TSR means that the TI will increase too.

Influence of the Pitch Angle
The pitch angle (θ) had little influence on the velocity and TI profiles in the far wake if the values were kept close to the designed condition.On the other hand, the wake profile was influenced by pitch angle values beyond the designed one.The MEXICO rotor designed condition (θ = −2.3°,U = 15 m•s −1 , and TSR = 6.6) would result in the best aerodynamic performance when the rotor operates under this specific condition.Part I of this research [56] simulated the MEXICO rotor for pitch angle values ranging from +2.3° to −3°, showing that the velocity deficit in the near wake was significantly higher for pitch angle values close to the designed condition.Particularly, a pitch angle of −3° would result in a velocity deficit three times higher than a pitch angle of +2.3°.These results are expected, since the axial induction of the rotor was higher for the designed condition because more energy was being extracted from the incident wind.For a value of +2.3°, the near wake velocity deficit was lower because of the lower rotor axial induction.In this work, three different values of pitch angle were tested (Figure 10), considering the same free-stream velocity and TSR conditions for all of them.The idea was to check the effect of the variation of the pitch angle on the far wake profile.All three pitch angle values tested were close to the designed condition (θ = −2.3°).The velocity profile (Figure 10a) remained the same at 10 diameters downstream the rotor for all the pitch angles values, whereas

Influence of the Pitch Angle
The pitch angle (θ) had little influence on the velocity and TI profiles in the far wake if the values were kept close to the designed condition.On the other hand, the wake profile was influenced by pitch angle values beyond the designed one.The MEXICO rotor designed condition (θ = −2.3• , U = 15 m•s −1 , and TSR = 6.6) would result in the best aerodynamic performance when the rotor operates under this specific condition.Part I of this research [56] simulated the MEXICO rotor for pitch angle values ranging from +2.3 • to −3 • , showing that the velocity deficit in the near wake was significantly higher for pitch angle values close to the designed condition.Particularly, a pitch angle of −3 • would result in a velocity deficit three times higher than a pitch angle of +2.3 • .These results are expected, since the axial induction of the rotor was higher for the designed condition because more energy was being extracted from the incident wind.For a value of +2.3 • , the near wake velocity deficit was lower because of the lower rotor axial induction.In this work, three different values of pitch angle were tested (Figure 10), considering the same free-stream velocity and TSR conditions for all of them.The idea was to check the effect of the variation of the pitch angle on the far wake profile.All three pitch angle values tested were close to the designed condition (θ = −2.3• ).The velocity profile (Figure 10a) remained the same at 10 diameters downstream the rotor for all the pitch angles values, whereas there was no significant variation between θ = −2.3• and θ = −3 • for the TI profile (Figure 10b).Still considering the TI profile (Figure 10b), the case of θ = −1 • showed little deviation from the designed condition θ = −2.3• .This means that the pitch angle may be disregarded for an optimization routine.At least in a preliminary analysis, the pitch angle of individual rotors could be set to the designed condition in order to have the best aerodynamic performance.This could be very important tackling such a complex problem of optimizing wind farm layout, since it is desired to reduce the associated number of variables as much as possible.Figure 10a shows that the velocity wake profile would not be severely affected by doing that, and Figure 10b shows that the effect on the TI profile would be limited to less than a 10% increase.It is important to emphasize again that pitch angle values considerably different than the designed condition would severely affect the wake by altering the velocity and turbulence wake profiles, as shown in part I of this research [56].

Influence of the Free-Stream Velocity
Increasing/decreasing the free-stream velocity value did not affect the magnitude (percentage) of the velocity deficit (Figure 11a).On the other hand, increasing the free-stream velocity value greatly affected the magnitude of the TI (Figure 11b).Consequently, it is important to consider variable free-stream velocity conditions to verify that the optimal wind farm layout solution is not sensitive to the variation of the velocity.Since the turbine components lifetime was closely related to the TI conditions, the variation of the free-stream velocity could be a critical factor to determine the payback of a wind farm.

Influence of the Free-Stream Velocity
Increasing/decreasing the free-stream velocity value did not affect the magnitude (percentage) of the velocity deficit (Figure 11a).On the other hand, increasing the free-stream velocity value greatly affected the magnitude of the TI (Figure 11b).Consequently, it is important to consider variable free-stream velocity conditions to verify that the optimal wind farm layout solution is not sensitive to the variation of the velocity.Since the turbine components lifetime was closely related to the TI conditions, the variation of the free-stream velocity could be a critical factor to determine the payback of a wind farm.
of the velocity deficit (Figure 11a).On the other hand, increasing the free-stream velocity value greatly affected the magnitude of the TI (Figure 11b).Consequently, it is important to consider variable free-stream velocity conditions to verify that the optimal wind farm layout solution is not sensitive to the variation of the velocity.Since the turbine components lifetime was closely related to the TI conditions, the variation of the free-stream velocity could be a critical factor to determine the payback of a wind farm.

Discussion
The MEXICO experiment tested specific operating conditions, providing experimental PIV (Particle Image Velocimetry) measurements for the velocity flow field in the near wake.Far wake measurements in a controlled (wind tunnel) was not viable because of technical constraints, which explains why the MEXICO experiment did not provide such measurements.Although it is impossible to compare the far wake data from the simulations performed in this work against experimental data

Discussion
The MEXICO experiment tested specific operating conditions, providing experimental PIV (Particle Image Velocimetry) measurements for the velocity flow field in the near wake.Far wake measurements in a controlled (wind tunnel) was not viable because of technical constraints, which explains why the MEXICO experiment did not provide such measurements.Although it is impossible to compare the far wake data from the simulations performed in this work against experimental data from the MEXICO rotor, the results of the velocity in the far wake can be verified by comparing other relevant works in the literature.Such comparisons show that the results of the CFD model developed in this work consistently agree with previous literature studies.It is worth mentioning that the following studies cited in this discussion neither computationally implemented the MEXICO rotor nor the exact same operating conditions than the ones implemented in this current work.Despite that, the comparisons presented in this discussion intend to verify an agreement with an acceptable range of consistency.
For instance, the range of velocity deficit found in this study agreed with other works in literature.Some of these studies specified operating conditions in which the results were obtained.Mittal et al. [54] analyzed operating conditions of TSR = 6 at x/D = 5, finding 40% of velocity deficit and a not completely symmetrical radial curve profile.Asymmetry was attributed to interaction with the tower.For the off-design condition of TSR = 3, the velocity deficit had a peak of 30% at x/D = 5.For the off-design condition of TSR = 10 showed a peak of 80% of the velocity deficit at x/D = 1, and 40% at x/D = 5.Although the velocity deficit values were slightly greater than the ones found in this work, there is consistency between the results found in this work (Figure 7), in which the estimated velocity deficit showed a value of approximately 25% for a TSR = 6.6 and 15% for a TSR = 4. Storey et al. [37] found that as the free-stream velocity increased and the TSR decreased (rpm maintained constant), the overall velocity deficit decreased.These results confirm the trend found in Figure 7, Figure 8, and Figure 9: A decrease of the TSR, results in a decrease of the velocity deficit.Additionally, Storey et al. [37] found that the shape and magnitude of the velocity deficit vary significantly with the wind speed and TSR.The expansion of the wake varies with wind speed, confirming the trend observed in this work in Figure 11b: The wake presents more pronounced expansion as the velocity varied from 10 m•s −1 to 15 m•s −1 .
Moreover, even though some of the works in the literature did not specify operating conditions regarding TSR, their data consistently agree with the results for the velocity deficit found in this work.Prospathopoulos et al. [22] considered a downstream spacing of 5 D between the turbines, finding a Energies 2019, 12, 1328 19 of 27 velocity deficit in the wake of 40% at 2.5 D and 30% at 3.5 D for the stable stratification case of the Energy Research Centre of the Netherlands (ECN) measurements.Those results are similar to the values found in this present work (Figure 7).Gundling et al. [7] modeled wake wind speed deficits for different wake models and compared them.The UWAKE (Free-vortex Wake with a BEM) model has a maximum velocity deficit of 67% at 4R for 5 m•s −1 , and a maximum velocity deficit of 81% at 7R for 10 m•s −1 .The FLOWYO (ALM and ADM) LES has a maximum velocity deficit of 70% at 3R for 5 m•s −1 , and a maximum velocity deficit of 83% at 6R for 10 m•s −1 .The wake deficit is similar for the FLOWYO and UWAKE, but little diffusion of the wake was found when using FLOWYO RANS.The diffusion in the wake is similar using UWAKE and FLOWYO LES, while the FLOWYO RANS had not enough turbulent eddy-viscosity produced by the ADM to result in a similar wake diffusion compared to FLOWYO LES and UWAKE.The HELIOS DES model has a velocity deficit of 58% at 4R for 5 m•s −1 , and 65% at 9R for 10 m•s −1 .All those values are within an acceptable range when compared to the values found in this work (Figures 7 and 8).Troldborg et al. [10] analyzed different turbulent inflow conditions for a FRM, an ALM, and an ADM approach.The ALM and the ADM showed the same results for velocity deficit in the wake.A maximum peak of approximately 60% was found at 2 R, which remained almost constant in the same value up to the 10 R analyzed.The FRM showed the same 60% velocity deficit at 2 R, but decreasing the peak value to approximately 50% at 10 R. Those results are similar to the ones found in this research for the wake characteristics (Figures 7 and 8).
Furthermore, a similarity in the shape/format of the wake profile was identified.Mo et al. [8] determined velocity profiles at the wake for several downstream positions.A near-symmetrical but not completely at the blades location, the vertical wake velocity profile had a clear W shape at 1 D and 2 D, and overall the velocity deficit decreased as the free-stream velocity increased from 5 m•s −1 to 15.1 m•s −1 .This was attributed to the state of the completed attached flow in the turbine blade for smaller velocities, and not for more extracted power from the incident wind.Wake shape was not well defined for the further downstream positions.The W shape of the velocity deficit curves is similar to the curve shape found in this present work in Figure 9a for TSR = 4.
The TI behavior showed an increase of its peak value as the TSR increased (Figure 9b).For the cases of wake interaction (Figure 6a,b), there was a considerable increase of TI in the wake of the second row.A verification analysis consistently agrees with the literature results.Troldborg et al. [38] analyzed two turbines case for wake interaction and found that a spacing of 7 D was large enough to allow the wake profile to reach a steady state after the second turbine.AbdelSalam et al. [55] found 65% velocity deficit peak for x/D = 2, 60% velocity deficit peak at x/D = 4, 50% peak velocity deficit at x/D = 6, and 30% velocity deficit peak at x/D = 8.Those results are within an acceptable agreement with the results found in this work .Wilson et al. [16] modeled ADM, ALM, and FRM (full rotor).The wake interaction case showed a strong interaction of wakes when spacing 5 D, which is not easy to compare with our study because the downstream distance implemented is not the same.For a single turbine case, the velocity deficit found by Wilson et al. [16] was slightly higher for ADM than ALM, and MR presented the highest velocity deficit in the wake.The TI was significantly higher for ADM and ALM when compared with FRM.
In regards to field experimental data, Barthelmie et al. [62] studied the influence of the downstream spacing between the turbine rows in the normalized power for the case of the Horns Rev offshore wind farm.Considering 8 m•s −1 and the 2 • sector and a downstream spacing of 7 D, the ratio between the output power of the second and first turbine row was approximately 58%.The output power ratio between the second and third row was 56%.For a downstream spacing of 9.4 D, the ratio between the output power of the second and first turbine row was approximately 70%.The output power ratio between the second and third row was 68%.For a downstream spacing of 10.5 D, the ratio between the output power of the second and first turbine row was approximately 75%.The output power ratio between the second and third row was 70%.Considering 8 m•s −1 and the 30-degree sector, the downstream spacing of 7 D had a ratio between the output power of the second and first turbine row of approximately 80%.The output power ratio between the second and third row was 79%.For a downstream spacing of 9.4 D, the ratio between the output power of the second and first turbine row was approximately 85%.The output power ratio between the second and third row was 80%.For a downstream spacing of 10.5 D, the ratio between the output power of the second and first turbine row was approximately 88%.The output power ratio between the second and third row was 83%.As stated at the beginning of the discussion, operating conditions as well as turbines were not the same of the ones analyzed by Barthelmie et al. [62], such that no direct comparison can be made.However, all those results consistently agree within acceptable levels with the results found in this work .In this work, a downstream distance (spacing) between rows was 10 D, and the estimated ratio between output power of the first and second row was within the range of 73% (for a TSR = 4) and 35% (TSR = 10).A TSR = 4 was a more descriptive operating condition for commercial turbines than a TSR = 10, since turbines are more likely to operate in higher TSR conditions (such as TSR = 10) and only in extreme events such as a wind gust.Figure 6 allows a comparison between the second and third row for a TSR = 6.6, showing a power ratio estimate of 56%.The verification against Barthelmie et al. [62] results is especially important since these are wake field data, showing that the CFD model developed in this work is consistent even when compared against atmospheric measurements.
The discussion of Figure 6 gives a good sense on why to simulate two turbines instead of just one.Figure 6 shows the wake interaction effect that the two-turbine side-by-side case can cause on the TI profile in the wake region between these turbines.The simulation of only one turbine would not allow the detection in the increase in TI levels in the referred region.Moreover, a major objective of this research was to create a computational tool to be implemented in future research aiming to improve wind farm land use by investigating the ratio between output power and area for both aligned and staggered configurations.In such a study, TI effects caused by wake interaction need to be pointed out when proposing improvements for land use.This explains the need for simulating at least two turbines.
There are different approaches to model wind turbines, and some of the most common ones are ADM, ALM, and ASM.The ALM represents the blade by a line, and the ASM represents the blade by a planar surface.The ADM does not model the blades surface, and the rotor is modeled by an infinitesimal disc that represents a discontinuity in pressure.All the actuator approaches (ADM, ALM, and ASM) require knowledge of tabulated lift and drag on the blades, and they require corrections for Coriolis, centrifugal, and tip effects when 2D airfoil data are used [3].According to Vermeer et al. [63], there are many reasons for 2D airfoil data to be corrected to better represent 3D cases.For instance, rotational effects limit the growth of the boundary layer at separation, resulting in increased lift for 3D cases in comparison with 2D cases.Moreover, the drag coefficient can largely differ for 3D cases in comparison with 2D because airfoil characteristics depend on the aspect ratio of the blade.Finally, airfoils under large temporal variations of the angle of attack present hysteresis that changes the static airfoil data.There are some critiques of actuator models.According to Sanderse et al. [3], the ASM requires more accurate airfoil data, as well as knowledge of pressure and skin-friction on the airfoil.According to Churchfield et al. [64], one critique of the ALM is that it does not model the surface of the blade, and in this way this technique is not capable of replicating finer flow features such as boundary layer and separation at high angle of attack.In spite of that, these techniques are still widely used because of the reduced required computational costs.
Previous works in the literature discussed the accuracy of ADM, ALM, and ASM models.Even though the referred models generally represent aerodynamic performance in a satisfactory way, there is a consensus about the need for improvements for properly representing wake flow field under more complex fluid flow configurations.Gundling et al. [7] found that the aerodynamic performance for the NREL Phase VI is well predicted at pre-stalling conditions by low complexity models (BEM and actuator approaches), presenting a similar level of accuracy as higher complexity methods (FRM CFD).The results for high and low fidelity models differ for stalled conditions, but they generally present a good agreement for other less complex situations.The aerodynamic performance at the transition regime was not well predicted by neither low complexity nor high complexity models, but high complexity models accurately predict performance at higher wind speeds when stalling is dominant at the blade surface.The use of a high complexity model (FRM CFD) with adaptative mesh accurately solved the far wake flow field up to a distance of 20 radii downstream of the rotor, but did not show remarkable benefits for performance prediction, in which only the near wake up to 1 radius downstream of the rotor was resolved.Troldborg et al. [10] found that the wake predicted by the ALM and the ADM have very close agreement for uniform inflow conditions, but there is a significant difference when compared with the full resolved rotor method.Additionally, the fully resolved rotor presented higher turbulence levels in the wake.At turbulent inflow, the three methods (ADM, ALM, and FRM) present a close agreement.Troldborg et al. [32] also found good agreement between computational and field data for the turbulence flow field and mean wake deficit.Réthoré et al. [13] compared improved ADM approaches with a CFD FRM, finding significant differences in the turbulence flow field at the wake.The blades and the nacelle from the FRM approach showed a production of turbulence several orders of magnitude higher than the turbulence produced by the ADM.Theunissen et al. [19] accurately performed power calculations using ADM but they found differences in comparison with experimental data for the wake velocity profiles.Storey et al. [37] showed that correct modeling of the ABL and turbulence inflow conditions are important to determine the stability of wind turbine wakes, which may suggest that FRM approach is more suitable for more accuracy on wake prediction.Laan et al. [39] showed that there can be improvements to the actuator approaches by introducing new methodologies to better represent the forces on the rotor.
Even though actuator approaches have still been widely applied in research because of its reduced computational costs, as Sanderse et al. [3] points out, there is need for a detailed study on the influence of exact blade geometry on far wake characteristics.In order to account for the influence of blades geometry on wake characteristics, the FRM CFD approach has been proved as a powerful tool for wind farm design.Unlike the ADM, ALM, or ASM, the CFD model of this work is a direct modeling approach of the rotor which considers the exact 3D blade geometry, including variable chord length, local twist angle, and pitch angle.The boundary layer was solved using 10 inflation layers with a ratio of 1.1 to ensure y+ < 1 next to the blade surface.Therefore, the CFD model of this work does not require the use of tabulated lift and drag data as does the ADM, ALM, and ASM.The capabilities of the model developed in this research allow the evaluation of wind turbine wakes interaction for multiple turbines while still using reduced computational resources.Even though further improvements of the modeling technique implemented in this research might be necessary to account the influence of downstream rows on the wake of upstream rows, the technique is promising towards reducing the computational costs for wind farms simulation.The model is a steady state approach (RANS) aiming to save computational time but can easily be improved to a transient simulation by using a sliding mesh.Moreover, the model of this work has been validated in part I [56] against experimental data for the near wake velocity field.This brings more confidence that the model can properly represent aerodynamics characteristics of the wake flow field.As pointed out by Rodrigo et al. [47], a more realistic description of the wake generation mechanisms in the near wake allows to understand and improve far wake models.
The method implemented in this work applied the outlet from upstream rows as the inlet of downstream rows, with the aim of studying wake interaction effects.This has been performed in this work in a one-way coupling, meaning that only the wake effects from an upstream row will be experienced by downstream rows.A possible effect on upstream rows coming from the interaction with a downstream row (e.g., a second row influencing the far wake of the first row) has not been considered in this work.This might be particularly important for the case of staggered wind farms, which are more densely packed than wind farms with aligned rows.In spite of that, the results presented in this manuscript considered spacing between rows (10 rotor diameters) large enough to dissipate the influence of downstream rows on the wake of an upstream row.A further improvement of the model developed in this research could include the development of a two-way coupling method capable of simulating the effects of downstream rows on the wake of an upstream row, which, as previously mentioned, might matter for the case of staggered farms.Aligned configurations will typically have greater spacing between rows, which potentially dissipates the referred wake interaction effects.
Inflation layers (Figure 2f) were implemented next to the blade surface, keeping y+ < 1 at this location for accurately characterizing boundary layer effects.In ANSYS Fluent, inflation layers can be automatically generated next to selected surfaces, regardless of the surface length (scaled or full-size blades).As specified in the methods section, 10 inflation layers were sufficient to keep y+ < 1 next to the prototype blade surface.In the case of a full-scale wind turbine which operate at higher Reynolds number in comparison with scaled prototypes, a different mesh resolution would be necessary to resolve boundary layer effects.The method implemented in this work would require a larger quantity of inflation layers to achieve mesh resolution small enough to keep y+ < 1 close to the blade surface of a full-scale turbine.The increase in the quantity of inflation layers would require an extra spent, making it more computationally expensive to keep a low y+ value for full-scale rotors.Even though the challenges previously discussed could require more computational costs, they would not prevent the method to be applicable for characterizing full-scale turbines.Moreover, wind turbine wake characteristics have been successfully studied using scaled prototypes such as the one implemented in this work (MEXICO rotor).Whale et al. [65] performed PIV measurements for an untwisted prototype with a flat plate airfoil profile, operating at TSR within 3 and 8 and Reynolds number within 6400 and 16,000.The study showed that wake behavior is insensitive to blade chord Reynolds number, as long as similarity of the TSR is maintained.This is the same hypothesis assumed in this current work (described in Section 2.4), where we assumed that wake similarity can be achieved by matching TSR values of scaled and full-scale turbines.Hossain et al. [61] performed wake measurements of a 500 mm four-bladed turbine using a PIV system, ultrasonic anemometers, and hot-wire anemometers.They performed the same measurements for a 1/10 scaled prototype aiming to study geometry similarity, finding that the wake decays at almost the same levels for the full-scale and the prototype turbines.Ivanell et al. [30] implemented an ADM model for a 0.2 m diameter turbine, showing that the power coefficient is independent of Reynolds number if the Reynolds number is greater than 1000.Additionally, they found that the Reynolds number does not affect the strength of the wake vortex, but only the radial vorticity distribution.Sturge et al. [60] performed an experiment using a wind turbine prototype in which the blade Reynolds number was two order of magnitude lower than the Reynolds number of common full-scale turbines.They presented an interesting discussion about scaled turbines, showing that the Reynolds number became less important for far wake modeling when the ADM was implemented.All these studies mentioned above are examples in the literature that suggest that wake characteristics can be understood by using scaled turbines.

Conclusions
In this work, a CFD model based on the MRF approach was developed to assess wind turbine far wake characteristics according to operating conditions typically experienced in commercial wind farms.The influence of the TSR and free-stream wind speed on wake characteristics such as velocity deficit and TI was discussed and compared with the existing literature on this topic.
This paper reviewed most of the wind turbine wakes studies and wind farm CFD techniques from the literature.Overall, we found that the existing literature studies use different turbulence modeling techniques, as well as CFD solvers with different assumptions and boundary conditions.The wake results vary according to the approach adopted in each work.A gap was identified in the literature review of this work, showing that there is a need for more development of CFD models capable of simulating a whole wind farm.The vast majority of the CFD studies simulate single turbines, and only a few of them simulate more than one rotor.The computational resources may be a limiting factor for that, representing one of the biggest challenges on wind farm computational modeling: The expensive computational resources required to simulate several rows in a wind farm.In order to address this need, this work presented a novel methodology to analyze wind turbine wakes interaction with relatively reduced computational resources.The technique had never been applied before in the context of wind farm numerical modeling.Even though multiple simulations are required for studying the interaction effect between upstream and downstream rows in a wind farm, this work successfully achieved a reduction in computational capabilities (processors) required to perform wake interaction simulations.This represents an advance for wind farm modeling, and many researchers could benefit using such techniques to improve wind energy CFD models.
The model presented in this work was previously validated in part I [56] with regards to near wake data.In part II, a verification of the model against other studies in the literature showed consistency in the wake results within acceptable levels.In regards to the velocity deficit and TI assessment, the values found in this work were similar to other CFD wake studies in the literature.This demonstrates the ability of the proposed CFD model in predicting wake characteristics, and this way the model is ready to be applied for determining the optimal spacing between turbines in a wind farm.The capability of the proposed CFD model showed to be consistent when compared with field data, kinematical models, and CFD results from the literature, showing similar ranges of wake deficit.
Further improvement of the model will include a transient approach modeling to determine wake characteristics according to variable rotor operating conditions.This will extend the capabilities of the proposed model by adding a more realistic modeling approach to derive the aerodynamic behavior of the turbine rows.Moreover, a FSI (fluid solid interaction) model would be relevant to determine how the structural behavior of the blades is affected by variable wind conditions.Although the deformation of the blades will have an impact on the blade fatigue lifetime, no study has previously shown that far wake aerodynamics is significantly impacted by the level of blade deformation.Furthermore, there is still room for improvement of the mesh layout in order to reduce even more the computational resources required for simulating wakes interaction effects.

Figure 1 .
Figure 1.(a) Physical domain with two rotors and boundary conditions; (b) front view of the physical domain, showing the MEXICO (Model Experiments in Controlled Conditions) rotor; (c) MEXICO rotor geometry, a three-bladed rotor with a 4.5 m diameter; and (d) lateral view of the wind turbine, showing the rotating reference frame.

Figure 2 .
Figure 2. (a) Details of the meshing process for the physical domain; (b) lateral view of the mesh showing internal details of the sphere of influence; (c) top view of the mesh showing internal details of the sphere of influence; (d) sectional plane showing details of the mesh on the blade surface; and (e) details of the mesh close to the blade surface, showing inflation layers.

Figure 2 .
Figure 2. (a) Details of the meshing process for the physical domain; (b) lateral view of the mesh showing internal details of the sphere of influence; (c) top view of the mesh showing internal details of the sphere of influence; (d) sectional plane showing details of the mesh on the blade surface; and (e) details of the mesh close to the blade surface, showing inflation layers.
(a) The velocities (10 m•s −1 and 15 m•s −1 ) correspond to values tested in the MEXICO experiment; (b) TSR within 4 and 10 is typically experienced in commercial wind farms; (c) the designed condition for the MEXICO rotor is U = 15 m•s −1 , TSR = 6.6, and pitch angle (θ) = −2.3• ; and (d) part I [56] of this research analyzed a positive value of 2.3 • , confirming that pitch angle values much different from the designed condition strongly influence the wake axial induction.

Figure 3 .
Figure 3. Axial time-averaged velocity contours for the two-turbine case, representing the first row of wind turbines.(a) Lateral view of the wake; and (b) top view of the wake.

Figure 3 .
Figure 3. Axial time-averaged velocity contours for the two-turbine case, representing the first row of wind turbines.(a) Lateral view of the wake; and (b) top view of the wake.

Figure 4 .
Figure 4. Time-averaged axial velocity contours for the two-turbine case in a hypothetical second row of wind turbines.(a) Lateral view of the wake; and (b) top view of the wake.

Figure 4 .
Figure 4. Time-averaged axial velocity contours for the two-turbine case in a hypothetical second row of wind turbines.(a) Lateral view of the wake; and (b) top view of the wake.

Figure 5 .
Figure 5. Turbulence Intensity (TI) contours of a second row of turbines, using a profile from a simulation from a first row of turbines.(a) Lateral view of the wake; and (b) top view of the wake.

Figure 5 .
Figure 5. Turbulence Intensity (TI) contours of a second row of turbines, using a profile from a simulation from a first row of turbines.(a) Lateral view of the wake; and (b) top view of the wake.

Figure 6 .
Figure 6.Wake interaction effect, showing wake data plots in the wake of the first and second rows for different lateral spacing in terms of rotor diameter (D): (a) Axial velocity for 2 D of lateral spacing; (b) TI for 2 D of lateral spacing; (c) axial velocity for 3 D of lateral spacing; (d) TI for 3 D of lateral spacing; (e) axial velocity for 4 D of lateral spacing; (f) TI for 4 D of lateral spacing.The spacing distances refer to hub rotor distances.

Figure 6 .
Figure 6.Wake interaction effect, showing wake data plots in the wake of the first and second rows for different lateral spacing in terms of rotor diameter (D): (a) Axial velocity for 2 D of lateral spacing; (b) TI for 2 D of lateral spacing; (c) axial velocity for 3 D of lateral spacing; (d) TI for 3 D of lateral spacing; (e) axial velocity for 4 D of lateral spacing; (f) TI for 4 D of lateral spacing.The spacing distances refer to hub rotor distances.
6 and U = 10 m•s −1 presents a velocity deficit peak of 25% at x/D = 3 and 17.25% at x/D = 6, which was 9% and 6.25% smaller than the values for U = 10 m•s −1 and TSR = 4.The values of velocity deficit for the case of U = 15 m•s −1 and TSR = 4 were the same of the case U = 10 m•s-1 and TSR = 4, and so were the other two cases (U = 10 m•s −1 TSR = 6.6, and U = 15 m•s −1 and TSR = 6.6) as suggests the self-similar theory.for different lateral spacing in terms of rotor diameter (D): (a) Axial velocity for 2 D of lateral spacing; (b) TI for 2 D of lateral spacing; (c) axial velocity for 3 D of lateral spacing; (d) TI for 3 D of lateral spacing; (e) axial velocity for 4 D of lateral spacing; (f) TI for 4 D of lateral spacing.The spacing distances refer to hub rotor distances.
1 and TSR = 4.The values of velocity deficit for the case of U = 15 m•s − 1 and TSR = 4 were the same of the case U = m•s-1 and TSR = 4, and so were the other two cases (U = 10 m•s − 1 TSR = 6.6, and U = 15 m•s − 1 and TSR = 6.6) as suggests the self-similar theory.

Figure 7 .
Figure 7. Velocity deficit for two different values of Tip Speed Ratio (TSR) and free-stream velocity.

Figure 7 . 28 Figure 8 .
Figure 7. Velocity deficit for two different values of Tip Speed Ratio (TSR) and free-stream velocity.Energies 2019, 11, x FOR PEER REVIEW 16 of 28

Figure 9 .
Figure 9. (a) Axial velocity profile at 10D (diameters) downstream the rotor in the wake of the first row; and (b) TI profile at 10D (diameters) downstream the rotor in the wake of the first row.Different line colors represent the TSR (or λ) of 4 (blue) or 10 (orange).

Figure 9 .
Figure 9. (a) Axial velocity profile at 10D (diameters) downstream the rotor in the wake of the first row; and (b) TI profile at 10D (diameters) downstream the rotor in the wake of the first row.Different line colors represent the TSR (or λ) of 4 (blue) or 10 (orange).

Figure 10 .
Figure10.Influence of the pitch angle (θ) on the: (a) velocity profile at 10D (diameters) downstream the rotor in the wake; and (b) TI profile at 10D (diameters) downstream the rotor in the wake.

Figure 10 .
Figure10.Influence of the pitch angle (θ) on the: (a) velocity profile at 10D (diameters) downstream the rotor in the wake; and (b) TI profile at 10D (diameters) downstream the rotor in the wake.

Figure 11 .
Figure 11.Influence of the free-stream velocity on the: (a)velocity-deficit at 10D (diameters) downstream the rotor in the wake; and (b) TI at 10D (diameters) downstream the rotor in the wake.

Figure 11 .
Figure 11. of the free-stream velocity on the: (a)velocity-deficit at 10D (diameters) downstream the rotor in the wake; and (b) TI at 10D (diameters) downstream the rotor in the wake.