Numerical and Experimental Study of Wake Redirection Techniques in a Boundary Layer Wind Tunnel

The aim of the present paper is to validate a wind farm LES framework in the context of two distinct wake redirection techniques: yaw misalignment and individual cyclic pitch control. A test campaign was conducted using scaled wind turbine models in a boundary layer wind tunnel, where both particle image velocimetry and hot-wire thermo anemometers were used to obtain high quality measurements of the downstream flow. A LiDAR system was also employed to determine the non-uniformity of the inflow velocity field. A high-fidelity large-eddy simulation lifting-line model was used to simulate the aerodynamic behavior of the system, including the geometry of the wind turbine nacelle and tower. A tuning-free Lagrangian scale-dependent dynamic approach was adopted to improve the sub-grid scale modeling. Comparisons with experimental measurements are used to systematically validate the simulations. The LES results are in good agreement with the PIV and hot-wire data in terms of time-averaged wake profiles, turbulence intensity and Reynolds shear stresses. Discrepancies are also highlighted, to guide future improvements.


Introduction
Wake redirection is currently one of the most promising wind farm control methodologies, which may lead to the improved power capture and reduced structural loading. The current research on wake redirection is very active, covering high-fidelity numerical simulations [14], scaled experiments in the wind tunnel [10], direct measurements in the field [1], reduced order models [2] and control methods [9]. In a previous paper [27], we used a three-pronged approach based on the combination of an analytical formulation, a LES-lifting-line CFD method [26] and wind tunnel scaled models to investigate two wake redirection methods, yaw misalignment (noted YM in this paper) and cyclic pitch control (noted CyPC). The analysis showed that, as already noticed by other authors, YM is more effective than CyPC in displacing the wake laterally. CyPC, however, was shown to have some effect on the speed of recovery of the wake, although the practical use of this concept is put in serious doubt by the large loads that are generated on the wind turbine.
In the present paper, we use additional wind tunnel measurements to further characterize the behavior of the CFD approach in the context of these two wake redirection techniques. Our long term goal is to develop the ability to accurately simulate our wind tunnel experiments by a validated numerical environment, and to upscale the validation to full scale. The present paper considers both integral rotor quantities such as power and thrust, as well as local flow measurements obtained by particle image velocimetry (PIV) and hot-wire probes. Scanning LiDARs are used to create a map of the wind tunnel inlet, to account for some spatial variability of the generated flow conditions, which is useful for a more accurate verification of the numerical results. The analysis highlights the general good qualitative agreement between experiments and simulations. The results also show that some local improvements are still necessary, highlighting some weakness of the numerical approach that motivate further refinements.

Simulation model 2.1. CFD formulation
The simulation model is implemented in the SOWFA code [12], a CFD simulation tool based on an incompressible solver within the OpenFOAM repository, employing the Boussinesq approximation to include buoyancy-driven flow. Large scale transient simulations are conducted using the PIMPLE time marching algorithm. The grid-collocated formulation computes quantities at cell centers. The interpolation of Rhie and Chow [24] is employed to remove numerical oscillations caused by pressure-velocity decoupling induced by grid collocation. The time marching scheme uses a second order backward implicit discretization. The momentum and potential temperature equations use a diagonal incomplete preconditioned bi-conjugate gradient linear solver, while pressure is solved using a geometric agglomerated algebraic multi-grid solver. Parallel computing is based on the Message Passing Interface (MPI) [13]. In addition, SOWFA uses the lifting-line method embedded in a large-eddy simulation (LES) environment, coupled with the aeroservoelastic simulator FAST [13].
Apart from the numerical methods used for wind farm simulations, other techniques are employed to perform more robust and accurate large-eddy simulations of the scaled wind farm facility. The immersed boundary formulation is used for the modeling of geometries, like nacelle and tower, within the control volume. The convection vector-field scheme uses a deferred correction Gamma-bounded interpolation method, which belongs to the family of categorized normalized variable diagram (NVD) schemes [18]. This Gamma scheme contains a pre-specified constant, β m , which allows one to limit the numerical diffusion and to minimize numerical dispersion. In general, a larger value of β m implies a lower dispersion and a higher diffusion, vice versa. Therefore, β m =0.35 is employed within the near wake to stabilize the simulations, since actuator line body forces and immersed boundary possibly increment the numerical stability issues, while β m =0.05 is used for the far wake to minimize numerical diffusion. This algorithm represents an important improvement with respect to normal unbounded and non-flux-limited central differencing schemes. In fact, it is capable of minimizing the numerical dispersion otherwise often observed at large yaw misalignment angles, where the wake is oblique to the grid. Moreover, this approach performs better than other low diffusive schemes such as the one of Vanleer in terms of limiting the numerical diffusion in the wake region, where small-scale eddies are dissipated due to such flux limiters [16]. A tuning-free Lagrangian scale-dependent dynamic approach is used for the Smagorinsky sub-grid scale modeling [23]. This approach improves the modeling of inhomogeneous flow fields such as in the case of wind turbine wakes, where first and second-order flow quantities vary significantly within the domain.

