Large Eddy Simulation of an Onshore Wind Farm with the Actuator Line Model Including Wind Turbine’s Control below and above Rated Wind Speed

: As the size of wind turbines increases and their hub heights become higher, which partially explains the vertiginous increase of wind power worldwide in the last decade, the interaction of wind turbines with the atmospheric boundary layer (ABL) and between each other is becoming more complex. There are different approaches to model and compute the aerodynamic loads, and hence the power production, on wind turbines subject to ABL inﬂow conditions ranging from the classical Blade Element Momentum (BEM) method to Computational Fluid Dynamic (CFD) approaches. Also, modern multi-MW wind turbines have a torque controller and a collective pitch controller to manage power output, particularly in maximizing power production or when it is required to down-regulate their production. In this work the results of a validated numerical method, based on a Large Eddy Simulation-Actuator Line Model framework, was applied to simulate a real 7.7 MNW onshore wind farm on Uruguay under different wind conditions, and hence operational situations are shown with the aim to assess the capability of this approach to model actual wind farm dynamics. A description of the implementation of these controllers in the CFD solver Caffa3d, presenting the methodology applied to obtain the controller parameters, is included. For validation, the simulation results were compared with 1 Hz data obtained from the Supervisory Control and Data Acquisition System of the wind farm, focusing on the temporal evolution of the following variables: Wind velocity, rotor angular speed, pitch angle, and electric power. In addition to this, simulations applying active power control at the wind turbine level are presented under different de-rate signals, both constant and time-varying, and were subject to different wind speed proﬁles and wind directions where there was interaction between wind turbines and their wakes.


Introduction
Wind energy has expanded worldwide in the last decade. The installed capacity of wind power has increased close to 10% in 2017 [1], multiplying by almost 4.5 times between 2008 and 2017. Wind energy has become one of the main energy sources in some countries, achieving penetration levels in power generation of more than 30% in countries like Denmark [2] and Uruguay [3]. This development has been possible because of technological improvements making wind turbines more efficient, larger in size (rotor diameter, hub height, and rated power), and more capable of contributing to the electric grid regulation by supplying ancillary services [4][5][6]. In [4], the authors analyzed the evolution of onshore wind turbines installed worldwide between 2005 and 2014, finding that the rotor diameter increased on average around 42.3% in that period with a clear trend towards longer blades. In the same way, the wind turbine's rated power has increased in that period, from 1.38 MW in 2005 to 2.20 MW in 2014 (59%). A similar trend can be observed when looking at the evolution of hub height. In [5], results, with mean power values obtained from ten minutes of data from the Supervisory Control And Data Acquisition system (SCADA). In [28] a subset of five wind turbines of an onshore wind farm consisting of 43 Siemens SWT-2.3-93 wind turbines is simulated using the ADM-NR, ADM-R, and ALM. The study is mainly focused on analyzing the wake characteristics of wind turbines operating subject to an ABL flow, comparing only the mean power ratio of two aligned wind turbines against SCADA data. In [23] the Lillgrund offshore wind farm, consisting of 48 Siemens SWT-2.3-93 wind turbines, was simulated using the ALM to represent wind turbine rotors. Mean power predictions computed for rows of wind turbines were compared with time-averaged power production from ten minutes of SCADA data. The numerical domain used to simulate the wind farm contains about 315 million cells. In [29] the ALM is implemented in a mesh-adaptative fluid dynamic solver, under an unsteady Reynolds-averaged Navier-Stokes-based turbulence modeling approach, where the Lillgrund wind farm was also simulated. For model validation, the focus was made on predicting the power losses along different rows of wind turbines by again comparing the results with ten minutes of power averages obtained from SCADA data.
The authors have been working on using the ALM to simulate actual wind farms to reduce the computational cost by making it more flexible with respect to the mentioned guidelines. In [30] the ALM with a torque controller was used to simulate a row of the offshore wind farm Horns Rev, consisting of 10 wind turbines, while in [31,32] a 7.7 MW onshore wind farm was simulated under different wind directions and wind speeds at hub height, comparing the power output with statistics computed with ten minute averages from the SCADA data. The main contribution of the paper is to go a step further by simulating the same 7.7 MW onshore wind farm under different wind conditions, including a torque controller and a collective pitch controller, comparing the results with 1 Hz of SCADA data from specific events and looking at absolute values of power output as well as the angular speed of the rotors and pitch angle. It should be stated that real mesoscale forcing is not included in the present simulations, like coupling our microscale model with a mesoscale model, so the SCADA data acts as a reference to assess the actual dynamic of this variables. Ongoing research is focused on accomplishing that. A similar approach, but only analyzing power production (relative mean value and standard deviation) for a 30 min period, is presented in [33], where an onshore wind farm, placed in complex terrain and containing 60 1.5 MW wind turbines, is simulated with the ALM, comparing relative power output of different rows of wind turbines with the same statistic computed with 1 Hz of SCADA data. It should be mentioned that a torque controller and a pitch controller are not considered in that work, keeping the angular speed of the rotor and the pitch angle constant. Besides this, ref [34] presents a survey of high fidelity methods to simulate wind farms where the most commonly used codes and the possible case studies to consider are listed. It should be noted that this article indicates that only three of the nine mentioned codes have implemented power control, two of them with the ALM, and in all three cases that is done by coupling the CFD code with an existing aero-servo-elastic model. Unlike this approach, the present article implements torque control and pitch control in the CFD code instead of coupling the code with another tool. It should be mentioned that another limitation of this work is that the main components of the wind turbine are modeled as stiff bodies, including the blades. Currently, work is being done to model the deformation of the blades due to the aerodynamic, gravitational, and inertial loads with the final goal of modeling the interaction between loads and dynamic deformation.
The paper is organized as follows: The next Section describes the numerical method, including the implementation of the Active Power Control (APC). Section 3 presents the validation case, describing the simulated wind farm and giving details of the simulation setup, such as the computational domain, the grid characteristics, and values considered for the power controllers parameters. In Section 4 the results are shown and discussed, focusing on the comparison between numerical and experimental data. The paper ends with conclusions and future work in Section 5.

