Wake Statistics of Different-Scale Wind Turbines under Turbulent Boundary Layer Inflow

Subscale wind turbines can be installed in the field for the development of wind technologies, for which the blade aerodynamics can be designed in a way similar to that of a full-scale wind turbine. However, it is not clear whether the wake of a subscale turbine, which is located closer to the ground and faces different incoming turbulence, is also similar to that of a full-scale wind turbine. In this work we investigate the wakes from a full-scale wind turbine of rotor diameter 80 m and a subscale wind turbine of rotor diameter of 27 m using large-eddy simulation with the turbine blades and nacelle modeled using actuator surface models. The blade aerodynamics of the two turbines are the same. In the simulations, the two turbines also face the same turbulent boundary inflows. The computed results show differences between the two turbines for both velocity deficits and turbine-added turbulence kinetic energy. Such differences are further analyzed by examining the mean kinetic energy equation.


Introduction
As wind energy grows as a main energy resource for the whole world, further research is still needed to reduce the cost of wind energy to keep it economically competitive [1]. Turbine wakes affect both the power production and operation and maintenance costs of wind energy. Wind tunnel experiments at meter-scale or even smaller wind turbines [2][3][4][5][6][7][8] and field measurements of subscale wind turbines (e.g., the SWiFT facility [9][10][11][12]) play a vital rule in understanding the dynamics of turbine wakes and provide valuable datasets for validating computational models. However, the sizes of these meter-scale turbines and subscale turbines are often much smaller than a full-scale wind turbine, which calls into question how well these small-scale wind turbines can represent full-scale wind turbines [13].
When designing a meter-scale or a subscale wind turbine, geometric, kinematic and dynamic similarities should be maintained to ensure their equivalence to a full-scale wind turbine. For a meterscale wind turbine, which can be about 1000 times smaller than a full-scale wind turbine, it is even where x i (i = 1, 2, 3) are the Cartesian coordinates, ξ j (j = 1, 2, 3) are the curvilinear coordinates, U i denotes the contravariant volume flux, u i are the velocity vector in Cartesian coordinates, J is the Jacobian of the geometric transformation, ξ i l represent the transformation metrics, g jk represents the contravariant metric tensor, ρ is the density, µ is the dynamic viscosity, p is the pressure, f l (l = 1,2,3) are the body forces from actuator type turbine models, and τ ij is the subgrid scale stress modeled using the dynamic subgrid scale model [26]. The governing equations are discretized in space using a second-order central differencing scheme and in time using a second-order accurate fractional step method. The pressure Poisson equation is solved using an algebraic multigrid acceleration along with GMRES solver. The momentum equation is solved using the matrix-free Newton-Krylov method. The VFS-Wind code has been validated extensively using laboratory and field measurements [21,24] and applied to utility-scale wind turbines [22,27,28].
The actuator surface models for turbine blades and nacelle developed in [25] are employed in this work. In the actuator surface model for turbine blades, the blade is represented with the actuator surface defined by chords at different radial locations. Forces computed using the blade element method are distributed on the actuator surface to represent the effect of blades on the incoming flow. Compared with the actuator line parameterization, the actuator surface can represent better the geometrical effect of the blade in the chordwise direction. A model for the nacelle is also necessary for accurately predicting turbine wakes, which affects the hub vortex and meandering in the far wake [29]. In this work the nacelle is represented using an actuator surface formed by the actual surface of the nacelle. The effects of nacelle on the incoming flow are modeled using distributed forces with the wall-normal component calculated by satisfying the no-flux boundary condition and the wall-tangent component calculated by specifying a friction coefficient.

