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

Abstract

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.

1. 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 meter-scale wind turbine, which can be about 1000 times smaller than a full-scale wind turbine, it is even challenging to ensure only geometric similarity (e.g., the rotor diameter and chord length) and the kinematic similarity (e.g., the tip speed ratio $\lambda =\mathsf{\Omega }R/U$, where Ω is the rotor rotating speed, R is the rotor radius and U is the incoming wind speed). As such, researchers often build a meter-scale wind turbine in a way that the power coefficient and thrust coefficient are similar to that of a full-scale wind turbine [6]. For the scale effects on turbine wakes, Howard and Guala [14] compared the velocity deficits of a meter-scale turbine with that of the University of Minnesota 2.5 MW EOLOS turbine and observed significant differences at $x/D=1.5$ and $x/D=2.5$ turbine downwind locations but relatively small differences at $x/D=3$ turbine downwind location, where D is the rotor diameter. Heisel, Hong and Guala [15] further compared the wake meandering from the meter-scale turbine and the EOLOS turbine and observed similar meandering frequencies related to the bluff body shear layer instability.

For a subscale wind turbine, which may be 3 or 4 times smaller than a full-scale wind turbine, the geometric similarity and the kinematic similarity can be ensured relatively easier compared to meter-scale wind turbines. However, because of the difference in Reynolds number, the dynamic similarity, namely the distributions of lift and drag coefficients along the blade for wind turbines, still cannot be guaranteed easily while keeping the geometric and kinematic similarities. Kelley et al. [16] proposed to loosely maintain the dynamic similarity by keeping the same dimensionless bound circulation along the blade by relaxing the constraints on the geometric similarity. Although the aerodynamics of a subscale wind turbine can be made closer to a full-scale wind turbine, an important difference is that the subscale wind turbine is located closer to the ground where the mean shear stress and the turbulence intensity change significantly with distance from the ground. However, how this difference affects the representation of a subscale wind turbine wake to that of a full-scale wind turbine is not clear yet. To address this issue, in this work we simulate a subscale wind turbine and a full-scale wind turbine, which are geometrically and kinematically equivalent, and are dynamically equivalent by applying the same lift and drag coefficients in turbine parameterizations, and under the same turbulent boundary layer inflow.

It is noted that this work is different from the studies in the literature on investigating the effects of inflow turbulence on turbine wakes. In the literature, the turbine wakes under different turbulent inflows, which are due to different ground roughness lengths, were investigated in [2,17] using wind tunnel experiments and large-eddy simulation, respectively. The effects of different inflow turbulence on turbine wakes, which is caused by different thermal stratifications, were studied using wind tunnel experiments and large-eddy simulation in [18,19], respectively. In [20], turbine wakes under different inflow turbulence caused by an upwind hill of different heights were simulated using large-eddy simulation. In [21], the coherent tip vortices of a utility-scale wind turbine were investigated for inflows of different turbulence intensities. In [22], the wake meandering of a utility-scale wind turbine was investigated for inflows of different turbulence intensities. All these studies were focused on the wake of one-size wind turbine. To the best of our knowledge, no studies have been carried out on comparing wakes from turbines of different sizes, which is critical as we use the knowledge from a subscale wind turbine to a full-scale wind turbine.

The rest of the paper is organized as follows: in Section 2 we briefly describe large-eddy simulation with actuator surface models for turbine blades and nacelle; we then compare the results from the two turbines of different sizes in Section 4; at last we draw conclusions from this work in Section 5.

2. Numerical Method

The Virtual Flow Simulator (VFS-Wind) code [23,24,25] is employed in this work to simulate wakes from the two turbines, which solves the incompressible Navier–Stokes equations in curvilinear coordinates shown as follows:

$$J\frac{\partial U^j}{\partial {\xi }^j}=0,$$

(1)