Numerical Method
In this section the numerical method is described, first presenting the Caffa3d flow solver. Following that we describe the wind turbine model, by means of the Actuator Line Model [17] to represent the rotor, and a conventional torque and pitch controller to follow a desired power reference at wind turbine level.

Flow Solver
A brief description of the CFD code is presented in this section, for further details please see [35,36]. Caffa3d is a second order accurate in the space and time finite volume code, parallelized with the Message Passage Interface (MPI) through domain decomposition. The mathematical model consists of the mass balance Equation (1) and momentum balance Equation (2) (written here only for the first Cartesian direction e 1 ) for a viscous incompressible fluid, including a generic passive scalar transport Equation (3) for scalar field φ with diffusion coefficient Γ. The balance equations are written for a region Ω, limited by a closed surface S, being n S the outward pointing normal unit vector.
where v = (u, v, w) is the fluid velocity, ρ is the fluid density, β is the fluid thermal expansion factor, T is the fluid temperature, and T re f a reference temperature. g is the gravity, p is the pressure, µ is the dynamic viscosity of the fluid, and D is the strain tensor. The generic transport Equation (3) is not used in the present work, nevertheless it allows one to implement further physical models like heat transport, turbulence models, etc. The above equations, which are expressed in their global balance form, with the finite volume method force the conservation of fundamental magnitudes, like mass and momentum in the the solving procedure [37].
The domain is divided in unstructured blocks of structured grids [37] using a collocated arrangement. The same block structure is used for parallelization through MPI. Each block consists of an orthogonal Cartesian grid or a curvilinear body fitted grid, while geometrical properties are always expressed in a Cartesian coordinate system. Flow properties are expressed in primitive variables in the same Cartesian coordinate system, like velocities for example. The immersed boundary method [38] has been implemented in the code, providing more geometrical flexibility.
A discrete approximation of each balance equation of the mathematical model, like (4) which is written for the u velocity component, is obtained by its discretization and linearization at each cell. In (4) the value of the variable at cell center P is related to the values at the six neighbors (W, E, N, S, T, B). Different implicit methods are implemented in the code to approximate the temporal term in the balance equations, like the three time level method and the Crank Nicolson method, while the convective term is discretized with an implicit lower order method and an explicit higher order deferred correction. For further details on the discretization please see [35,36].
An outer iteration loop is performed within each time step, solving sequentially the linear system of each balance equation by an inner iteration loop employing a block structured variant of the Stone-SIP solver algorithm [39] in order to solve and couple the balance equations (please see Figure 1). The outer loop is repeated until the desired level of convergence is achieved. In this work, the Large Eddy Simulation (LES) technique is used to solve the balance equations. To accomplish that, different subgrid scale models are implemented in the code, ranging from the standard Smagorinsky model [40] including different damping functions depending on the surface roughness ( [37,41] for smooth and rough surfaces respectively), and dynamic strategies like the dynamic Smagorinsky model [42], the dynamic mixed Smagorinsky model [43], and the scale-dependent dynamic Smagorinsky model [44].The dynamic models different averaging schemes are included in the code.