Specifics of the Employed Turbine Design
A full-scale wind turbine and a subscale wind turbine are studied in this work as shown in Figure 1. The diameters of the two turbines are 80 m and 27 m, respectively. The full-scale turbine T80 is located 93.3 m above the ground, while the subscale turbine T27 is located 31.5 m above the ground. The T27 turbine is designed by Kelley et al. in [30] (Design A). The T80 turbine is designed by simply scaling the T27 turbine, such that the radial distributions of the twist angle and the normalized chord length (as shown in Figure 2) are the same for the two turbines. It is noticed that the T27 turbine has the same rotor diameter and hub height as the SWiFT turbine but with different designs as shown in Figure 2 and Table 1. The scaling ratio is chosen in a way the size of the full-scale turbine is similar to that of a land-based wind turbine. The radial distribution of the airfoil shapes is shown in Table 1. For the present cases, the tip speed ratio is 9. The power coefficient, C P = P 0.5ρπR 2 U 3 (where P is the power, R is the rotor radius, D is the rotor diameter, and U is the incoming wind speed), is 0.47 and 0.49 for the turbine T80 and T27, respectively. The thrust coefficient, C T = T 0.5ρπR 2 U 2 (where T is the thrust), is 0.65 and 0.68 for the turbine T80 and T27, respectively. The incoming wind speed U employed for calculating C P and C T is obtained by averaging the time-averaged velocity over a disk of radius R parallel with the rotor rotating plane and located 1D upwind of the turbine. The same lift and drag coefficients are employed for both turbines to ensure dynamic similarity and to simplify the problem such that the only two differences are the turbine size and the inflow.    [30]. The information in parentheses shows the airfoil shape distribution of the SWiFT turbine [31].

Case Setup
The sizes of the computational domain are L x × L y × L z = 1.20δ × 2.09δ × δ and L x × L y × L z = 0.405δ × 2.09δ × δ with the corresponding numbers of grid nodes N x × N y × N z = 1501 × 361 × 282 and N x × N y × N z = 1501 × 351 × 256 for T80 and T27 turbine simulations, respectively, where x, y and z represent the downwind, lateral and vertical directions, respectively, and δ = 1000 m is the boundary layer thickness. The width (L y ) and the height (L z ) of the computational domain are set being the same as those of the precursory simulation. In the downwind direction, the length of the computational domain (L x ) is 15 rotor diameters for both cases. The turbine is located 2D from the inlet plane. The length of the computational domain, although it is much shorter than its width, is typical for simulations of stand-alone wind turbines [4,22], and will not affect incoming large-scale coherent structures as they are generated from a precursory simulation. In the downwind direction, the grid nodes are uniformly distributed with grid spacing D/100. In the other two directions, the grids are uniformly distributed in the near turbine region (|y − y t | < 0.75D and z < 2.2D, where y t is the turbine coordinate in the lateral direction) with the grid spacing ∆h = D/100 and gradually stretched to boundaries of the computational domain. Employing the actuator surface model requires higher spatial resolutions. A grid of spacing D/160 was employed in [21], where the same actuator surface model is applied to predict the coherent tip vortices of a utility-scale wind turbine. In [25], on the other hand, a grid of spacing D/40 was shown being enough for the same actuator surface model to accurately predict the wake statistics of a hydrokinetic turbine, e.g., the velocity deficit and turbulence kinetic energy, which are also of interest in this work. Based on this study and the other work for a utility-scale wind turbine [27], we believe that the resolution of the employed grid is enough for the quantities of interest in this work. The sizes of time step are ∆tU h /D = 7.9 × 10 −4 and 5.9 × 10 −4 for T80 and T27 cases, respectively, where U h and D are their own incoming downwind velocity at hub height and diameter, respectively.
The roughness length of the ground is k o = 0.01 m, which represents the field characteristics at the SWiFT site. The thermal stratification is neutral. The free-slip boundary condition and the logarithmic law for rough walls are applied at the top and bottom of the computational domain, respectively. The periodic boundary condition is applied in the lateral direction. The same inflow is applied in the two simulations, which is generated from a precursory simulation. In the precursory simulation, the computational domain is L x × L y × L z = 6.28δ × 2.09δ × δ. The grid nodes are uniformly distributed in the horizontal directions, while are stretched in the vertical direction with the height of the first off wall grid cell ∆z 1 = δ/1000. Periodic boundary conditions are applied in the horizontal directions. At the ground and top boundary of the domain, boundary conditions the same as those in the wind turbine simulations are employed. The size of time step is ∆tU h /δ = 8 × 10 −4 . The statistics of the inflow generated from the precursory simulation is shown in Figure 3. In the precursory simulation, the velocity fields on a y − z plane at every time step are saved. Since the spatial and temporal resolutions of the precursory simulation are different from those of the turbine simulations, linear interpolations are employed in both space and time in order to apply the obtained inflow on the inlet plane of the wind turbine simulation.
In both cases, simulations are first carried out to a full developed state, and then further advanced in time for 100 and 70 rotor revolutions to compute time-averaged quantities for T80 and T27 cases, respectively. It is noted in Section 4 that the profiles of turbulence kinetic energy at 1D rotor upwind are not smooth enough. Although a similar number of rotor revolutions for temporal averaging has been employed in a hydrokinetic turbine case [32], further averaging in time is probably needed especially to take into account the low frequency variations caused by incoming large eddies. However, it will require about 1000 rotor revolutions [27] or even more, which is time-consuming and not realistic for this study. Taking the T80 case as an example, as one rotor revolution needs about 3.5 h wall clock time by using 640 compute cores, simulating 1000 rotor revolutions will take about 5 months wall clock time. As shown in Section 4, the profiles of turbulence kinetic energy and other quantities in the turbine wake are fairly smooth, so we believe that the main conclusions drawn from this work are not affected by the length of temporal averaging.

