Comparative CFD study of the effect of the presence of downstream turbines on upstream ones using a rotational speed control system

The effect of a downstream turbine on the production of a turbine located upstream of the latter is studied in this work. This is done through the use of two CFD simulation codes, namely OpenFOAM and EllipSys3D, which solve the Navier-Stokes equations in their incompressible form using a finite volume approach. In both EllipSys3D and OpenFoam, the LES (Large Eddy Simulation) technique is used for modelling turbulence. The wind turbine rotors are modelled as actuator disks whose loading is determined through the use of tabulated airfoil data by applying the blade-element method. A generator torque controller is used in both simulation methods to ensure that the simulated turbines adapt, in terms of rotational velocity, to the inflow conditions they are submited to. Results from both simulation codes, although they differ slightly, show that the downstream turbine affects the upstream one when the spacing between the turbines is small. This is also suggested to be the case looking at measurements performed at the Lillgrund offshore wind farm, whose turbines are located unusually close to each other. However, for distances used in today's typical wind farms, this effect is shown by our calculations not to be significant.


Introduction
As larger wind turbine farms are being built, much effort is put into understanding the dynamics of the wake development downstream of turbines. However, little is known about how a turbine located downstream from another turbine affects the latter. The IEC norm 61400-12 [1] suggests that, for performance test purposes, a meteorological measurement mast should be located from two to four turbine diameters (D) upstream of a given turbine to avoid blocking effects from the turbine, and it recommends to use a distance of 2.5D. Medici et al. [2] recently investigated the blockage effect of a turbine in terms of wind velocity distribution upstream from the turbine using experimental wind tunnel data as well as actuator line computations. The experimental data suggested an important effect up to more than three D upstream from a given turbine, while the computations showed a smaller blockage effect. Differences were obtained between the experiments and computations, and concluded on the need to perform further work on this issue. To the author's knowledge, however, no work has been done on studying the effect from the blockage of a given turbine on the actual production of an upstream turbine. This issue is suggested to be studied here. This is done by performing CFD simulations using two different codes, i.e. EllipSys3D and OpenFoam, into which turbine rotors are modelled as actuator discs, and whose results will be compared. Measurements acquired at the Lillgrund offshore wind farm are also presented to shed light on this issue.
To perform such a study, it was important to ensure that the operational conditions of the simulated turbines are as realistic as possible. This regards, for example, the adequate representation of the rotational velocity of the discs, which is needed to determine the distribution of the angle of attack along the blades, which in turn is used to calculate their loading. This task, which is not straightforward for turbines operating in the wake of other turbines, is performed in this work using a generator torque controller that is implemented in the two codes below rated wind speed. To the author's knowledge, it is the first time that the proposed control strategy is implemented in a CFD method with turbine rotors modeled as actuator discs, although it has been implemented together with the actuator line method (see e.g. [3]), as well as in Blade Element Momentum based codes [4]. Therefore, some details are given about this strategy below, followed by comparative results concerning the effect of downstream turbines on upstream ones in terms of production.

Modelled turbine characteristics
The turbine under consideration in this study is the 2.3MW Siemens SWT93, installed in the Lillgrund offshore wind farm south of Sweden. It is a pitch-controlled variable speed wind turbine, whose blades are twisted and tapered. As the geometry of the blade of the SWT93 machine is unknown to the authors, geometric data is scaled down from the NREL 5MW model turbine presented in [4], and the corresponding airfoil data are used. The thrust and power curves of the actual and modelled turbines were compared and found to be reasonably close to each other.