Actuator Line Model
To represent wind turbine rotors and their interaction with the wind flow in the simulations, the Actuator Line Model (ALM) has been implemented in the code [45,46]. A brief description of this model is presented below, for further information about its validation and application cases please see [30][31][32][47][48][49][50][51].
In the ALM, instead of representing the wind turbine blades by their exact geometry and solving the blades' boundary layers, the wind turbine rotor is represented as a body force field. A line that rotates with the angular speed of the wind turbine rotor is used to represent each blade. Each line is divided in radial sections where the aerodynamic forces are computed according to Equation (5).
where ρ is the air density, V rel is the relative velocity, c is the chord length, C L and C D are the lift and drag coefficient respectively, e L is a unit vector in the direction of the lift force, e D is a unit vector in the direction of the drag force, and dr is the length of the radial section. Equation (5) requires the chord length and twist angle of the blades to be known, as well as their aerodynamic properties (lift and drag coefficients) in each radial section. The chord and twist angle as a function of the rotor radius are obtained from the wind turbine model, while the lift and drag coefficients are computed from tabulated data of the airfoils used. As mentioned above, the fluid flow perceives the presence of the wind turbine rotor as a body force field by adding a source term in the momentum balance equations. This term is computed as the product of the aerodynamic forces obtained in each radial section and a projection function (please see Figure 1). The projection function smeares the aerodynamic forces onto the computational domain, avoiding numerical instabilities. In this work a Gaussian smearing function Equation (6) is used, taking into account the distance between each grid cell center and each radial section, for three different directions: Normal to the rotor plane, tangential to the blade, and along the blade (n, r, t respectively). This smearing function is consistent with using anisotropic grids and allows the desired ratio of smearing factor and spatial resolution for numerical stability and accuracy to be kept.
The presence of the tower and nacelle are considered through drag coefficients in a similar approach to [16]. At each time step of the simulation, the rotational speed is obtained from the rotor momentum balance Equation (7) , where M gen is the generator torque, ω is the angular speed of the rotor at the current time step, and I is the sum of the inertias of the drive train (rotor, shaft, gearbox, and generator), referred to as the low-speed side. ∆t corresponds to the time step size and ω t−1 accounts for the rotational speed of the previous time step. Finally, the aerodynamic and electric power are calculated as the rotational speed multiplied by the shaft aerodynamic torque and the generator torque, respectively.
A conventional torque and collective pitch controller was implemented in the code based on the NREL 5 MW reference wind turbine [52]. It consists of two control systems which work independently, the generator torque (M gen ) and the blade-pitch controllers, which are described in the following sections.

Generator Torque Controller
At the below-rated rotor speed range, referred to as region 2, the objective of the torque controller is to maximize the power capture [53]. When the rotor speed is below 0.9 × ω 0 , being ω 0 the rated speed, the equation that governs the generator torque controller is: ω accounts for the instantaneous rotor angular speed and K gen is a constant that optimizes the power extraction from the wind and which depends on the aerodynamics and geometrical characteristics of the rotor, defined as in Equation (11), where ρ is the air density, A is the area swept by the rotor, R is the rotor radius, CP opt is the optimal power coefficient, and λ opt is the optimal tip speed ratio. The last two values where computed with the Blade Element Momentum (BEM) method.
At the above-rated wind-speed range, region 3, the goal of the controller is to ensure that the power output is at the desired level, for example at rated power or a given fraction of it, regardless of the fluctuations that the wind speed or the rotor speed may have. The conditions to be at region 3 are either that (A) ω is greater than ω 0 , or (B) the pitch angle is greater than a minimum value θ R3min . At region 3, M gen is computed according to Equation (12). M gen has an upper bound of Sat Factor × M rated , being Sat Factor a saturation factor and M rated the rated torque to Equation (13).
When none of the conditions defining region 3 are fulfilled, and ω is between 0.9 × ω 0 and ω 0 , then M gen is computed as a linear interpolation between M rated and M optimal , according to Equation (14), where M optimal is computed as K × (ω 0 × 0.9) 2 . This region is called 2.5 and is a transition between regions 2 and 3.

