Large-Eddy Simulation of Waked Turbines in a Scaled Wind Farm Facility

The aim of this paper is to present the numerical simulation of waked scaled wind turbines operating in a boundary layer wind tunnel. The simulation uses a LES-lifting-line numerical model. An immersed boundary method in conjunction with an adequate wall model is used to represent the effects of both the wind turbine nacelle and tower, which are shown to have a considerable effect on the wake behavior. Multi-airfoil data calibrated at different Reynolds numbers are used to account for the lift and drag characteristics at the low and varying Reynolds conditions encountered in the experiments. The present study focuses on low turbulence inflow conditions and inflow non-uniformity due to wind tunnel characteristics, while higher turbulence conditions are considered in a separate study. The numerical model is validated by using experimental data obtained during test campaigns conducted with the scaled wind farm facility. The simulation and experimental results are compared in terms of power capture, rotor thrust, downstream velocity profiles and turbulence intensity.


Introduction
Wind farms are collections of wind turbines, often operating in close proximity of one another. During operation, several major phenomena take place within a wind farm: there is an interaction between the atmospheric boundary layer and each individual wind turbine (as well as a more global interaction between the boundary layer and the whole wind farm), and an interaction among upstream and downstream wind turbines through their wakes. The wake produced by the upstream wind turbines, coupled with the atmospheric flow, has a significant influence on the performance of downstream operating machines. In fact, the wake-affected inflow results in lower power output and increased loading for the downstream wind turbines when compared to an isolated condition. A number of control strategies are currently being investigated to optimize the operation of wind farms, including power derating and wake deflection [10].
Our group has developed scaled wind turbine models, operated in a large boundary layer wind tunnel [4]. The machines are actively controlled, and optionally aeroelastically scaled, and are designed to generate realistic wakes notwithstanding the low Reynolds numbers at which they operate. This unique facility, which has been previously used within other research projects [5][6][7], enables experiments to study aeroservoelasticity, turbine wakes, machine-to- machine interactions, and wind farm control for power maximization and load mitigation. The facility is highly instrumented, allowing for the collection of a wide range of high quality data, both regarding flow conditions and the mechanical response of the machines. A parallel on-going research effort of our group aims at developing a digital copy of this experimental facility, including the wind turbines and the wind tunnel itself. State-of-theart numerical methods are used to model the wind tunnel generated flow, the wind turbines (including blades, nacelle, tower), and their shed wakes interacting with the other machines and the wind tunnel. Initial efforts towards the validation of this simulation tool for the single wind turbine case were described in previous work [21]. This paper describes recent improvements in the capabilities of the numerical model for the single wind turbine case, and it also considers multiple waked setups in a number of operating conditions, including the effects of changes in axial induction (power derating).

Numerical model 2.1. CFD formulation
The current simulation model is developed within SOWFA [8], a simulation tool based on a standard incompressible solver in the OpenFOAM repository. The current version of the code uses an actuator line method embedded in a large-eddy simulation (LES) environment, coupled with the aeroservoelastic simulator FAST [13]. An immersed boundary (IB) formulation [11] is used to model the wind turbine nacelle and tower. A blended connective scheme algorithm, coupled with the Gamma [12] and center differencing schemes, is also implemented in the solver to minimize numerical diffusion (Gamma is used at near wake, while center differencing scheme is used at far wake). The employed numerical methods are described in detail in Ref. [22].