Results
As no measurements are available for the employed turbine designs, we attempt to validate the employed computational setup by comparing the simulated velocity deficit ∆U with the measurements at the SWiFT site considering the T27 turbine and the SWiFT turbine are of the same size and have comparable blade designs, in which the velocity deficit ∆U is defined as where U is the time-averaged downwind velocity at different downwind locations, U in is the time-averaged incoming downwind velocity (which is taken at 1D upwind the turbine). As in [33], a slight offset of 0.1D is imposed in the negative y direction to compensate for the wake deflection observed in the measured data. It is seen in Figure 4 that the simulation velocity deficit profiles show an overall good agreement with the measurements considering the complex wind and turbine operating conditions in the field and different turbine designs. It is also noticed the lateral velocity deficit profiles from the T80 case and the T27 case are very similar at the considered turbine downwind locations. The measured data are digitized from [33]. Details on the measurements can be found in [34]. Red lines: T80; Blue lines: T27; Circles: measurements.
We then examine the contours of the instantaneous downwind velocity from the two cases in Figure 5. It is seen that the wake remains an annular shape until about 2D downwind of the T80 turbine, while becomes unstable immediately downwind of the T27 turbine. At further downwind locations, the flow structures of the T80 turbine's wake remain quite coherent with its center biased towards the ground, which, on the other hand, looks chaotic and meanders significantly in the vertical direction with its center shifting above the centerline of the rotor for the T27 turbine. Next, we examine the time-averaged quantities from the two cases. In Figure 6 we compare the time-averaged downwind velocity from the two cases. It is seen that the velocity deficits from the T27 case are higher for the upper parts of the wake when x/D < 2, which are similar for the T80 case. At far wake locations, the wake center from the T80 case is below the hub height, while above the hub height for the T27 case. In Figure 7 we examine the turbulence kinetic energy (TKE) from the two cases. One similar observation from the two cases is that the high TKE region are in the top shear layer for both cases. One significant difference between the two cases is that the high TKE region persists to about 7D downwind the turbine for the T80 case, which only persists to about 4D turbine downwind for the T27 case. To qualitatively show the differences between the two cases, we examine the time-averaged velocity deficit (Equation (3)) and the turbine-added TKE profiles. The turbine-added TKE (∆k) is defined as ∆k = k − k in (4) where k is the TKE at different downwind locations, k i n is the TKE of the inflow (which is taken at 1D upwind of the turbine for the present cases). It is seen in Figure 8a that incoming time-averaged downwind velocity profiles normalized using the corresponding length and velocity scales are quite similar, with one difference that the incoming velocity is slightly lower for the lower part of the turbine for the T27 case. At x/D = 2, the velocity deficit profiles are similar between the two cases with slightly higher velocity deficits from the T27 case. Moving to further downwind locations, the center of the wake gradually shifts to the upper part of the wake, which is above the upper tip of the turbine at x/D = 10. In Figure 9 we compare the turbine-added TKE from the two cases. First, we examine the normalized TKE at the 1D upwind of the turbine in Figure 9a. Because of different turbine sizes and hub heights, it is seen that the TKE of the inflow in the turbine region are quite different between the two cases, that the inflow TKE for the T27 case is higher than the T80 case for locations above the hub height while slightly lower for locations below the hub height. At x/D = 2, the vertical distributions from the two cases are similar, in that two peaks exist within the turbine top and bottom shear layers, respectively, with the one on the top about 1.5 times higher. At this location, it is also noticed that the normalized ∆k from the T27 case is higher at almost all z locations. Starting from x/D = 4, one observation is that the bottom peak of ∆k disappears for both cases. One major difference is that the top peak of ∆k still exists for the T80 case, which becomes flat and with a much wider high ∆k region for the T27 case. It is also noticed that the maximum ∆k from the T80 case is much higher than that from the T27 case. These observations are consistent with the lower incoming TKE for the upper region of the T80 case allowing the wake to maintain coherent helical wake vorticity structures for longer distance downstream than for the same region of the T27 case, where the higher incoming turbulence likely drives these structures to mix in to more evenly distributed vorticity sooner. To explore the reason for different distributions of velocity deficits and turbine-added TKE, we examine the mean kinetic energy (MKE) equation integrated over y − z plane, which is shown as follows: where MC, PT, TC, DF, TP, DP are the convection of the MKE by the mean flow, transport terms due to mean pressure, turbulence fluctuations, the diffusion term, the negative of turbulence production term (which transfers energy from the mean flow to turbulence), and the dissipation term, respectively. The expressions for the various terms in Equation (5) are given as follows: where y t is the coordinate of the rotor center in the lateral direction, R is the rotor radius, S ij = ∂ u i /∂x j + ∂ u j /∂x i /2 is the strain rate tensor. The turbulence convection term TC is further decomposed into three components for the contributions from three directions as follows: where First, we show in Figure 10 different terms in Equation (5) for both cases. It is seen that the mean convection (MC) term is balanced with the pressure transport term in the near wake region where the pressure recovers to the ambient pressure by extracting mean kinetic energy from the wake. In the far wake region, the MC term is mainly balanced with the turbulence convection (TC) term. The above observations are similar to the wake of a model wind turbine located downwind of a three-dimensional hill [20]. The major differences between the T80 and the T27 cases are observed in the TC and TP terms, of which the magnitudes are higher in the near wake for the T27 case.
In Figure 11 we compare the three components of the TC term between the T80 and T27 cases. It is seen that the mean kinetic energy losses are due to the downwind component of the turbulence convection term TC x term when x/D < 4 and x/D < 2 for the T80 and T27 cases, respectively. At further downwind locations the effect of the TC x term on the MKE budget is negligible. On the other hand, the other two components of the turbulence convection terms TC y and TC z contributes positively to the MKE budget at almost all downwind locations except in the region immediately downwind the turbine. The TC y term from the T27 case is higher than that from the T80 case for all downwind locations. On the other hand, the TC z term from the T27 case is almost the same as that from the T80 case for x/D < 2, while lower than that from the T80 case for further downwind locations. That MKE entrainment from the top is lower than that from two sides explains the upward shift of the wake center observed in Figure 8 for the T27 case.  (6)); Line with circles: PT term (Equation (7)); Line with stars: TC term (Equation (8)); Line with squares: DF term (Equation (9)); Line with diamonds: TP term (Equation (10)); Line with triangles: DP term (Equation (11)). Different terms are normalized using 1 2 DU 3 h . Figure 11. The three components of the turbulence convection (TC) term for (a) TC x (Equation (13)), (b) TC y (Equation (14)) and (c) TC z (Equation (15) To better understand the differences in TC y and TC z terms between the two cases, we examine the averaged Reynolds stress and downwind velocity on the control surface from the two cases in Figures 12 and 13, in which different terms are defined as follows: As seen in Figure 12a, the averaged downwind velocity on the two sides of the surface are nearly the same for the T80 and T27 cases. The averaged Reynolds stress term u v side from the T27 case, on the other hand, is larger than that from the T80 case as shown in Figure 13a, which is the key reason the TC y term the T27 case is larger than that from the T80 case (shown in Figure 11b). On the bottom of the control surface, the averaged downwind velocity from the T27 case is larger than that from the T80 case (as shown in Figure 12b), while the u w bot from the T27 case is smaller than that from the T80 case from 2D to 5D rotor downwind and close to zero at further downwind locations, which is the same as the T80 case (shown in Figure 13b). This makes that the differences in entrainment from the bottom surface are insignificant between the two cases especially at far wake locations. On the other hand, the averaged downwind velocity on the top surface is significantly lower for the T27 case at all downwind locations as shown in Figure 12c, which is caused by the upward shift of the wake center for the T27 case as shown in Figure 8. Meanwhile, the u w top from the T27 case is also significantly lower than that from the T80 case at downwind locations x/D > 4. Therefore, this explains why the TC z term from the T27 case is lower than that from the T80 case at downwind locations x/D > 2 (shown in Figure 11c).
To explore the reason for the different downwind variations of the TKE shown in Figure 9, we compare the turbulence production term from the two cases in Figure 14. It is observed that the magnitude of the TP term from the T27 case is significantly larger than that from the T80 case when x/D < 2, while is similar to that from the T80 case at further downwind locations. Higher TP term indicates more energy is transferred to TKE from MKE. This is consistent to what we observed in Figure 9 that the TKE from the T27 case is significantly higher than that from the T80 case in the near wake region.    (10)). Red line: T80; Blue line: T27. The TP term is normalized using 1 2 DU 3 h .