Blade-Pitch Controller
This controller operates only at region 3 and its objective is to regulate the generator speed and thus, the rotor angular speed at the rated operation point. At region 2, the pitch-angle of the blades is fixed at its minimum value, θ R2min , to optimize the power extraction from the wind. At regions 2.5 and 3, the rotor-collective blade-pitch-angle values are computed using a proportional-integral (PI) gain-scheduled control on the speed error (∆ω) between the current rotor speed and the rated speed (ω 0 ), as shown in Equation (15).
where θ is the pitch angle and K I and K P account for the proportional and integral gains, respectively. The integral term accounts for the accumulated error over time and in the simulations it is computed by simply adding ∆ω to its previous value in each time step. By linearizing Equation (15), the PI controlled rotor-speed error responds as a second-order system with a natural frequency ω ϕn and damping ratio ζ ϕ . This PI controller ensures that ω fluctuates around its reference value, and thus the active power output fluctuates around the active power rated value, P rated . Based on [52], K I and K P are computed as described in Equation (16).
The blade-pitch sensitivity, ∂P ∂θ , is an aerodynamic property of the rotor that depends on the wind speed, rotor speed, and blade-pitch angle. For the above-rated wind speed range, it can be estimated by applying the BEM method and considering the rated rotational speed and the optimal pitch angle, which imply rated power.

Active Power Control
If the wind speed is high enough, a way to track a power reference signal given by the Transmission System Operator (TSO) can be to de-rate the wind turbines from their rated power. A method to accomplish this is to apply an upper bound to the generator torque of Sat Factor × M DR , instead of Sat Factor × M rated . M DR is computed considering the de-rating command (DR) required by the TSO in Equation (17).
This mode is simply an extension of the controllers presented in the previous Sections 2.2.2 and 2.2.3, where the wind turbine, when it is operating at region 3, limits its power production to the required value when the wind speed is high enough, regulating the aerodynamic torque with the blade-pitch controller. When there is not enough wind to generate the required power, the wind turbine operates at region 2, optimizing the power extraction from the wind by regulating the rotor speed with the generator torque controller.
In the simulations, DR can be a constant signal or vary over time but it needs to be specified for each wind turbine individually, prior to the execution of the simulation. A wind plant global controller is planned to be implemented soon so as to coordinate the operation of the individual wind turbines, accounting for the interactions of wakes based, for example, on the available power of each turbine, as in [54] or [55].

Libertad Wind Farm
As a case study we simulated a 7.7 MW onshore wind farm called 'Libertad', which has been operating since August 2014 by Ventus (www.ventusenergia.com). It is located in the south of Uruguay At selected periods of time, data acquired by the SCADA system of the wind turbines was recorded on a 1 Hz frequency basis, which was used for validating the simulations. Ten-minutes mean and standard deviation of several signals are available for the whole period of operation of the wind farm. The mean free wind velocity and direction was computed from the ten-minutes averages measured at the meteorological mast during the same time periods. In this work we compared results from SCADA data and the simulation of the following four signals from each wind turbine: Electric power, wind speed, rotor angular speed, and blade pitch angle.