Immersed boundary method
An IB method is employed to model the turbine geometry, as this approach does not require a complex body-conformal mesh, which usually implies a tedious process and significant effort in the present structured-grid context [15]. Furthermore, the immersed boundary method can usually achieve an adequate local resolution with minimum overall mesh density, and eliminate the risk of obtaining deteriorated local cells. The adopted IB approach uses a direct imposition of the boundary conditions (discrete forcing approach). To introduce the characteristics of this method and explain its advantages over the alternative indirect imposition method, the principles of both are briefly reviewed next. The present formulation is derived based on Mittal et al. [15] and Bandringa's work [1]. The incompressible governing equations are given as where − → u is the fluid velocity, p the pressure, ρ the density and µ the dynamic viscosity. The governing equations are first discretized without considering the immersed boundary and are written as where U * = ( − → u , p) specifically indicates the solution of the velocity field without considering the effects of the immersed boundary, while L is the operator representing the discretized Navier-Stokes equations. Next, the immersed boundary condition (Γ b ) is imposed on the flow-field U = U Γ , resulting into Eq. (2) being reformulated as Here [I] is the identity matrix and δ is the Dirac delta function, which allows for the formula to switch between Eq. (2) and U = U Γ . However, the cell center of a specified volume mesh location (vector), noted − → x k , in general does not coincide with the immersed boundary point, − → x ib . The two numerical methods -direct and indirect -differ in the way they address this problem. The indirect method tries to replace δ with a smear function d, which writes The forcing term resulting from this smooth distribution implies that sharp edges of the geometry cannot be exactly represented. Therefore, the forcing term will smear in space, which is an undesirable effect when one wants to resolve the surface boundary layer of the geometry at high Reynolds numbers.
On the other hand, the direct imposition approach preserves the sharpness of the geometry. No smeared force will be generated in the neighboring cells, but instead an interpolation technique is used to treat the IB cell center using its neighboring cells. The current solver uses a two-dimensional polynomial function with second order accuracy [11]. The Dirichlet immersed boundary condition imposed at location k is formulated as where all the coefficients C 0 , C 1 , C 2 , C 3 and C 4 are computed using weighted least squares on the extended stencil [11]. x k and x ib represent x coordinate (scalar) of the specified points. With this interpolation technique, the direct imposition method is able to locate the boundary surface, identify each individual solid and fluid cell, and assign proper boundary conditions to the fluid-solid interface. This way, the sharpness of the geometry surface is preserved and no extra stability constraints are introduced. With the sharp representation of the IB, one can directly impose boundary conditions and a wall model on the surface boundary, which typically results in an improved quality of the solution for high Reynolds number viscous flows.

Computational setup
The computational setup is similar to the one described in Ref. [22], and it is shown in figure 1. Zone 1 represents the background mesh, with isotropic grid spacing of ∆x = ∆y = ∆z = 0.08 m. Zones 2 and 3 represent progressively denser grids, the latter having a cell size of 0.01 m. The wind turbines, whose diameter D is 1.1 m, are aligned with each other in the streamwise direction (i.e. there is no lateral offset), and longitudinally spaced by 4.1 rotor diameter. Two layouts are considered, comprising either two or three aligned wind turbines. During the experiments, each wind turbine operated under its own pitch and torque control system, and therefore it adjusted automatically to the local wind conditions. On the other hand, the simulated wind turbine rotors spun at the constant angular velocity measured in the experiments. It should be noticed that this has also the effect of limiting the impact of errors in the inflow conditions for downstream wind turbines, which clearly deviate from the real ones on account of inaccuracies in the modeling of the wakes. The closed-loop operation of the simulated wind turbines would amplify such errors, and it will be the subject of specific investigations once the validation and tuning of the present simulation environment has been completed.
A static velocity inflow map, measured experimentally using scanning LiDARs [20], is imposed at the inlet of the computational domain to represent the non-uniform inflow generated by the wind tunnel. This inflow map offers a more realistic input compared to a uniform inlet velocity, as described more in detail in Ref. [22]. In fact, it was observed that matching the inflow non-uniformity ultimately increases the precision of the numerical simulations. In particular, the simulated wake profiles more closely match the experimental data. This is a nontrivial   Table 1 reports the power (P ) and thrust (T ) for a single wind turbine setup at four different power ratings, namely 1.0, 0.975, 0.95, 0.925. The machine was operating in region II (partial load condition) at a wind speed of 5.97 m/sec measured at hub height, 3D in front of the rotor disk, while the air density was approximately 1.19 kg/m 3 . The wind speed was employed to normalize the velocities measured within the wake. Due to the inflow non-uniformity, the rotor-disk average wind speed was however slightly lower, i.e. approximately 5.83 m/sec [22]. Derating is accomplished by holding the rotor angular velocity constant above the wind speed that makes the power production exceeding its derated level at optimal pitch setting, and then using the pitch controller to maintain the power production as the wind speed increases. This means that reducing the demanded power requires either pitching the blades to feather while also slightly decreasing the rotor angular velocity. The model includes the presence of the nacelle and tower, although this has only a very marginal effect on the rotor performance.
Results show that power and thrust are in general in reasonable agreement. However, for increasing derating the error in the power increases, while the error in the thrust decreases. This is probably due to the accuracy of the airfoil polars. In the present simulations, the airfoil lift and drag are computed at each instant of time based on the local angle of attack by interpolating tabulated values at Reynolds 75000 and 90000. Future efforts will try to improve the quality of the polars by identification techniques [3], as these have clearly a crucial importance in a LES-lifting-like formulation.