Numerical models
Two numerical CFD codes are used in this work to model wind turbine aerodynamics. Making comparisons between results from two different codes is aimed at validating the assumptions on which each of them is based, in order to develop simulation tools that are as accurate as possible. The first one is EllipSys3D [5], a CFD code which solves the Navier-Stokes equations in their incompressible form using a finite volume approach. The second code is OpenFOAM [6], an open source CFD platform which also uses a finite volume approach to solve the Navier-Stokes equations. In both EllipSys3D and OpenFoam, the LES (Large Eddy Simulation) technique is used for modelling turbulence. There, the dynamically important eddies are explicitly solved whereas the impact of smaller eddies is modelled (through the use of a subgrid-scale (SGS) model). A different SGS model is used here in the simulations performed with EllipSys3D (Ta Phuoc [7]) and OpenFoam (Smagorinsky [8]). In Ellipsys, convective terms are solved using a mixture of third order QUICK (10%) and fourth order central difference scheme (90%), whereas in OpenFoam a second order QUICK scheme is used.
The same grid is used in both methods, see Fig. 1. It consists of a Cartesian grid of dimensions 17R × 17R × 30.4R respectively in the x(horizontal), y(vertical) and z(axial) directions, with a respective number of nodes of 64×64×256, with R the rotor radius. The grid is equidistant in the region of interest (shown in grey in Fig. 1) where the turbines are placed and where the wake calculations are performed. The resolution in the inner equidistant region is 0.1R, allowing 21 nodes on the diameter of the actuator disc. In order to save grid points, the grid is stretched towards all boundaries as seen in Fig. 1 where the xy (left figure) and yz (right figure) planes of the grid are shown, respectively. The equidistant region size is 3.1R × 3.1R × 24R. It is centered in the xy plane, as seen in Fig. 1(left), and it begins at a z value of 6.4R. The resolution of 0.1R was found to be satisfactory in a recently performed study, see Nilsson et al. [9]. In the latter, a grid sensitivity study showed that the use of a finer resolution and larger simulation domain had negligible impact on the production results. The total size and expansion of the domain relatively to the region of interest are also similar to the one used in this study.

Streamwise direction Vertical direction
Wind turbine 1 6

Vertical direction
3.1R

17R
Crosswise direction Figure 1. Sketch of the grid used in the simulations. Axes are non dimensionalized with respect to the rotor radius. The grey zone is the equidistant region, outside which the grid is stretched. In both methods, the wind turbine rotors are modelled as actuator disks. The actuator discs are centered in the xy plane, as seen in Fig. 1. As can also be observed in this figure, the upstream turbine is located 8R from the inlet, inside the equidistant region, and the downstream turbine is located different distances ∆S behind the first one. In the actuator disc method, tabulated airfoil data are used to determine the loading on the rotors by applying the blade-element method. In Ellipsys3D, forces are calculated on a local polar mesh representing the disc and later projected back on the cartesian mesh. Eleven points were used on the polar mesh radially on the actuator disc from the rotational axis to the edge of the disc, while 81 points were used in the azimuthal direction. The reader is referred to [10] for more information regarding this procedure. A Gaussian spreading of the forces was used in the axial direction to avoid the occurrence of wiggles in the velocity distribution around the disc, with a standard deviation equal to 0.2R/ √ 2 [10]. In OpenFoam, where the actuator disk model was recently implemented, the forces are calculated directly at the center of the cartesian cells contained in the actuator disc following the idea developed by Ammara [11]. No special treatment for the Cartesian cells located at the edge of the actuator disc is yet implemented, neither is a spreading of the forces axially around the actuator disc. In both numerical methods, a fixed value inflow was applied at the inlet. The boundary conditions were the following: farfield at the top and bottom, cyclic on the sides, and convective at the outlet. Neither the effect from shear nor the effect of atmospheric turbulence was modelled in this work.

Measurement data
Measurement data acquired at the Lillgrund offshore wind farm will be used to investigate the effect of downstream turbines on upstream ones. Characteristics associated with this farm's turbines are given in section 2. The inflow sectors considered here, i.e, 120±2.5 • and 222±2.5 • , are associated with turbines spacings of respectively 6.6R and 8.6R, which is very small compared to other typical offshore wind farms. For each of these wind directions, a line of 8 undisturbed turbines positioned normally to the incoming wind is available for study. An incoming velocity of 8±0.5m/s is considered, while the turbulence intensity was measured to be equal to approximately 5% for both inflow sectors. For more information about layout of the farm as well as the method used to filter the experimental data in the Lillgrund wind farm, the reader is referred to [9].