Numerical Setup
The computational domain size is 3.00 km × 1.50 km × 0.75 km in the streamwise, spanwise, and vertical directions respectively. It is uniformly divided into 144 × 128 grid cells of size ∆x = 20.8 m and ∆y = 11.7 m in the streamwise and spanwise directions respectively, while a stretched grid containing 80 grid cells in the vertical direction, and a size that ranges between 3 m and 21.3 m. These cell's dimensions imply a resolution of R/2.4, R/4.3 in the streamwise and spanwise directions respectively, being R the rotor radius, while in the vertical direction, 19 grid nodes cover the rotor diameter. Regarding boundary conditions used, at the outlet a zero-velocity gradient was imposed and at the surface a wall model based on the log law was used to compute the stress while periodic boundary conditions were used in the lateral boundaries. The Crank-Nicolson scheme, with a time step of 0.25 s, was used to advance in time. The sub-grid scale stress was computed according to the scale dependent dynamic Smagorinsky model, with a local averaging scheme.
The inflow conditions were obtained from precursor simulations to simulate an ABL like wind flow, using a domain of the same size and resolution than the one used for the main simulations, but without considering the topography and the wind turbines. The ABL wind flow was generated considering periodic boundary conditions at the east and west boundaries, with a constant pressure gradient as a forcing term, and it was run until the wind speed reached statistical convergence. The time evolution of a cross plane of the precursor simulation was considered as a boundary condition at the inlet of the main simulations, as it is usually done for this type of simulation.
For the ALM implementation, the chord length and twist angle distributions were taken from [31], and are depicted in Figure 3a. The airfoil used was the NACA 63-415 for the entire blade. Lift and drag coefficients were also obtained from [31] and the corresponding lift and drag coefficients are presented in Figure 3b. The used smearing factors, normalized by the rotor diameter, are ε n = 0.35, ε r = 0.2, and ε t = 0.17, in the normal, radial and tangential direction respectively. The inertia of the drive-train was estimated based on [57] considering the rotor diameter, the material it is made of, and the wind turbine rated power, obtaining a value of I = 1.72 × 10 7 kgm 2 . Regarding the torque controller, K gen was estimated from the turbines ten-minutes SCADA data, by computing the mean angular speed and mean electric torque at wind speed bins, considering only the below rated range of wind speed. The obtained value was 4.38 × 10 5 Nm ( rad s ) 2 , which differs 25% from the theoretical value in Equation (11), obtained based on the BEM method. In the simulation, the SCADA-obtained value of K gen was used as we think it was a better representation of the actual operation of the wind turbine based on real operational data.
It should be mentioned that the above numerical setup, including a mesh sensitivity study, has been used by the authors in previous studies, with a particular focus on computing the mean power production of each wind turbine at a below rated wind speed, comparing the results with averages obtained after filtering ten-minutes SCADA data according to the simulated wind speed and wind direction. For further details please see [31,32].
The PI blade-pitch controller implementation requires two parameters to be defined, the natural frequency and damping, and the blade-pitch sensitivity. The natural frequency value used is the same as the one considered in [52], but we chose a lower damping value than what was recommended, aiming for the controller response to be similar to what we observed in the 1 Hz of SCADA data. We chose the values of ω ϕn = 0.7 rad/s and ζ ϕ = 0.1. The blade-pitch sensitivity, ∂P ∂θ , was computed as the difference in aerodynamic power (δP) by considering small variations of the pitch angle (δθ) around its optimal value divided by δθ. As suggested in [52], we invoked the frozen-wake assumption in the BEM method, in which the induced wake velocities are held constant while the blade-pitch angle is perturbed. To validate this methodology, we also computed it for the 5 MW reference wind turbine, considering the published data in [52]. It was observed that the pitch sensitivity varied nearly linearly with blade pitch angle, so a linear fit could be computed, where θ K is the blade-pitch angle at which the pitch sensitivity doubled from its value at (θ = 0). Based on SCADA data, the pitch-angle was limited by saturation values of [−2 • , +90 • ], and its maximum rate of change was set to 8 • /s.