Wake properties
The quantitative verification of the wake characteristics is essential for the validation of simulation tools, as the matching of rotor quantities does clearly not guarantee an accurate modeling of the wake. For the non-derated case, figure 2 shows the time-averaged  Circles are used to report the experimental values, obtained by triple hot-wire probes [20]. First of all, it is important to note that the wake deficit measured immediately downstream of the rotor agrees with 1D momentum theory expectations, since the BEM-based axial induction factor of the rotor for the non-derated case is approximately 0.31. Moreover, CASE 1 performs better than CASE 2 in terms of velocity profile for all three downstream positions. This highlights the importance of including nacelle and tower in the simulation models (as already noted by other authors [18]), especially for the present scaled wind turbine for which these two components are relatively bigger than in a multi-MW wind turbine. Indeed, the sum of the frontal area of the nacelle and of the portion of tower located within the rotor swept area is 0.037A, with A = π D 2 4 , while it is 0.023A for the 5MW NREL wind turbine [14]. Although this parameter is larger for the G1, we expect that the effects of nacelle and tower on wake evolution might not be negligible even for typical multi-MW wind turbines.
The mean-over-the-rotor velocity difference ∆U between simulation and experiment is -0.4%, 2.9% and 2.9% for CASE 1 at 3D, 4D and 7D, respectively, while for CASE 2 it is 6.9%, 17.5% and 6.4% at the same respective positions. Moreover, the wake profiles of CASE 1 are closely correlated with the experimental curves for all three downstream positions. Matching is also good at 3D, a position in the near wake. Reference [22] also shows results at 0.56D, a position far closer to the rotor, where results are reasonable but not as accurate as here. While it is still very challenging to achieve high accuracy very close to the rotor (LES-lifting-line doesn't allow fully capturing the near wake behavior), the present results show that relatively good results can be obtained further downstream, where mixing generates a more elementary, smooth and quasi-Gaussian shape of the wake speed deficit.
An accurate modeling of turbulence intensity at the near wake by LES is still difficult to achieve, with some good quality results shown by Ref. [17] and [16]. For the present work, figure 2 in the bottom row reports the time-averaged TI profiles. The mean-over-the-rotor-disk difference ∆T I between simulation and experiment is -18.8%, -17.3% and -25.9% for CASE 1 at 3D, 4D and 7D, respectively, while ∆T I for CASE 2 is -52.1%, -41.2% and -41.3% at the same respective positions. Clearly, the accuracy of the mean TI for CASE 1 is significantly superior to CASE 2, despite still having a relatively large margin of error with respect to the experimental results. Here again the inclusion of nacelle and tower in the model proves to be of great importance.
It could appear surprising that a good match between experimental and numerical wake deficits is observed, while the prediction of the turbulence intensity is unsatisfactory. However, the following considerations could potentially explain this behaviour. First, unlike Pesudo- spectral and Energy-conserving discretization method, the adopted finite-volume method cannot guarantee zero numerical dissipation [2], which can only be limited by employing a fine spatial and time grids, as well as by using dedicated differencing schemes. As described in section 2.1, the low-diffusion Gamma scheme is used in the present LES framework. The results shown in this section were obtained by manually tuning the parameter β m , a weighing factor for the Gamma Convection bounded scheme that controls the amount of numerical dissipation [12]. By decreasing β m , one may observe the appearance of numerical oscillations in proximity of the immersed boundaries.These may increase the tendency of the simulation to crash and induce artificial fluctuations that contaminate the velocity profile whose accuracy, in turn, drops significantly. Increasing β m , on the other hand, increases the upwind-like behavior of the scheme, which results in more turbulent kinetic energy being dissipated, while leaving the velocity profile largely unaffected. Indeed, the numerical dissipation is mainly affecting second order quantities, like TI, while first order quantities, such as velocity, are less influenced by the intrinsic numerical dissipation. Second, in this paper a spatially non-uniform but temporally steady inflow map was used as inlet condition for the simulations. Such inflow differs from experimental inlet flow, which is characterized by a 2% of turbulence intensity. Such a difference can also be observed in figure 2, where the simulated TI outside the wake is almost null. This implies that the experimental wake recovery should be slightly faster than the simulated one, which in turn means that the experimental wake deficit should be slightly lower than the simulated one. On the other hand, the intrinsic numerical diffusion of the adopted numerical method makes the wake to recover slightly faster than what is observed in reality. This means that, on one side, the lower simulated TI implies a stronger wake deficit, while the intrinsic numerical diffusion implies a lower wake  figure 1). As shown earlier in section 4.1, the wake profile and velocity magnitude for the single wind turbine case appears to be of a good quality. Therefore, it is expected that the inflow conditions for the donwstream machine in the two wind turbine case will also be accurate. Table 2 reports the power P and thrust T for wind turbine No.2. In the non-derated condition for CASE 1, there is a reasonable agreement for P with respect to the experimental value, with only an about 4% difference. Such accuracy for power is non-trivial to obtain, as its absolute value is much lower than that of wind turbine No.1 and there is also a strong sensitivity of this result to even small variations in the wake shed by the upstream machine. On the other hand, CASE 2 shows a much larger discrepancy in power, which is about 20% larger than the experimental value. In fact, neglecting the nacelle and tower from the model will cause a substantial change in the downstream flow mixing, resulting in the strong differences in the velocity profile and TI shown in figure 2. The table also illustrates that the accuracy of the results decreases with increasing derating of wind turbine No.1. The reason of this behavior is still under investigation. On the other hand, the results show that the model is able to appreciate the expected power increase in wind turbine No.2 that follows from derating wind turbine No.1. Figure 3 shows the wake properties in terms of time-averaged velocity profile and TI downstream of the two aligned wind turbines. Similarly to the wake analysis for the single wind turbine configuration, results for CASE 1 are well correlated with the experimental values. The mean velocity difference between simulation and experiment across the rotor disk, ∆U, is 1.5%, 2.7% and 2.6% for CASE 1 respectively at 8D, 10D and 11D downstream of wind turbine No.1, while ∆U for CASE 2 is -8.1%, -10.1% and -8.4% at the same positions. Unlike the single wind turbine setup, the mean wind speeds for CASE 1 are consistently greater than for CASE 2, with the latter showing a larger wake deficit than the former at all three downstream positions. This might be due to a lack of flow mixing when nacelle and tower are not present in the model, and to the larger power capture of wind turbine No.2. The good comparison of wake properties at 11D (slight higher wake recovery), a distance  The behavior of TI for the double wind turbine configuration is well predicted in terms of its shape and mean value across the rotor. Indeed, the effects produced by the intrinsic numerical diffusion become less dominant as the upcoming flow becomes turbulent, as shown by Troldborg et al. [19]. Since the upstream wind turbine adds substantial turbulence to the flow impinging on the second wind turbine, the results for the double wind turbine setup appear to be better than in the single wind turbine case by at least 10% on average.