Actuator line method
The LES formulation is augmented with an actuator line method. The complete geometry of the blades is not resolved in the present LES formulation, due to its excessive computational costs. Instead, blades are represented by a given number of points on an actuator line. Aerodynamic forces are calculated based on the local angle of attack by means of look-up tables storing the local airfoil lift and drag coefficients. The resulting forces acting at each blade section are then smoothly mapped to a body force field by a Gaussian projection characterized by a width parameter ε. The determination of the ε parameters is dependent on several factors. According to Ref. [13], ε should account for cell size, chord width and other geometrical parameters, but cannot be lower than a minimum value to avoid spurious oscillations of the velocity field. A comprehensive study of the effect of the projection width on the simulation quality is reported in Ref. [20].

Immersed boundary method
An immersed boundary (IB) formulation [17] is used to model the wind turbine nacelle and tower. The IB method is employed to avoid the use of surface conforming meshes to resolve the shape of the body [22]. The present IB approach, based on a discrete forcing method, uses a direct imposition of the boundary conditions. This way, the sharpness of the body shape is preserved and no extra stability constraints are introduced. In fact, with this approach boundary conditions and wall models can be directly imposed on the IB surfaces, yielding a significant improved solution quality for higher Reynolds number viscous flows [3], such as the present case of wind turbines operating in an atmospheric boundary layer. Details of the formulation are described in Ref. [28].

Experimental methods
Experiments were conducted in a 36 m × 16.7 m × 3.84 m boundary layer wind tunnel at Politecnico di Milano (POLIMI) [7]. Tests were performed with a scaled wind turbine model, shown in figure 1, together with the reference frame adopted to express the flow velocity components, which has a rotor diameter (D) of 1.1 m and a rated rotor speed of 850 RPM, and which has already been used within other research projects [9][10][11]. The model is characterized by a realistic aerodynamic performance, both at the airfoil and rotor levels, and its wake is characterized by shape, deficit and recovery that are in good accordance with those of full-scale machines. Moreover, the model features active individual pitch, torque and yaw control that, together with a comprehensive onboard sensorization of the machine (including measures of shaft and tower loads), enables the testing of modern control strategies similar to the ones described in [5] and references therein. The wind tunnel can be operated with different configuration of spires and pyramids at the inlet of the test section, to realize varying vertical shear and turbulence intensity conditions. Hot-wire probes, LiDAR (Light Detection And Ranging) and PIV (Particle Image Velocimetry) are used to perform high quality measurements of aerodynamic quantities, which will then be used for the validation of the LES-lifting-line model.