Simulation Cases and Results
In this section we present the results of the simulations, focusing on the dynamic of the variables involved in the power controller. The numerical results were compared to 1 Hz of SCADA data acquired at three specifics periods. These periods were selected to show the operation of the wind turbines under various conditions, above and below rated speed, with and without wake interaction as well as with a de-rate operation. It is worth mentioning though that in this work we are not trying to reproduce the exact events that occurred on those specific dates and that we are taking the SCADA data as a reference to show the operation under real conditions. One of the reasons that makes the comparison difficult is that we only have information about the atmospheric conditions at the location from a single meteorological mast, which has cup anemometers at three different heights (please see Figure 2b). Figure 4 shows the temporal evolution of the meteorological variables measured at the meteorological mast and the wind turbines on the day of the event presented in the Section 4.1.
With the meterological mast information it is possible to estimate the vertical profile of the wind speed and turbulence intensity at a given point in the wind farm, but it is not possible to reproduce the same conditions that actually happened. Present and future work is and will be dedicated to accomplish that by coupling the Caffa3d code with the Weather Research & Forecasting (WRF) model. Another fact is that the data acquired by the meteorological mast of the wind farm is recorded on a ten-minutes basis, which is a far larger temporal scale than the one we are considering in the simulations. This means that changes in the atmospheric conditions that may occur within those ten minutes are not recorded by the anemometers, but those changes will imply the controllers to act accordingly. Nevertheless, even if the comparisons between the simulations and 1 Hz of SCADA data are not precise enough because of the mentioned reasons, it is the author's opinion that it is better to present them, rather than presenting only the numerical results. In the second part of this section we present numerical results about active power control at the wind turbine level. To evaluate the effect of wake interaction in the power control, we simulated the operation of the wind turbines under several de-rate signals, comparing two wind directions. Results of the individual and total power production of the wind farm are shown.

Operation without Wake Interaction at above Rated Wind Speed
In this section, the wind farm subject to an ABL like inflow condition, from a direction without wake interaction and wind speed above rated is analyzed. Figure 5 depicts the temporal evolution of 1200 s simulation and SCADA data of WT1 and WT2. The SCADA data corresponds to a period starting on the 5 March 2018 at 15.44 h, when the average wind speed at hub height and wind direction were 13 m/s and 264 • respectively, according to the meteorological mast measurements (please see Figure 4). The wind direction considered at the simulation was 250 • , but it is worth noting that in none of these directions (250 • and 264 • ) were the wind turbines affected by the wake of another wind turbine. The shown signals are wind speed, rotor angular speed, pitch angle, and electric power. The simulation velocity signal corresponds to the value at the cell located at the inlet of the domain, at hub height and upstream of the corresponding wind turbine position. The SCADA signal is the velocity measured at the wind turbine nacelle.
Looking at the time series, the pitch signals take values above the minimum (−2 • ) during the entire 1200 s period. This shows that the pitch-controller was operating to regulate the angular speed of the rotor, which in turn was fluctuating around its rated value in both the SCADA and simulation signals. As a consequence, the power output of both wind turbines was being controlled around its rated value, as desired. Large peaks were observed in the electric power signal from the SCADA of WT2, both above and below rated speed, which lasted for around 10 s. These peaks may be associated to a sudden change in wind speed, as the pitch angle acted accordingly prior and after the occurrence of these peaks. This in not captured in the simulation, where the wind speed did not change in that way over time. The electric power peaks could also be caused by the operation of electric components (e.g., the inverter or the generator), which we did not consider in our simulation.  Figure 6 shows the histograms of the same signals of Figure 5, for WT1, WT2, and WT3. The vertical dashed lines represent the temporal mean value of each signal. We obtained good agreement between the mean values of the simulation and the SCADA data. For the wind speed and pitch histograms, a somewhat different shape can be observed, being narrower and with a higher peak in the simulation. This can be explained by the fact that, in the simulation, the turbulence intensity was lower than in the real event, 8% and 12% respectively. Nevertheless, the general behavior was well captured, showing the capacity of the numerical framework to model the wind turbine's operation, including its controller.