Wake properties
Moreover, by comparing TI for CASE 1 and CASE 2 at 8D, a position that is located 3.9D downstream of wind turbine No.2, it appears that neglecting nacelle and tower from the model dramatically underestimates TI by over 40%. This provides an explanation for the consistent over-estimation of the velocity deficit seen for CASE 2, as discussed earlier in this section.

4.3.
Three aligned wind turbines 4.3.1. Integral rotor quantities The three wind turbine configuration is considered only in the non-derated condition. Table 3 reports power capture and thrust for all three operating wind turbines. The power capture for the two front wind turbines shows a slight difference with respect to the values in Table 1 and 2. This is due to a slight different inflow speed in the wind tunnel during the two sets of experiments. In fact, since each wind turbine was operated by its own closed-loop controller, a variation of inflow speed will result in the wind turbines operating at a different set point. The power capture of wind turbine No.3 for CASE 1 shows a 17% difference with respect to the experimental value. This can be attributed to a lack of background Similarly to the two aligned wind turbines cases, we observed that the accuracy of the results decreases as the derating of the upstream wind turbine increases.  Figure 4 shows normalized instantaneous U x and TI fields at a streamwise vertical plane for CASE 1 and CASE 2. The plots reveals the dramatic difference in the flow evolution which results from including or not nacelle and tower in the model. In fact, nacelle and tower generate a substantial amount of turbulence in the field, which in turn induces an enhanced downstream flow mixing. Plots of Reynolds stress, vorticity and turbulent kinetic energy all demonstrate a similar dependence on the presence of nacelle and tower. At the same time, these model features also cause considerable changes in the air flow velocity in the wakes, creating large discrepancies in power output and loading. Furthermore, the TI plots also indicate that turbulence starts to accumulate after approximately 2D downstream of wind turbine No.2 and 3, whereas it is relatively low immediately downstream of the rotors. As expected, TI is higher downstream of the blade tip and root regions than for the rest of the blade span. The rotor appears to be initially damping, just immediately downstream from it, the incoming TI, which then bursts up again further downstream.