Strategy for the implementation of a controller
It is intended to model a system that makes it possible for a turbine whose aerodynamics is simulated in a CFD model to adapt its rotational velocity to the conditions in which it finds itself. We focus here on the operation of the turbine below rated power, so that the pitch of the modelled blades does not change. What is referred to as a controller in this work is therefore more exactly related to a strategy aimed at modelling the natural response of the system in order to simulate how the turbine would react in real operation. An actual controller per se will in a later step be implemented, that will rely on controlling the pitch of the blades so that the turbine operation can be simulated above rated wind speed.
The strategy suggested here relies on the calculated aerodynamic torque on the actuator disc for a given incoming wind and rotational velocity, and its comparison with the torque produced by the generator at this rotational velocity. No use is made of tabulated data that would for example provide the rotational velocity of a turbine as a function of the wind incoming to the wind turbine. The latter method would imply choosing a position where the undisturbed velocity (i.e., without induction from the turbine under consideration) for the given turbine would be determined, or performing simulations without the turbine in consideration in order to know the wind conditions prevailing at this position without the effect of the induction from this turbine. Both of the latter strategies can be difficult, especially under turbulent inflow conditions, and are avoided by the present method, which is also closer to modelling real-life operation of a variable speed turbine.
In the strategy for determining the rotational velocity of the turbines, let us suppose that at a given time step, the rotational velocity of a given turbine is Ω(t), and the wind velocity distribution at the rotor is V (t). The CFD codes will, with these data, calculate the aerodynamic torque T aero on the actuator disc. The torque produced by the generator (T gen ) is a function of the rotational velocity and can be determined from the generator torque curve (i.e., the generator torque as a function of Ω, see Fig. 2). If there is a difference between the generator torque and the aerodynamic torque, the system will experience a rotational acceleration (or deceleration), following Newton's second law for rotational bodies: T aero − T gen = (I rotor + I gen ) ∆Ω ∆t (1) where I rotor is the rotor inertia, I gen the generator inertia (a value of 10 6 kgm 2 is used in this work for the sum of these two moments of inertia), and Ω the rotational velocity in rad/s. After a time ∆t, Ω should reach a new value Ω(t + ∆t)=Ω(t)+ ∆Ω ∆t ∆t. At this new time step t + ∆t, Ω(t + ∆t) as well as the incoming wind velocity at this time, V (t + ∆t), will be used in the CFD codes to calculate a new value of the aerodynamic torque. The new value of the generator torque corresponding to Ω(t + ∆t) will be determined from the generator torque curve. If there is a difference between these two values, we will obtain a rotational acceleration following Eq.1. This rotational acceleration will be used to determine Ω value at the next time step, and so on as time goes forward. This principle that we will refer to as generator control is used below rated power. As mentioned above, in order to use this strategy, a generator torque curve is needed. As it was not available for the turbine that is simulated here, it was created. This was done following the method presented in Jonkman et al. [4] together with an in-house BEM code based on the book by Hansen [12]. The resulting curve is shown in Fig. 2. The optimal tip-speed ratio and corresponding power coefficient were determined with the BEM code for a zero pitch angle to be given respectively by 7.74 and 0.486. The cut-in and rated rotational velocities were respectively 6 and 16rpm, and rated wind velocity was 11m/s. An efficiency of 86% was used for  Figure 2. Generator torque curve created for the SWT93 wind turbine the generator and gearbox, linking the generator torque to a torque value on the upwind side of the rotor. Fig. 2 is composed of 5 control regions, as explained in Jonkman et al. [4]. Region 1 is before the cut-in wind speed, where the generator torque is null. Region 2 is the control region where power capture is optimized. The goal in this region is to operate the machine at the optimal tip speed ratio. The generator torque should then be proportional to the square of the generator speed, following the definition of the power coefficient. In region 3, which corresponds to rotational velocity values equal to and above rated value, power should be kept constant through the use of a pitch control system (which is outside the scope of the present article as mentioned earlier), making the generator torque inversely proportional to the generator speed (let us note that a gearbox ratio of 1 is used to facilitate calculations, so that the generator speed is the same as the rotor speed). Region 1 1/2 is a linear transition between regions 1 and 2. In this region, aerodynamic torque from the wind is used to accelerate the rotor for start-up. This region is said to be used in order to have a lower limit on the generator speed as to limit the operational speed range of the turbine. As in [4], it is extended from cut-in rotational velocity to 30% above this value. Region 2 then extends from the end of region 1 1/2 to the start of region 2 1/2. Region 2 1/2 is a transition region between regions 2 and 3. If operation at an optimum tip speed ratio were maintained too far to the right in region 2, the rotor speed would become too high and reach values above rated rotational speed for torque values below rated torque. The point at which region 3 starts is a known value, which corresponds to rated Ω and power. Therefore, a steep linear region has to be used from the end of region 2 to the beginning of region 3. This rotational velocity value at the beginning of region 2 1/2 is found using an analogy with an asynchronous generator, with a generator-slip percentage of 10%. Figs. 3 and 4 show respectively the rotational velocity and power as a function of time from the start of the simulation obtained for the first upstream turbine, as well as the turbine located downstream from the latter at different distances, for simulations performed with EllipSys3D and OpenFoam. The incoming flow velocity is constant at 10m/s, and the initial rotational velocity given to the two turbines is 9 RPM. It can be seen that a certain time is required for the rotational velocity of the two turbines to settle towards a value that then stays almost constant. The turbine located downstream from the first turbine takes more time to settle to such a value, as the flow passing through the first turbine has to reach the downstream turbine  Figure 3. RPM as a function of time for turbine 1 alone, and turbine 2 located at different distances from turbine 1. Calculations performed with EllipSys (top) and OpenFoam (bottom), V 0 =10m/s before the latter gets updated information about the condition in which it will be operating in the later stages of the simulation. As the distance of the downstream turbine increases, this time also does so. Moreover, because of the inertia of the rotating system, the response of the system not immediate. As expected, the downstream turbine, being submitted to a lower axial velocity, will rotate at a lower rotational velocity than the upstream one. It will in fact aim towards operating at optimal tip speed ratio, performing a coupling with the generator torque. The turbines are clearly seen to adapt to the environment in which they find themselves, which is the goal of the controller implemented here. It is shown to be capable of modeling the turbine  The aerodynamic power calculated at the two turbines is shown in Fig. 4. It is determined on the actuator discs by multiplying the aerodynamic torque produced by the forces acting tangentially on the disc cells by the rotational velocity of the discs. Let us recall that the latter depends on the coupling performed in the control mechanism between the aerodynamic and generator torques. The power is seen to follow the same trend as the rotational velocity. One particular thing to mention is that it is decreasing for the downstream turbine as the distance from the upstream turbine increases, although the difference is rather small between the turbines that are located from 5R downstream from the first turbine. These results suggest that the flow velocity decreases over a certain distance downstream from the first turbine, after which it reaches a plateau. A great distance behind the turbine is then expected to be needed for the flow to recover to its original state. Wake losses are not seen to decrease within the distances travelled by the flow which are tested here. Let us remind that the simulations performed here do not include neither shear nor atmospheric turbulence. It is expected that introducing these features would contribute to breaking up the vortex structures in the wake and decrease the wake effects that the downstream turbine are subjected to, as well as decrease the effect from the latter on the upstream turbine. Differences are obtained between the rotational speed and power reached by the turbines in EllipSys and OpenFoam, although both codes predict rotational velocities and power values that are coherent with each other, comparing with Fig. 2. The differences are being investigated at this point, and might be due among other things to the absence of smearing for the actuator disc in the OpenFoam calculations. The latter is also expected to be the source of the small oscillations observed in the power predicted by OpenFoam, related to wiggles obtained near the actuator disc in the axial velocity profile.