$$\frac{1}{J}\frac{\partial U^i}{\partial t}=\frac{{\xi }_l^i}{J}\left(-\frac{\partial (U^j{u}_l)}{\partial {\xi }^j}+\frac{\mu }{\rho }\frac{\partial }{\partial {\xi }^j}\left(\frac{g^{jk}}{J}\frac{\partial u_l}{\partial {\xi }^k}\right)\right.\left.-\frac{1}{\rho }\frac{\partial }{\partial {\xi }^j}\left(\frac{{\xi }_l^{j}p}{J}\right)-\frac{1}{\rho }\frac{\partial {\tau }_{lj}}{\partial {\xi }^j}+\frac{1}{\rho }{f}_l\right),$$

(2)

where $x_i$ ($i=1,2,3$) are the Cartesian coordinates, ${\xi }^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, ${\xi }_l^i$ 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 ${\tau }_{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.

3. Test Cases and Computational Details

3.1. 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. Distribution of airfoil shapes along the radial direction [30]. The information in parentheses shows the airfoil shape distribution of the SWiFT turbine [31].

Section	r/R	Shape
1	$0\le r/R<0.125(0\le r/R<0.038)$	Circle (Circle)
2	$0.125\le r/R<0.225(0.038\le r/R<0.231)$	DU 40 (NACA 63-235)
3	$0.225\le /R<0.375(0.231\le r/R<0.615)$	DU 35 (NACA 63-224)
4	$0.375\le /R<0.475(0.615\le r/R<0.923)$	DU 30 (NACA 63-218)
5	$0.475\le /R<0.675(0.923\le r/R<1)$	DU 25 (NACA 63-214)
6	$0.675\le /R<0.825$	NACA ${63}_6-621$
7	$0.825\le /R\le 1$	NACA ${63}_6-618$

. 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=\frac{P}{0.5\rho \pi 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=\frac{T}{0.5\rho \pi 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.

3.2. Case Setup

The sizes of the computational domain are $L_x\times L_y\times L_z=1.20\delta \times 2.09\delta \times \delta$ and $L_x\times L_y\times L_z=0.405\delta \times 2.09\delta \times \delta$ with the corresponding numbers of grid nodes $N_x\times N_y\times N_z=1501\times 361\times 282$ and $N_x\times N_y\times N_z=1501\times 351\times 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 $\Delta 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 $\Delta t{U}_h/D=7.9\times {10}^{-4}$ and $5.9\times {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\times L_y\times L_z=6.28\delta \times 2.09\delta \times \delta$. 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 $\Delta z_1=\delta /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 $\Delta t{U}_h/\delta =8\times {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.

4. 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 $\Delta 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 $\Delta U$ is defined as

$$\Delta U=U_{in}-U$$

(3)

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.

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 ($\Delta k$) is defined as

$$\Delta 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 $\Delta 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 $\Delta k$ disappears for both cases. One major difference is that the top peak of $\Delta k$ still exists for the T80 case, which becomes flat and with a much wider high $\Delta k$ region for the T27 case. It is also noticed that the maximum $\Delta 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:

$$0=MC+PT+TC+DF+TP+DP,$$

(5)

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:

$$MC=-{\int }_{y_t-R}^{y_t+R}{\int }_{z_h-R}^{z_h+R}\langle u_j\rangle \frac{\partial \left(\langle u_i\rangle \langle u_i\rangle /2\right)}{\partial x_j}dzdy$$

(6)

$$PT=-{\int }_{y_t-R}^{y_t+R}{\int }_{z_h-R}^{z_h+R}\frac{\partial \left(\langle p\rangle \langle u_j\rangle /\rho \right)}{\partial x_j}dzdy$$

(7)

$$TC=-{\int }_{y_t-R}^{y_t+R}{\int }_{z_h-R}^{z_h+R}\frac{\partial \left(\langle u_i^{\prime }{u}_j^{\prime }\rangle \langle u_i\rangle \right)}{\partial x_j}dzdy$$

(8)

$$DF=2{\int }_{y_t-R}^{y_t+R}{\int }_{z_h-R}^{z_h+R}\frac{\partial \left(\left(\nu +{\nu }_t\right)S_{ij}\langle u_i\rangle \right)}{\partial x_j}dzdy$$

(9)

$$TP=2{\int }_{y_t-R}^{y_t+R}{\int }_{z_h-R}^{z_h+R}\langle u_i^{\prime }{u}_j^{\prime }\rangle \frac{\partial \langle u_i\rangle }{\partial x_j}dzdy$$

(10)

$$DP=-2{\int }_{y_t-R}^{y_t+R}{\int }_{z_h-R}^{z_h+R}\left(\nu +{\nu }_t\right)S_{ij}\frac{\partial \langle u_i\rangle }{\partial x_j}dzdy$$

(11)

where $y_t$ is the coordinate of the rotor center in the lateral direction, R is the rotor radius, $S_{ij}=\left(\partial \langle u_i\rangle /\partial x_j+\partial \langle u_j\rangle /\partial x_i\right)/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:

$$TC=T{C}_x+T{C}_y+T{C}_z,$$

(12)

where

$$T{C}_x=-{\int }_{y_t-R}^{y_t+R}{\int }_{z_h-R}^{z_h+R}\frac{\partial \left(\langle u_i^{\prime }{u}^{\prime }\rangle \langle u_i\rangle \right)}{\partial x}dzdy,$$

(13)

$$T{C}_y=\left({\int }_{z_h-R}^{z_h+R}\left(-\langle u_i^{\prime }{v}^{\prime }\rangle \langle u_i\rangle \right)dz\right){|}_{y_t-R}^{y_t+R},$$

(14)

$$T{C}_z=\left({\int }_{y_t-R}^{y_t+R}\left(-\langle u_i^{\prime }{w}^{\prime }\rangle \langle u_i\rangle \right)dy\right){|}_{z_h-R}^{z_h+R}.$$

(15)

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 $T{C}_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 $T{C}_x$ term on the MKE budget is negligible. On the other hand, the other two components of the turbulence convection terms $T{C}_y$ and $T{C}_z$ contributes positively to the MKE budget at almost all downwind locations except in the region immediately downwind the turbine. The $T{C}_y$ term from the T27 case is higher than that from the T80 case for all downwind locations. On the other hand, the $T{C}_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.

To better understand the differences in $T{C}_y$ and $T{C}_z$ terms between the two cases, we examine the averaged Reynolds stress and downwind velocity on the control surface from the two cases in Figure 12 and Figure 13, in which different terms are defined as follows:

$$U_{side}=\left(\left({\int }_{z_h-R}^{z_h+R}\langle u\rangle dz\right){|}_{y_t-R}+\left({\int }_{z_h-R}^{z_h+R}\langle u\rangle dz\right){|}_{y_t+R}\right)/2,$$

(16)

$$U_{bot}=\left({\int }_{y_t-R}^{y_t+R}\langle u\rangle dy\right){|}_{z_h-R},$$

(17)

$$U_{top}=\left({\int }_{y_t-R}^{y_t+R}\langle u\rangle dy\right){|}_{z_h+R},$$

(18)

$${\langle u^{\prime }{v}^{\prime }\rangle }_{side}=\left(\left({\int }_{z_h-R}^{z_h+R}\langle u^{\prime }{v}^{\prime }\rangle dz\right){|}_{y_t-R}+\left({\int }_{z_h-R}^{z_h+R}\langle u^{\prime }{v}^{\prime }\rangle dz\right){|}_{y_t+R}\right)/2,$$

(19)

$${\langle u^{\prime }{w}^{\prime }\rangle }_{bot}=\left({\int }_{y_t-R}^{y_t+R}\langle u^{\prime }{w}^{\prime }\rangle dy\right){|}_{z_h-R},$$

(20)

$${\langle u^{\prime }{w}^{\prime }\rangle }_{top}=\left({\int }_{y_t-R}^{y_t+R}\langle u^{\prime }{w}^{\prime }\rangle dy\right){|}_{z_h+R}.$$

(21)

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 ${\langle u^{\prime }{v}^{\prime }\rangle }_{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 $T{C}_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 ${\langle u^{\prime }{w}^{\prime }\rangle }_{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 ${\langle u^{\prime }{w}^{\prime }\rangle }_{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 $T{C}_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.

5. 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.