Conclusion and outlook
In this paper, the results from a LES-lifting-line simulation formulation were discussed with reference to the experimental measurements obtained on scaled wind turbine models operating in a wind tunnel. Three different setups were considered: a single isolated machine, as well as double and triple aligned and fully waked configurations.
From the comparison of power capture, thrust, wake profile and other flow quantities, it is concluded that the results of the simulation model are in general in a good agreement with the measured data. The quality of the results, however, appears to be decreasing for increasing derating of the front machine, an issue still needs to be properly investigated.
The paper also considered the effects of the inclusion of nacelle and tower in the model. In all cases, results always consistently indicate the crucial importance of modeling the effects of nacelle and tower, as their absence has a strong effect on the mixing of the flow. The nacelle and tower of the scaled models are in general larger than the ones of a real wind turbine, because of the need to house all necessary mechanical and electronic components, including pitch, torque and yaw actuators and various sensors. Because of this, the effects of neglecting nacelle and tower from the model might be more pronounced in this than in a full scale case. Similar conclusions on the importance of a complete model have also been reported by other authors [18]. A parallel study is considering a similar validation activity for cases of higher ambient turbulence. Turbulence generating devices were used in the wind tunnel to generate a turbulent sheared profile, in which the wind turbines were operated in different configurations and at different operating points, as done in the present work. A large-eddy simulation of the turbulence generation will be used to produce a turbulent wind field (same way as in the experiment), which will then be used as input for the subsequent LES-lifting-line simulations of the interacting wind turbines. The study of these conditions, in addition to the low turbulence case considered here, will be used to further the understanding of the present simulation framework, with the final goal of obtaining a digital model of the scaled wind farm facility with predictive capabilities.