Inflow measurements
Previous experimental testing in the POLIMI facility has revealed that the inflow generated by the wind tunnel is non-uniformly distributed, with variations in flow velocity of up to 6%. The inhomogeneity of the flow is non-negligible, and it will influence the operation of the wind turbine models and of the associated flow-field. To address this problem, instead of imposing a uniform inflow for the inlet boundary condition, a non-uniform inflow was used to more accurately represent the incoming flow generated in the wind tunnel. Properly capturing the inflow velocity requires the full cross-sectional wind map to be recorded upstream of the wind turbines, rather than sampling a single point (or several scattered points in space) at hub height. The measured experimental wind map can then be implemented as the inlet boundary condition to the CFD simulation.
In order to obtain a reliable high resolution wind map, two synchronized scanning LiDARs were employed. LiDARs are instruments equipped with a steering laser beam that fires rapid pulses and, by measuring the frequency shift of the backscattered light from the aerosols present in the air, are able to record the line-of-sight component of the wind speed. In this work, LiDAR data was obtained from the collaborative efforts of ForWind-Oldenburg, TUM, Technical University of Denmark and POLIMI [25]. Figure 2 shows the time-averaged velocity scanned inflow map. The figure clearly shows the lack of uniformity of the inflow, probably due to the presence of a discrete number of fans and of the structural arrangement of the tunnel just upstream of the test section. The black circle indicates the location of the rotor, which appears to be placed in correspondence with a local peak of the flow velocity. wake. In order to achieve a higher spatial resolution of the velocity field, the measurement area was divided into several smaller windows with slight overlapping areas between them. A rapid scanning of the entire measurement area was achieved by the use of an automated traversing system, moving both the laser and the cameras. In particular, the two cameras, equipped with a 1952 × 1112 pixel array, were connected to a metallic arm and mounted on a 4 m × 2 m double axis traversing system. The Nd:Yag double pulsed laser, with a 200 mJ output energy, was also mounted on a traversing system to move the laser sheet simultaneously with the measurement window. The measuring windows were divided into 32 × 32 pixel interpolation areas, which resulted in an approximately 15 mm spatial resolution. For each measuring window, 200 pairs of images were acquired (per camera) without any phase lock (time-averaged flow field). Additional details concerning the PIV instrumentation are given in [29] 4. Computational setup The computational domain, as shown in figure 3, has a height of 3.84 m height, equal to the height of the wind tunnel test section. The inlet is located 3.3 m upstream of the wind turbine, in the same plane where the LiDAR inflow map is measured. The streamwise length is 10D. The lateral dimension is reduced to 4 m, to limit the domain size and hence the computational cost. Multiple simulation tests have shown that the reduction of the lateral dimension has no noticeable effect on the results. After determining the size of the computational domain, grid resolution requirements are investigated through a mesh independence study. Three levels of refinement are used: zone 1 defines a cube-shaped background mesh (∆x = ∆y = ∆z) with a cell size of 0.08 m, zone 2 is the first level of refinement with a cell size of 0.04 m, followed by zone 3 that uses a size of 0.01 m. The inflow velocity measured by LiDARs is used at the domain inlet, to capture the lack of uniformity of the wind tunnel generated flow. Apart from this, the inflow of the numerical simulation is steady and non-turbulent. This creates a difference between the numerical and experimental flow conditions, because the latter has a turbulence intensity equal to about 2%. This should be kept in mind when comparing numerical and experimental results, because the different ambient turbulence conditions may have some effects.
Slip boundary condition is applied at the left and right walls, while a wall model equivalent to a no-slip condition is used for the floor and ceiling of the tunnel. A wall model is also used for the immersed boundary surface that models the nacelle and tower. Neumann type boundary condition for the pressure at all walls is also imposed. A modified version of the mesh partitioning algorithm was developed, which tries to allocate the same number of cells to all processors without cutting the mesh at the immersed boundary surfaces. A fixed time-step of 0.0002 s is used for baseline and two CyPC setups, while 0.0003 s is used for YM setup. The time-step equates to a maximum Courant number C max of 0.17 and a displacement of 0.99 cells per time-step at the tip of the blades for the baseline case. In all cases, the machines are operating in the partial load region at an averaged wind speed of 5.83 m/sec across the rotor, while the rotating speed varies from 780 to 840 RPM due to different redirection control methods. The Initial blade pitch were kept to be constant for all 4 setups. The ratio ε ∆x between the Gaussian width and the cell size has been kept equal to 2.2 for the baseline and two CyPC simulations, while it has increased to 3.4 for 20 deg yaw misalignment case. Most simulations were performed on an Intel Haswell-EP computing cluster. With these computational resources, 5.2 s of physical time (≈26,000 time steps) require about 13,440 cpu hours. The solution converges to steady state after about 4 s of physical time.