Results and discussion
On the other hand, power and RPM values obtained from the EllipSys3D calculations for the downstream turbines are not seen to oscillate around their mean value, as might be expected for a turbine located in the wake of another turbine. This is expected to follow from the fact that the aerodynamic torque used in Eq.1 is calculated at each time step over the entire actuator disc, most probably resulting in a smoothing effect. Indeed, the wind speed might vary spatially and timewise over the whole rotor, but as all cells of the actuator disc are used in the calculation of the total aerodynamic torque, the latter might not change in an important way from one timestep to the next. It will also be interesting in a future work to introduce a turbulent inflow at the inlet of the domain with increasing turbulence intensities and study the corresponding response of the system, as well as studying the effect of the implementation of a smearing of the forces around the actuator disc in OpenFoam, which is being performed. Fig. 5 presents the aerodynamic power produced at the first turbine when a second turbine is located at different distances downstream from the former. The power produced by the first turbine seems to be influenced by the turbines located closely downstream from this turbine. The EllipSys3D calculations predict that a turbine located 2R or 3R downstream is seen to make the power produced by the first turbine to decrease by about 5.9 and 3.4% respectively. The OpenFoam calculations predict in this case a decrease of about 4.6 and 2.5%. The smaller effect from the downstream turbine as calculated from OpenFoam is expected to be due to the fact that no axial spreading of the forces is used for the moment for its actuator disc implementation, as opposed to EllipSys3D. Precisely, by applying a convolution to the forces of the actuator disc in EllipSys3D, the region where the forces act is extended. In consequence, an increased influence of the downstream turbine on the upstream rotor is expected. The small difference observed is an indicator that using such a spreading, which solves numerical problems associated with wiggles in the computation of the velocity field, would consist in a reasonable compromise. The effect from the turbines located further downstream seems significantly smaller, which suggests that downstream turbines do not have a significant effect on upstream turbines considering the spacings considered in today's typical wind farms.
Measurement data acquired at the Lillgrund offshore wind farm can allow us to see if experimental data show the same trend as predicted by our simulation methods. Table 1 shows the power averaged over the undisturbed (upstream) first row turbines in the two main wind directions studied here, for which downstream turbines are located at two different distances, as mentioned above. Unfortunately, production values acquired in the absence of downstream turbines are not available, so that only production values in the presence of downstream turbines can be looked at. We can see a significant difference in the mean power values associated to each  Figure 5. Power as a function of time for turbine 1 when turbine 2 is located different distances downstream from turbine 1. Calculations performed with EllipSys (top) and OpenFoam (bottom), V 0 =10m/s spacing, where a smaller spacing is associated to a lower production. One should however be careful with these results to which some experimental uncertainties are associated, as discussed in [9]. They were also acquired at a different mean velocity as the one used in the simulations, and in the presence of atmospheric turbulence. Nonetheless, the same trend as observed in the simulations is seen in these experimental results, i.e., a decrease in turbine spacing results in a decrease in the production of the upstream turbine. While the spacing used in Lillgrund is smaller than the one used in typical offshore wind farms, it is shown here to be helpful in investigating the effect of downstream turbines on upstream ones.  Table 1. Power production averaged over the eight front row (upstream) turbines in Lillgrund [13] for an incoming mean wind velocity of 8±0.5m/s and turbulence intensity of approximately 5% for the two main wind directions associated to two turbine spacings. The uncertainty on the production of each front row turbine is first determined from calculating the standard error of the mean. The uncertainty on the mean production averaged over the eight front row turbines is then estimated for the two extreme cases where the individual uncertainties are either fully correlated (uncertainty 1, worst case scenario) or uncorrelated (uncertainty 2, best case scenario).

Conclusion
The effect of the presence of a turbine located downstream from another turbine on the production of the latter was studied using two different CFD simulation codes, i.e., EllipSys3D and OpenFoam, into which the turbines were modelled as actuator discs. The rotational velocities of the discs were determined throughout the simulations using a newly implemented generator torque controller based on comparing the generator and aerodynamic torques below rated wind speed, which was seen to perform as expected. Differences were obtained in the results produced by the two simulation codes, and more work will be performed in the future in comparing and validating both codes. The results obtained from the two codes however lead to the same conclusions. They suggest that the production of an upstream turbine is affected by the presence of a closely located downstream turbine. A decrease in power of up to 5.9% was observed from the EllipSys3D calculations for a turbine located 2R downstream. However, this effect decreased rapidly with increasing turbine spacing, and was found not to be critical for typical spacings used in today's wind farms (typically over 14 R for offshore wind farms). It is also expected to decrease in the presence of atmospheric turbulence (not modelled here) that would, by contributing to break down the wake behind the downstream turbine, decrease its blockage effect. Experimental results from the offshore Lillgrund wind farm, whose turbines are located untypically close to each other, also suggested a decrease of the influence of downtream turbines on upstream ones with increasing turbine spacing.