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 107, 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 agreemen<sup>t</sup> 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 agreemen<sup>t</sup> with measurements and LES in comparison to the original k–<sup>ε</sup>. 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–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 <sup>m</sup>·s<sup>−</sup>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.

#### *1.2. 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.

#### **2. Methods: Wind Farm CFD Modeling**

## *2.1. 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:

$$
\lambda = \frac{\boldsymbol{\omega} \cdot \boldsymbol{R}}{\boldsymbol{\mathcal{U}}\_{freestream}}.\tag{1}
$$

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):

$$TI = \frac{\sigma\_{U}}{\overline{\mathcal{U}\_{frossstream}}} \tag{2}$$

## *2.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.

**Figure 1.** *Cont*.

**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.

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

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.

**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.

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 × 105. 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.

#### *2.3. 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.

## *2.4. 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–61], supporting the assumptions made in this work.

#### *2.5. 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 <sup>m</sup>·s<sup>−</sup><sup>1</sup> and 15 <sup>m</sup>·s<sup>−</sup>1, pitch angle (θ) within −1◦ and <sup>−</sup>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: (a) The velocities (10 <sup>m</sup>·s<sup>−</sup><sup>1</sup> and 15 <sup>m</sup>·s<sup>−</sup>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 <sup>m</sup>·s<sup>−</sup>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.