Operation with Wake Interaction at Rated Wind Speed
In this section, a second reference case corresponding to the 15 April 2019, between 11.35 and 12.10 h, is considered. During this period the wind direction was between 120 • and 125 • , where WT1 and WT2 were affected by the wake of WT3 (please see Figure 2). This was noticed in the SCADA data of the affected wind turbines by a reduction of their power production and a wind speed deficit at their location. The free wind speed at hub height was between 9 m/s and 10 m/s, being 10 m/s the rated value. Furthermore, during the entire period WT1 was de-rated to 80% of its rated power, that is 1530 kW. It should be mentioned that, even though in the real event the wind direction was changing, we represented this situation by computing two independent simulations subject to constant wind directions, one at 120 • and the other at 125 • . As pointed out in [58], small changes in wind direction may have large impacts on power production, affecting the wind turbine's operation. Figure 7 shows the histograms of wind speed, angular speed of the rotor, pitch angle, and electric power of the four wind turbines, comparing 1 Hz of SCADA data and the results of the two simulations. For this case, good agreement between the simulations and SCADA data could be observed for WT3 and WT4, which were the ones that not affected by the wake of other turbines. Furthermore, for that reason, notice that the results of the two simulations were similar for these wind turbines. In both simulations the pitch angle of these wind turbines was frequently in its minimum value, showing that the wind turbines were operating at region 2, trying to optimize their power production. Consistently, the rotational speed of their rotors was below its rated value as well as their power output.
For WT1 and WT2, the results of the simulation with wind direction 120 • were quite different to the results of the simulation with wind direction 125 • . This is expected, since there was a higher wake interaction in these wind turbines for wind direction 125 • than for wind direction 120 • . Due to the wake impingement the wind speed seen by the wind turbines was reduced, WT2 was operating at region 2, keeping its pitch angle at its minimum value and adjusting the generator torque to maximize the power capture looking to operate at an optimal tip speed ratio. The rotational speed of the rotor was below its rated value and its power output was lower than the one of WT3. When looking at WT1, there were some differences in its operation, mainly because WT1 was de-rated. As expected, its pitch angle was a bit larger than the minimum value, while its rotational speed was larger than the one of WT2 and close to rated value, particularly in the simulation with wind direction 120 • . The wind speed in both simulations was under-estimated with respect to the SCADA data, presenting differences in the mean values of 13% and 37% with 120 • and 125 • simulations, respectively. Regardless of this, for WT2 the mean values of the SCADA data of electric power, rotor speed, and pitch angle fell between the mean values of the simulations, which shows that on average the operation of this wind turbine was well reproduced.

Operation with Wake Interaction at Below Rated Wind Speed
The third case also corresponds to the 15 April 2019, between 15.35 and 16.20 hrs. During this period the wind direction was also between 120 • and 125 • , but the free wind speed at hub height was about 7 m/s, which is below rated for all wind turbines. In the same way as presented in Section 4.2, wind direction 125 • and wind direction 120 • were simulated to consider the effect of the wind direction on the operation of each wind turbine. Figure 8 depicts the histograms of wind speed, rotational speed of the rotor, pitch angle, and electric power of the four wind turbines. As presented in Section 4.2, the operation of WT3 and WT4 is well captured when looking at their rotor's angular speed, pitch angle, and power production, as these wind turbines were not affected by the wake of other wind turbines. A main difference with respect to the results presented in Figure 7 is that, as the inflow wind speed was below rated, each wind turbine was operating at region 2 keeping its pitch angle at its minimum. The rotational speeds of the rotor of WT3 and WT4 were quite similar, as these wind turbines were subject to similar inflow conditions without the presence of an upstream wind turbine. When looking further downstream, WT1 and WT2 were operating with a lower angular speed than WT3 and WT4, producing less electric power. Again, WT2 was more affected at a wind direction of 125 • than a wind direction of 120 • , with less rotational speed and power output in the former than in the latter.