Power and thrust
Four operating conditions are considered: a baseline condition (constant pitch setting, no misalignment), cyclic pitching with two different phase angles, and a yaw misalignment of 20 deg. For CyPC, the blade is pitched according to θ i = θ 0 + θ c cos(ψ i + γ), where θ 0 is the collective pitch constant, θ c the 1P pitch amplitude, ψ i the blade azimuth angle (counter-clockwise looking upstream, and null when the blade is pointing vertically up along the tower), and γ is the phase angle (clockwise looking upstream). Two values for γ are considered: 52 and 270 deg. The wind tunnel was operated without spires at the test chamber inlet, resulting in low-turbulence (≈2%) flow conditions. Table 1 shows the wind turbine power and thrust at the four operating conditions. A good agreement between the simulations and experiments is observed for the baseline case, with differences of 3% and 6% for power and thrust, respectively. The match is less accurate for the CyPC cases. For both cases, power is underestimated by the numerical simulations, while thrust is underestimated in one case and overestimated in the other. Several possible effects may help explain the larger errors in the CyPC cases. First, unsteady airfoil aerodynamic effects (including dynamic stall) are completely neglected, as the present lifting-line formulation only uses look-up-tables of steady lift and drag coefficients. This could be improved by using the suitable unsteady aerodynamics model for the airfoils. Some of these models, such as the Beddoes-Leishman [21], are available within the current LES framework. However, this model requires the definition of several airfoil-dependent governing parameters that were unknown to us at the time when simulations were performed, and would need to be specifically calibrated. Second, the CyPC experiments were conducted at a lower RPM than the  [10], are tabulated for two values of the Reynolds number (namely, Re=90000 and Re=75000), and they are interpolated during the simulation based on the local instantaneous conditions [21]. However, the lower Reynolds coefficients may be less accurate than the higher ones, which might induce a larger discrepancy in the results in this case. This might be improved by tuning the aerodynamic coefficients, as suggested in Ref. [6]. The same table shows also good results for the misalignment case. However, various numerical experiments showed a significant sensitivity of the results to the Gaussian width parameter ε. In fact, the results of the table were obtained with a parameter ε 54% higher than the one used for the baseline case, as addressed in section 4. This problem deserves further careful investigation. Here again, a possible reason might be the lack of modeling of unsteady aerodynamic effects at the level of the airfoils, effects that increase for increasing misalignment. Clearly, finding a specific value of the Gaussian width parameter to correct for this deficiency of the method is not a sensible option, and it is certainly not suggested as a general solution. The analysis of the U x component at 0.56D downstream of the rotor is particularly interesting, as this is a position where few results have been previously reported. The images show that the use of CyPC has a strong effect on the wake structure, leading to marked unsymmetrical shapes. Measurements are missing from two areas left and right of the rotor disk, where, due to the close proximity of the measuring plane with the wind turbine, part of the nacelle (which is of a white color) was on the background, leading to bad correlation between the images. A comparison between experimental and numerical results shows that there is in general a good qualitative agreement, and that the main distortion effects caused by CyPC are reasonably captured. The difference ∆U i between experiment and simulation velocities was evaluated at several locations within the rotor swept area (black dot circles in figure 4). Based on these quantities, the average velocity difference ∆U and its root mean square ∆U RMS were calculated. ∆U is equal to -1.36%, 2.69% and 5.05% for the baseline, CyPC γ=52 and 270 deg cases, respectively, while ∆U RMS is 0.43 m/sec, 0.79 m/sec and 0.72 m/sec for the same three cases. Positive value of ∆U means that the average simulated velocity is higher. Figure 5 shows the wake plots in the far wake region at 6D. ∆U in this case is 4.9%, 9.4% and 4.9% for the baseline, CyPC γ=52 and 270 deg cases, respectively. ∆U RMS is 0.33 m/sec, 0.46 m/sec and 0.44 m/sec for the same three cases. ∆U RMS are slightly improved at 6D compared to the near wake for all three cases, while ∆U at 6D show in general higher wake recovery. Lower ∆U RMS implies a higher correlation between LES and PIV wake profile. In fact, compared with the near wake region, the flow in the far wake is further mixed and less structured, exhibiting a velocity profile that is typically well approximated by a single Gaussian. Lastly, it should be remarked that CyPC leads to a faster recovery of the wake than in the baseline case, as already noticed in Ref. [27]. In principle, this could be of interest for wind farm control, although the very large loads exerted on the rotor by the use of CyPC severely limits the possible applicability of this control concept.
Next, the hub-height velocity profile is considered. Figure 6 reports the time-averaged streamwise velocity profiles at hub height for the three considered operating conditions. In the near wake (0.56D), matching for the baseline case is reasonable, except the region close to the hub. Here, the effects of the nacelle and tower are not captured accurately. A visualization  of the numerical results seems to indicate a strong tower shedding. Accuracy might be affected by a lack of background turbulence in the simulations. In fact, as previously explained, the wind tunnel flow is characterized by a turbulence intensity of 2%, while the numerical simulations are conducted in uniform non-turbulent inlet conditions. For the CyPC case, matching is much poorer, with significant differences both close to the hub and especially at the tip of the low speed region. As previously noted, differences might be due to limits of the lifting-line formulation, and/or because of a lack of grid resolution. The situation is improved for the 6D case, as expected. Here good results are obtained for the baseline case, while small discrepancies are visible for the CyPC conditions. Note that the wind speed velocities at the left and right ends of the rotor are not identical, which, as previously explained, is a phenomenon due to the non-uniform inflow conditions of the wind tunnel.