Conclusions
In this work we investigated the wake from a full-scale turbine (T80) and a subscale turbine (T27) using the VFS-Wind code to carry out large-eddy simulation with actuator surface models for turbine blades and nacelle. The T80 turbine is 2.96 times larger than the T27 turbine in terms of both rotor diameter and hub height. In both cases, the same turbulent inflow is employed. The same lift and drag coefficients are also employed for the two turbines. The key differences between the inflows for the two cases are the different distributions of wind shear and turbulence intensity across the rotor caused by different rotor diameters and different hub heights. The computed results show differences between the two cases for both velocity deficits and turbine-added turbulence kinetic energy. It is observed that the wake center of the T27 turbine shifts upwards in the far wake of the turbine (e.g., more than 0.5D above the hub height at 10D downwind the rotor), which remains at hub height for the T80 turbine. The maximum turbulence kinetic energy in the wake of the T27 turbine is higher than that from the T80 case in the near wake region (e.g., more than 20% higher at 2D downwind the rotor), but decreases rapidly and becomes lower than that from the T80 turbine at further downwind locations (e.g., about 50% lower at 8D downwind the rotor). We explore the reason for these differences by examining the budgets for the mean kinetic energy. Compared with the T80 case, we found that the MKE entrainment for the T27 case is higher for the lateral component of the turbulence convection term, while lower for the vertical component of the turbulence convection term. This explains the upward shift of wake center for the T27 case. We also examined the turbulence production term. It is observed that the magnitude of the turbulence production term from the T27 case is significantly higher than that from the T80 case in the near wake region, which explains the higher turbulence kinetic energy in the near wake for the T27 case. Both the inflow turbulence and the size of the turbine relative to the incoming eddies may cause these observed differences in the wake characteristics. Further study (e.g., space-time correlation study [35]) is needed to probe in more depth the underlying cause for the differences between the two turbines of different scales. A systematic study on the wakes from turbines of different scales for different turbulent inflows is also needed to further investigate the phenomenon observed in this work.

Conflicts of Interest:
The authors declare no conflict of interest.