*3.2. Case Setup*

The sizes of the computational domain are *Lx* × *Ly* × *Lz* = 1.20*δ* × 2.09*δ* × *δ* and *Lx* × *Ly* × *Lz* = 0.405*δ* × 2.09*δ* × *δ* with the corresponding numbers of grid nodes *Nx* × *Ny* × *Nz* = 1501 × 361 × 282 and *Nx* × *Ny* × *Nz* = 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 (*Ly*) and the height (*Lz*) 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 (*Lx*) is 15 rotor diameters for both cases. The turbine is located 2*D* 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* − *yt*| < 0.75*D* and *z* < 2.2*D*, where *yt* 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 <sup>Δ</sup>*tUh*/*<sup>D</sup>* = 7.9 × <sup>10</sup>−<sup>4</sup> and 5.9 × <sup>10</sup>−<sup>4</sup> for T80 and T27 cases, respectively, where *Uh* and *D* are their own incoming downwind velocity at hub height and diameter, respectively.

The roughness length of the ground is *ko* = 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 *Lx* × *Ly* × *Lz* = 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*<sup>1</sup> = *δ*/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 <sup>Δ</sup>*tUh*/*<sup>δ</sup>* = <sup>8</sup> × <sup>10</sup>−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 1*D* 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.

**Figure 3.** Statistics of inflow employed in the simulation for (**a**) mean downwind velocity and (**b**) standard deviation of velocity fluctuations, where *Ub* is the bulk velocity.

#### **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 Δ*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

$$
\Delta \mathcal{U} = \mathcal{U}\_{\text{in}} - \mathcal{U} \tag{3}
$$

where *U* is the time-averaged downwind velocity at different downwind locations, *Uin* is the time-averaged incoming downwind velocity (which is taken at 1*D* upwind the turbine). As in [33], a slight offset of 0.1*D* 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.

**Figure 4.** Comparison of the simulated velocity deficit profiles with measurements at the SWiFT site. 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 2*D* 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.

**Figure 5.** Instantaneous downwind velocity on the x-z plane passing through the rotor center for (**a**) T80 and (**b**) T27, respectively.

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.

**Figure 6.** Time-averaged downwind velocity on the x-z plane passing through the rotor center for (**a**) T80 and (**b**) T27, respectively.

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 7*D* downwind the turbine for the T80 case, which only persists to about 4*D* turbine downwind for the T27 case.

**Figure 7.** Turbulence kinetic energy (TKE) on the x-z plane passing through the rotor center for (**a**) T80 and (**b**) T27, respectively.

*Energies* **2020**, *13*, 3004

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

$$
\Delta k = k - k\_{\rm in} \tag{4}
$$

where *k* is the TKE at different downwind locations, *kin* is the TKE of the inflow (which is taken at 1*D* 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.

**Figure 8.** (**a**) Vertical profiles of time-averaged downwind velocity at 1*D* upwind of the turbine; (**b**–**f**) Vertical profiles of time-averaged velocity deficit at different locations downwind of the turbine. Red lines: T80; Blue lines: T27.

In Figure 9 we compare the turbine-added TKE from the two cases. First, we examine the normalized TKE at the 1*D* 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.

**Figure 9.** (**a**) Vertical profiles of turbulence kinetic energy at 1*D* upwind of the turbine; (**b**–**f**) Vertical profiles of turbine-added turbulence kinetic energy at different locations downwind of the turbine. Red lines: T80; Blue lines: T27.

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,\tag{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:

$$\mathcal{MC} = -\int\_{y\_l - R}^{y\_l + 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} dz dy \tag{6}$$

$$PT = -\int\_{y\_t - R}^{y\_l + R} \int\_{z\_h - R}^{z\_h + R} \frac{\partial \left( \left< p \right> \left< u\_j \right> / \rho \right)}{\partial x\_j} dz dy \tag{7}$$

*Energies* **2020**, *13*, 3004

$$T\mathbb{C} = -\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} dz dy \tag{8}$$

$$DF = 2\int\_{y\_l - R}^{y\_l + 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} dz dy \tag{9}$$

$$TP = 2\int\_{y\_l - R}^{y\_l + R} \int\_{z\_h - R}^{z\_h + R} \langle u\_i' u\_j' \rangle \frac{\partial \langle u\_i \rangle}{\partial x\_j} dz dy \tag{10}$$

$$DP = -2\int\_{y\_l-R}^{y\_l+R} \int\_{z\_h-R}^{z\_h+R} \left(\nu + \nu\_l\right) S\_{i\bar{j}} \frac{\partial \langle u\_i \rangle}{\partial x\_{\bar{j}}} dz dy \tag{11}$$

where *yt* is the coordinate of the rotor center in the lateral direction, *<sup>R</sup>* is the rotor radius, *Sij* <sup>=</sup> *∂ui*/*∂xj* + *∂uj*/*∂xi* /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 = TC\_x + TC\_y + TC\_{z\prime} \tag{12}$$

where

$$T\mathbb{C}\_{\mathbf{x}} = -\int\_{y\_l - R}^{y\_l + R} \int\_{z\_h - R}^{z\_h + R} \frac{\partial \left( \langle u'\_l u' \rangle \langle u\_l \rangle \right)}{\partial \mathbf{x}} d\mathbf{x} dy\_r \tag{13}$$

$$T\mathbb{C}\_{\mathcal{Y}} = \left( \int\_{z\_h - R}^{z\_h + R} \left( - \langle u\_i' v' \rangle \langle u\_i \rangle \right) dz \right) \Big|\_{y\_t - R}^{y\_l + R} \tag{14}$$

$$T\mathbb{C}\_z = \left( \int\_{y\_l - R}^{y\_l + R} \left( - \langle u\_i' w' \rangle \langle u\_i \rangle \right) dy \right) \Big|\_{z\_h - R}^{z\_h + R}. \tag{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 *TCx* term when *x*/*D* < 4 and *x*/*D* < 2 for the T80 and T27 cases, respectively. At further downwind locations the effect of the *TCx* term on the MKE budget is negligible. On the other hand, the other two components of the turbulence convection terms *TCy* and *TCz* contributes positively to the MKE budget at almost all downwind locations except in the region immediately downwind the turbine. The *TCy* term from the T27 case is higher than that from the T80 case for all downwind locations. On the other hand, the *TCz* 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.

**Figure 10.** Budget of mean kinetic energy for (**a**) T80 and (**b**) T27, respectively. Line with crosses: *MC* term (Equation (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 <sup>1</sup> <sup>2</sup> *DU*<sup>3</sup> *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)), respectively. Red line: T80; Blue line: T27. Different terms are normalized using <sup>1</sup> <sup>2</sup> *DU*<sup>3</sup> *h*.

To better understand the differences in *TCy* and *TCz* 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:

$$\mathcal{U}\_{\rm side} = \left( \left( \int\_{z\_h - R}^{z\_h + R} \langle u \rangle dz \right) \Big|\_{y\_t - R} + \left( \int\_{z\_h - R}^{z\_h + R} \langle u \rangle dz \right) \Big|\_{y\_t + R} \right) / 2,\tag{16}$$

$$
\Delta I\_{b\nu t} = \left( \int\_{y\_t - R}^{y\_t + R} \langle \mu \rangle d\mathcal{y} \right) \Big|\_{z\_h - R} \tag{17}
$$

$$
\Delta I\_{top} = \left( \int\_{y\_t - \mathbb{R}}^{y\_t + \mathbb{R}} \langle u \rangle dy \right) \Big|\_{z\_h + \mathbb{R}'} \tag{18}
$$

$$
\langle \langle \boldsymbol{u}' \boldsymbol{v}' \rangle\_{\text{side}} = \left( \left( \int\_{z\_h - R}^{z\_h + R} \langle \boldsymbol{u}' \boldsymbol{v}' \rangle d\boldsymbol{z} \right) \Big|\_{y\_l - R} + \left( \int\_{z\_h - R}^{z\_h + R} \langle \boldsymbol{u}' \boldsymbol{v}' \rangle d\boldsymbol{z} \right) \Big|\_{y\_l + R} \right) / 2,\tag{19}
$$

$$
\langle \langle u'w' \rangle\_{\rm bot} = \left( \int\_{y\_t - R}^{y\_t + R} \langle u'w' \rangle dy \right) \Big|\_{z\_h - R} \tag{20}
$$

$$
\langle \langle u'w' \rangle\_{top} = \left( \int\_{y\_t - R}^{y\_t + R} \langle u'w' \rangle dy \right) \Big|\_{z\_k + R}. \tag{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 *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 *TCy* 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 2*D* to 5*D* 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 *TCz* 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.

**Figure 12.** The averaged downwind velocity on the surface of the control volume for (**a**) on the two side surface of the control volume (Equation (16)), (**b**) on the bottom surface of the control volume Equation (17) and (**c**) on the top surface of the control volume Equation (18).

**Figure 13.** The averaged Reynolds stress on the surface of the control volume for (**a**) on the two side surface of the control volume (Equation (19)), (**b**) on the bottom surface of the control volume Equation (20) and (**c**) on the top surface of the control volume Equation (21).

**Figure 14.** The turbulence production (TP) term (Equation (10)). Red line: T80; Blue line: T27. The TP term is normalized using <sup>1</sup> <sup>2</sup> *DU*<sup>3</sup> *h*.