5.2.2.
Active yaw control Next, we consider the wake velocity profiles for the 20 deg YM case.
As previously mentioned, the ε parameter was adjusted to better approximate the experimental results, while the other parameters remained identical to the baseline case. Figure 7 shows the numerically computed time-averaged streamwise velocity at three downstream locations; no PIV measurements are available in this case. The wake structure exhibits a characteristic 'kidney' shape, in accordance with the observations of other authors [4,15].   Figure 8 reports the first order (velocity profile) and second order (turbulence intensity TI and Reynolds shear stress u v ) flow quantities at 4D for the baseline case (top row of plots) and YM case (bottom row of plots). Experimental measurements were obtained with triple hot-wire probes [25], while simulations are reported with and without the presence of nacelle and tower. The baseline wake profile with nacelle and tower included correlates well with the measurement data. It appears that the presence of the nacelle and tower in the model is indispensable for an accurate solution (∆U = 2% in this case), as already noticed by other authors [8]. Although the situation is not as good for the YM case, matching with the experimental measurements is still quite reasonable (∆U = 4%). In addition, the prediction of the displacement of the wake center, computed as the location of the minimum of the Gauss fit to the wake, appears to be highly accurate and essentially identical to -0.26D for both the simulation and experiment.
Similar effects can be observed for the two second order quantities, TI and u v . Regarding the former quantity, the plots show that ambient TI, as seen at the end of the scan line, is nearly 0% for the simulation case and 2% for the experiments, as previously noted. Although there is more scatter for TI than for u v , the plots show once again the importance of including tower and nacelle in the model to improve accuracy.

Conclusion and outlook
This paper has presented the comparison of the results of a LES-lifting-line solver with wind tunnel measurements obtained with scaled wind turbine models. The analysis has considered a baseline flow-aligned case, the use of cyclic pitching to affect the wake, and the misalignment of the wind turbine with the incoming wind. Comparisons included integral rotor quantities as power and thrust, and local flow quantities measured by PIV and hot-wire probes. Measurements were performed in the near wake region, and after vortex breakdown in the far wake. Two scanning LiDARs were used to generate a map of the velocity at the inlet of the test section, to account for the non-uniformity of the wind tunnel generated flow. Based on the results shown in the paper, it appears that the CFD formulation is able to capture both the integral and local quantities with reasonable accuracy. In general, the baseline case exhibits the best results, while both CyPC and YM are qualitatively acceptable but not always very precise.
Overall, results appear to be very encouraging, although there are several effects that may potentially further improve the level of accuracy of the current numerical model. The LES framework can be improved on three fronts: using more accurate airfoil polars, which may be achieved with the existing polar-identification procedures, and introducing unsteady effects in the airfoil model. Results presented herein showed in a very clear manner the importance of including the nacelle and tower in the CFD model. However, the same results also show that the accuracy of such components is still not sufficient, and improvements may be obtained by the use of denser grids and possibly an improved formulation. Finally, the present simulation did not match the ambient turbulence of the experiment, which might be the source of additional discrepancies. For higher turbulence cases (which are closer to full scale operating conditions), eddies in the wind tunnel flow are generated by the use of spires at the inlet. For such conditions we have already developed a CFD model that can be used as a precursor of the lifting-line simulations. In the future, we will work on expanding that method to also account for the nonuniformity of the mean velocity, a problem that was considered here using a LiDAR-measured inflow map.