Active Power Control at Wind Turbine Level
We performed simulations considering several de-rate values for each wind turbine and another wind direction, 147 • where WT3 is clearly affected by the wake of WT4 (see Figure 2b). Figure 9a,b show the operation of WT3 and WT4 respectively, at 147 • , with span-wise averaged velocity of 9 m/s at hub height, considering 6 de-rate commands for both wind turbines: 100%, 90%, 80%, 70%, 60%, and 50%. In the case of WT3 it is clearly noticed that the wind turbine did not reach the required de-rated power, except for some moments when DR = 50% and DR = 60%, when there was enough available power on the wind flow. Furthermore, it can be noticed that the power production of WT3 decreased with higher de-rate values of WT4, as there was less available power due to its wake. The rotor speed was regulated and the pitch angle was kept constant at its minimum value most of the time. In the case of WT4, the required power was reached in all de-rate levels, although for DR = 100%, 90%, and 80% there were periods where the angular speed did not reach rated value, and thus the generator torque controller operated at region 2, regulating rotor speed to optimize the power extraction from the wind. Again, for lower de-rate commands, higher pitch angles were adopted.
The gray zone at the beginning of the time series represents the first 200 s, and indicates the period affected by transient effects, such as the development of the wakes, due to the sudden inclusion of the wind turbines in the simulations. This can be clearly observed in the angular speed, pitch and power signals of WT3, due to the wake propagation of WT4. Figure 10 compares results of WT3 operation at 250 • and 147 • with wind speed at the inlet and hub height of 12 m/s, showing the histograms of rotor speed, pitch angle, and electric power. Two constant de-rate commands were considered: 100% (equivalent to rated power) and 60% . As there was enough wind velocity to reach rated power, for both wind directions, the required levels of electric power could be accomplished, while the rotor speed oscillated around its rated value in both cases. Regarding the pitch angle, higher values were observed at wind direction 250 • than at 147 • because in the latter case, WT3 was at the wake of WT4 with less available power. Furthermore, higher pitch values were adopted for the 60% de-rate command than for the 100%, for the same wind direction (e.g., 250 • ). Figure 10. WT3 operation. Inflow condition of the simulations: Wind direction 250 • and wind direction 147 • (without and with wake interaction respectively) and wind speed above rated. Two constant de-rated signals applied, 100% (top) and 60% (bottom). Figure 11 depicts the total wind farm power output, as well from the individual wind turbines, comparing three de-rate levels, for wind direction 147 • and wind speed of 9 m/s at hub height at the inlet. The de-rate commands were equally distributed to each wind turbine, regardless of their capacity to reach the desired power output. Only in the case of DR = 50% was the required level reached most of the time. In the other simulations, the power output was not reached mainly because WT3 did not reach its required power as it was strongly affected by WT4 wake. These results show the need to implement a global wind farm controller, with the purpose of assigning de-rate commands based on the available power of each wind turbine, particularly when there are wind turbines operating in the wake of others. Different types of global controllers were studied in [54,55], where they evaluate the effect on the total power output, considering an open-loop or a closed-loop global controller, and it is planned to implement those in Caffa3d. Figure 11. Power production of the wind farm (top) and of each individual wind turbine. Inflow condition of the simulation: Wind direction 147 • (wake interaction between wind turbines) and wind speed below rated. Three constant de-rated signals applied on each wind turbine.

Conclusions and Future Work
A Large Eddy Simulation framework with the Actuator Line Model was used to simulate the operation of a 7.7 MW onshore wind farm. The simulations considered different ABL wind profiles at the inlet, with hub height velocities below and above the rated one, and also considered two different wind directions, one of them with strong interaction between wakes and wind turbines. A closed-loop collective-pitch controller and a torque controller were implemented in the CFD code based on the one presented in [52]. The results of the simulations were compared to high frequency data (1 Hz) from the wind farm SCADA system, obtaining good agreement, both regarding the mean values and the temporal evolution of the signals. The ALM was shown to be capable of simulating with a moderately computational cost, compared with a FRS, the wind flow through wind turbines and wind farms. The ALM proved to be a complementary tool to the models currently applied by the industry, such as the BEM method, as it was envisioned by [7].
The response of the controller subject to different de-rate values, both constant and time-varying, was evaluated. We showed that when there was enough wind power in the flow, the individual wind turbines could follow the required signal, but when it came to the total wind farm power output, it failed to accomplish what was required because of the interaction between wakes and turbines. These results show the necessity of implementing a global controller, which takes into account the available power of each individual wind turbine, which we plan to work on in the near future. In addition to this, work is currently underway to couple the mesoscale model Weather Research and Forecasting (WRF) with the Caffa3d code so as to complement the available microscale measurements, looking to simulate an actual event. A structural model is also planned to be coupled to the LES-ALM approach presented.
Besides, the flow solver, including the wind turbine modules, is being expanded to use GPU computing platforms, as presented in [59], employing a dual CUDA/OpenCL syntax on top of the coarse MPI parallelization. Tests performed have shown speedups of up to 30x with respect to the same simulation executed using CPU platforms. Funding: This research was funded by the Agencia Nacional de Investigación e Innovación -Fondo Sectorial de Energía (ANII FSE).

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.