2.1. Flow solver

The turbulent flow was solved using the LES module of the Virtual Flow Simulator code (VFS-Wind). The capability of the employed code for simulating turbine wakes has been extensively validated using wind tunnel (Yang et al. Reference Yang, Sotiropoulos, Conzemius, Wachtler and Strong2015b) and field measurements (Yang, Pakula & Sotiropoulos Reference Yang, Pakula and Sotiropoulos2018; Yang et al. Reference Yang, Milliren, Kistner, Hogg, Marr, Shen and Sotiropoulos2021). The code is open source and available at GitHub. (https://github.com/SAFL-CFD-Lab/VFS-Wind.) The governing equations are the filtered incompressible Navier–Stokes equations in curvilinear coordinates shown as

(2.1)$$\begin{gather} J \frac{\partial U^{i} }{\partial \xi^{i} } =0, \end{gather}$$

(2.2)$$\begin{gather}\frac{1}{J}\frac{\partial U^{i}}{\partial t} = \frac{\xi^{i}_l}{J}\left( -\frac{\partial }{\partial \xi^{j}} \left(U^{j}u_l\right)+\frac{\mu}{\rho}\frac{\partial }{\partial \xi^{j}}\left(\frac{g^{jk}}{J}\frac{\partial u_l}{\partial \xi^{k}}\right) -\frac{1}{\rho}\frac{\partial}{\partial \xi^{j}}\left( \frac{\xi^{j}_l p}{J}\right)-\frac{1}{\rho}\frac{\partial \tau_{lj}}{\partial \xi^{j}} +f_l\right), \end{gather}$$

where $i,j,k,l=\{1,2,3\}$ are the tensor indices and $\xi ^{i}$ are the curvilinear coordinates related to the Cartesian coordinates $x_l$ by the transformation metrics $\xi ^{i}_l = \partial \xi ^{i}/\partial x_l$. Additionally, $J$ denotes the Jacobian of the geometric transformation, $U^{i} = (\xi ^{i}_l / J)u_l$ is the contravariant volume flux with $u_l$ as the velocity in Cartesian coordinates, $\mu$ denotes the dynamic viscosity, $\rho$ is the fluid density, $g^{jk} = \xi ^{j}_l\xi ^{k}_l$ are the components of the contravariant metric tensor, $p$ is the pressure and $f_l$ are body forces introduced by the actuator-type wind turbine model. In the momentum equation, $\tau _{ij}$ is the subgrid-scale (SGS) stress, which is computed using the dynamic Smagorinsky model (Smagorinsky Reference Smagorinsky1963; Germano et al. Reference Germano, Piomelli, Moin and Cabot1991). The governing equations are discretized on a structured curvilinear grid. A second-order accurate central differencing scheme is used for the spatial discretization with a second-order fractional step method (Ge & Sotiropoulos Reference Ge and Sotiropoulos2007) for the temporal integration. The momentum equation is solved with the matrix-free Newton–Krylov method (Knoll & Keyes Reference Knoll and Keyes2004). The pressure Poisson equation, for satisfying the continuity constraint, is solved using the generalized minimal residual (GMRES) method with an algebraic multi-grid acceleration (Saad Reference Saad1993).

2.2. Wind turbine parametrization method

Geometry-resolving simulations of wind turbine wakes are extremely expensive, because the characteristic length of wind turbine wakes is often more than two orders of magnitude larger than the thickness of the boundary layer over the blades. (The length scale of wind turbine wakes is proportional to the rotor's diameter $D$; the thickness of the boundary layer on the blade is approximately 1 cm for a turbine of $D \approx 100$ m operating in region II, i.e. when the wind velocity is between the cut-in and the rated speeds Yang & Sotiropoulos Reference Yang and Sotiropoulos2018) For this reason, the aerodynamics of wind turbines is often parametrized using actuator disk (Chattot Reference Chattot2014), actuator line (Sorensen & Shen Reference Sorensen and Shen2002) or actuator surface (Shen, Zhang & Sorensen Reference Shen, Zhang and Sorensen2009; Yang & Sotiropoulos Reference Yang and Sotiropoulos2018) models. In this work, the class of well validated actuator surface (AS) models for turbine blades and nacelle proposed by Yang & Sotiropoulos (Reference Yang and Sotiropoulos2018) was employed. The AS model represents each rotor blade with a simplified two-dimensional (2-D) surface defined by the chord length and the twist angle at each radial location. The lift and drag forces $\boldsymbol {L}$ and $\boldsymbol {D}$ at each radial location were determined using the tabulated airfoil data using the local instantaneous relative incoming velocity as follows:

(2.3)\begin{equation} \boldsymbol{L} = \tfrac{1}{2}\rho C_{L}c|V_{ref}|^{2} \boldsymbol{e_{L}} \end{equation}

and

(2.4)\begin{equation} \boldsymbol{D} = \tfrac{1}{2}\rho C_{D}c|V_{ref}|^{2} \boldsymbol{e_{D}}, \end{equation}

where $c$ is the chord length, $V_{ref}$ is the instantaneous incoming velocity relative to the rotating blade at the computing point averaged over the chord length, $\boldsymbol {e}_{L}$ and $\boldsymbol {e}_{D}$ are the unit directional vectors of the lift and drag forces, and $C_{L}$ and $C_{D}$ are the lift and the drag coefficients defined in 2-D airfoil tables as a function of Reynolds number and the angle of attack. Corrections including the three-dimensional stall delay model (Du & Selig Reference Du and Selig1998) and the tip loss correction (Shen et al. Reference Shen, Mikkelsen, Sørensen and Bak2005) were applied. With the computed $\boldsymbol {L}$ and $\boldsymbol {D}$, the body force in (2.2) was calculated by uniformly distributing the forces along the chord as follows:

(2.5)\begin{equation} \boldsymbol{f} = (\boldsymbol{L}+\boldsymbol{D})/c. \end{equation}

A smoothed discrete delta function (Yang et al. Reference Yang, Zhang, Li and He2009) was employed for transferring quantities between the actuator surface and the background grid nodes for solving the flow.

2.3. Simulation set-up

In this section, we present the set-up of the cases simulated in this work, as summarized in table 1. Four different yaw angles $\gamma$, i.e. $\gamma = 0^{\circ }, 10^{\circ }, 20^{\circ }, 30^{\circ }$ (one extra case of $\gamma = -30^{\circ }$ was also conducted to check the influence of negative yaw angle in Appendix D.), defined as the misalignment between the inflow and the rotor axis (shown in figure 1), were considered for three turbulent inflows with different levels of turbulence intensity. For inflow 1, the effects of different turbine operation conditions (tip speed ratios) were examined for different fixed yaw angles. For inflows 2 and 3, simulations were conducted under the same turbine operating conditions to investigate the influence of the inflow turbulence on the wake. For each inflow, an extra case without the wind turbine was also simulated to provide the reference turbulent boundary layer flow. All the simulations were conducted with the same set-up, i.e. the same boundary and initial conditions, and the same time step to ensure the inflow turbulent structures were synchronized for the simulated cases to facilitate the one-to-one comparison of instantaneous flow fields. Details on the employed wind turbine, the computational domain and the boundary conditions, and the turbulent inflow generation are provided in the following subsections.

Figure 1. Schematic of the computational set-up employed for simulating the wake of a yawed wind turbine.

Table 1. Summary of simulation cases: $\gamma$ denotes the yaw angle, $\tilde {\lambda }$ denotes the TSR defined with the inflow velocity projected along the rotor's axis.

2.3.1. Wind turbine and operation conditions

We simulated the three-blade horizontal-axis Clipper Liberty 2.5 MW wind turbine located at the EOLOS wind energy research field station at the University of Minnesota, USA. The rotor diameter was $D=96$ m, the hub was at $z_{hub} = 80$ m, and the nacelle had a cuboid-like shape with its dimensions being approximately equal to $5.3 \ \text {m} \times 4.7 \ \text {m} \times 5.5 \ \text {m}$. The tower was conical with diameters varying from 3.0 m at the top to 4.1 m at the bottom. It should be noted that the tower was not modelled in this work, which may affect the top-down asymmetry of the wake of a yawed wind turbine as reported by Howland et al. (Reference Howland, Bossuyt, Martínez-Tossas, Meyers and Meneveau2016). More information about the EOLOS wind turbine can be found in the papers by Hong et al. (Reference Hong, Toloui, Chamorro, Guala, Howard, Riley, Tucker and Sotiropoulos2014) and Chamorro et al. (Reference Chamorro, Lee, Olsen, Milliren, Marr, Arndt and Sotiropoulos2015).

In each simulation, the rotor rotated at a fixed modified TSR with respect to the project velocity $U_d$ normal to the rotor sweeping plane, as follows:

(2.6)\begin{equation} \tilde{\lambda} = \frac{\varOmega R}{U_d}, \end{equation}

where $\varOmega$ is the angular velocity of the rotor, $R = 48$ m is the radius of the rotor and

(2.7)\begin{equation} {U_d} = U_\infty \cos \gamma \end{equation}

is the inflow velocity projected in the rotor's axis direction.

The scenario we study herein is ideal, where the rotor rotates based on the modified TSR $\tilde {\lambda }$, which cannot be easily realized in real-life applications. Defining the TSR based on the velocity projected to the axial direction of the rotor was for the convenience of explaining how we control the rotational speed of the rotor in this work. It does not guarantee optimal energy conversion of a yawed wind turbine (Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2019). Three different modified TSRs $\tilde {\lambda } = \{7,8,9\}$ were employed to vary the thrust and power coefficients to investigate the effects of operation conditions on wake characteristics of yawed wind turbines.

2.3.3. Turbulent inflow

The turbulent inflows were generated from precursory LESs with a larger computational domain ($L^{'}_x \times L^{'}_y \times L^{'}_z = 62D \times 46D \times 1\ \text {km}$) to capture large-scale eddies in the atmospheric boundary layer (Wang & Zheng Reference Wang and Zheng2016; Liu, Wang & Zheng Reference Liu, Wang and Zheng2019). Periodic boundary conditions were applied in the horizontal directions. The top boundary was set as free-slip. The lower boundary was set as no-slip and modelled using the logarithmic law for rough walls, as in the turbine simulations. Three different ground roughness lengths were considered, i.e. $z_0 = 0.00016\ \text {m},0.001\ \text {m}, 0.01\ \text {m}$ for inflow 1, inflow 2 and inflow 3, respectively. The characteristics of the inflows are shown in figure 2. In the precursory simulations, the velocity fields on a $y-z$ plane were saved at each time step and then applied at the inlet of the turbine simulations. Linear interpolations were carried out if necessary.

Figure 2. Characteristics of the turbulent inflows: (a) the time-averaged streamwise velocity profile; and (b) the streamwise turbulence intensity profile. The turbulence intensity at hub height for inflows 1, 2 and 3 are 6.9 %, 7.9 % and 9.6 %, respectively. The horizontal dash–dotted lines represent the turbine hub height and the dotted lines represent the upper and the lower limits of the rotor.

3.1. Variation of the thrust and power with yaw angles

Before probing into the flow fields, the variations of the computed thrust (normal to the rotor sweep plane) and power (calculated from the shaft torque) with yaw angles are compared in figure 3. The results were normalized by the value at $\gamma = 0^{\circ }$ at each corresponding modified TSR $\tilde {\lambda }$. All results agreed well with the prediction obtained using the axial momentum theory (Burton et al. Reference Burton, Jenkins, Sharpe and Bossanyi2011), which states that the thrust $T$ and power $P$ scale with $\cos ^{2} (\gamma )$ and $\cos ^{3} (\gamma )$, respectively. This relation has also been observed in other numerical and experimental studies (Krogstad & Adaramola Reference Krogstad and Adaramola2012; Bartl et al. Reference Bartl, Mühle, Schottler, Sætran, Peinke, Adaramola and Hölling2018). With this relation, the thrust and power coefficients ($\tilde {C}_T$ and $\tilde {C}_P$), which are defined with respect to $U_d$ and shown as

(3.1)$$\begin{gather} \tilde{C}_T = \frac{T}{\frac{1}{2}\rho A U_d^{2}} =\frac{T}{\frac{1}{2}\rho A \left(U_\infty\cos\gamma\right)^{2}} , \end{gather}$$

(3.2)$$\begin{gather}\tilde{C}_P = \frac{P}{\frac{1}{2}\rho A U_d^{3}}= \frac{P}{\frac{1}{2}\rho A \left(U_\infty\cos\gamma\right)^{3}}, \end{gather}$$

are independent of the yaw angles and are equal to those in the non-yawed cases, because of $T \propto \cos {\gamma }^{2}$ and $P \propto \cos {\gamma }^{3}$. The values of $\tilde {C}_T$ and $\tilde {C}_P$ for different TSR are shown in table 2. It is worth noting that the yaw angle independence of $\tilde {C}_T$ and $\tilde {C}_P$ may result from the control strategy employed in the present work, i.e. the modified TSR $\tilde {\lambda }$ was fixed. However, if another control strategy is employed, the thrust $T$ and the power $P$ may follow $\cos ^{a} \gamma$ and $\cos ^{b} \gamma$, respectively, with $a<2$ and $b<3$ (Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2019; Howland et al. Reference Howland2020; Liew, Urbán & Andersen Reference Liew, Urbán and Andersen2020). Making the comparison in figure 3 is just to verify the implementation of this idealized rotor control rather than indicating that the power and thrust should vary in that way in real-life applications.

Figure 3. Influence of yaw angle $\gamma$ on the rotor's thrust and power: (a) the normalized thrust $T(\gamma )/T(0)$ (symbols) compared with $\cos ^{2} \gamma$ (solid line); and (b) the normalized power $P(\gamma )/P(0)$ (symbols) compared with $\cos ^{3} \gamma$ (solid line). Here, $T(0)$ and $P(0)$ are the thrust and the power at $\gamma =0$, respectively.

Table 2. Modified thrust and power coefficients at different modified TSRs.

3.2. Time-averaged velocity fields

In this subsection, we examine the influence of yaw angle on the time-averaged wake velocity. To illustrate such an influence, we first plot, in figure 4, the streamwise and the transverse velocities $U$ and $V$ on the hub-height plane for cases with $\tilde {\lambda } = 8$ (the cases with other $\tilde {\lambda }$ have a similar pattern). The red thick dashed lines and the red dotted lines plot the wake centreline ($y = Y_C(x$)) and the boundary ($y = Y_C(x) \pm R_{1/2}(x$)) obtained by fitting the time-averaged velocity deficit using the Gaussian fit at each downstream location ($x$), as follows,

(3.3)\begin{equation} \Delta U(x,y) = U_\infty - U(x,y) = \Delta U_C(x) \exp{\left(-\frac{(y-Y_C(x))^{2}}{2 S^{2}(x)}\right)}, \end{equation}

where $U(x,y)$ is the streamwise velocity and $S(x)$ is the standard deviation of the Gaussian distribution to be fitted with the velocity profile. A schematic of these definitions is illustrated in figure 5. The wake half-width is defined as $R_{1/2}(x) = \sqrt {2\ln 2} S(x)$, which gives the distance from $Y_C(x)$ to the position where $\Delta U = (\Delta U_C)/2$.

Figure 4. The time-averaged velocity field behind wind turbines (at hub-height plane $z=z_{hub}$): (a,c,e,g) streamwise velocity; and (b,d,f,h) the transverse velocity. Panels in each column are ordered by increasing yaw angle . The red dashed lines represent the centrelines of the streamwise velocity deficit. The red dots denote the wake width where the streamwise velocity deficit is half that on the centreline. The blue dashed line and blue dots denote the wake centreline and width defined by the transverse velocity. The black solid lines illustrate the rotors.

Figure 5. Schematic and and key variables for describing the wake behind a yawed wind turbine on the hub-height plane, where $U_\infty$ is the inflow velocity, $U_d$ is the inflow velocity component in the rotor's axis direction, $\gamma$ is the yaw angle, $T$ is the thrust on the rotor and $R = 0.5D$ is the rotor's radius. Here, $U$ and $V$ denote the streamwise and the transverse velocity, respectively.

First, figure 4 shows that the streamwise velocity $U$ behaves similarly in the yawed and non-yawed cases. As seen, $U$ is symmetric with respect to the wake centreline. The overall patterns of $U$ computed at different yaw angles are very similar, although the magnitude of the velocity deficit ($\Delta U$) and the wake width ($R_{1/2}$) gradually reduce with the increase of yaw angle owing to the reduction of the streamwise thrust on the rotor $T_x = T \cos {\gamma }$ and the projected rotor width $\tilde {R} = R \cos {\gamma }$.

In contrast, the transverse velocity $V$ behaves differently in yawed and non-yawed cases, as seen in the second column of figure 4. In the region close to the rotor, a diverging transverse velocity is observed for all yaw angles, which is related to the expansion of the streamtube at the rotor. In this region, it is difficult to distinguish the yaw-induced transverse velocity and that induced by the wake expansion. Differences are observed in the far wake: in the non-yawed case, the transverse velocity is close to zero and is symmetrical to the wake centreline; when the rotor yaws, on the other hand, a transverse velocity in the negative $y$ direction is observed and its magnitude increases with yaw angle. More importantly, the transverse velocity is observed as being asymmetric with respect to the wake centreline $Y_C(x)$ defined by the streamwise velocity deficit. That the transverse velocity resides only on one side of the wake centreline (behind the leading edge of the yawed rotor) is in agreement with the wind tunnel measurement of Bastankhah & Porté-Agel (Reference Bastankhah and Porté-Agel2016). Owing to these differences observed in the transverse velocity $V$ as compared with the streamwise velocity $U$, it is expected that the wake quantities defined by $U$ and $V$ should be treated differently. For this reason, we plot in addition the wake centre position $Y_C^{V}(x)$ and the wake width of $R_{1/2}^{V}(x)$, where the transverse velocity $V$ is reduced to half of the maximum. As seen, $Y_C^{V}(x)$ (the blue dash lines) approximately overlaps with the streamwise wake boundary (the red dots).

In the following subsections, we will examine the similarity of these wake characteristics, using the different characteristic length and velocity scales presented in Appendix A.

3.3. Similarity of wake centreline $Y_C(x)$

In figure 6, we examine the similarity of the wake centreline $Y_C(x)$ for different yaw angles. The length scale employed to normalize $Y_C(x)$ is

(3.4)\begin{equation} Y_N = D \tilde{C}_T \cos^{2} \gamma \sin \gamma, \end{equation}

which is proportional to the transverse thrust $T_y$ as shown by Jiménez et al. (Reference Jiménez, Crespo and Migoya2010) and Shapiro et al. (Reference Shapiro, Gayme and Meneveau2018). Figure 6(a) shows that the wake centreline deflections increased with yaw angle. For the same yaw angle, it was observed that $Y_C(x)$ increases with $\tilde {\lambda }$ at the far wake especially for the $\gamma =30^{\circ }$ case. However, the differences of $Y_C$ between cases of different $\tilde {\lambda }$ are small because of the relative small difference in the thrust coefficients $C_T$. In figure 6(b), it is observed that the wake centrelines of different yaw angles overlap with each other after being normalized by $Y_N$. This good scaling shows that the length scale ((3.4)) well describes the similarity of the wake centreline despite the complex dynamics of yawed turbine wakes. Once the wake centreline is known for one yaw angle, it can be generalized with the proposed length scale $Y_N$ to predict the wake centrelines at other yaw angles if the similarity observed in this work also exists in realistic conditions taking into account the effects of thermal stability, wind speed shear, wind direction, wake interaction, etc.

Figure 6. Wake centreline $Y_C(x)$: (a) normalized by the rotor's diameter $D$; and (b) normalized by the length scale $Y_N = D \tilde {C}_T\cos ^{2} \gamma \sin \gamma$.

3.4. Similarity of streamwise velocity deficit $\Delta U$

In this section, we examine the similarity of the quantities related to the streamwise velocity deficit. The characteristic velocity employed for scaling is

(3.5)\begin{equation} U_N = U_\infty \left(1-\sqrt{1-\tilde{C}_T\cos^{2} \gamma }\right), \end{equation}

which represents the streamwise velocity reduction caused by the rotor (Burton et al. Reference Burton, Jenkins, Sharpe and Bossanyi2011; Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2016; Shapiro et al. Reference Shapiro, Gayme and Meneveau2018). In figure 7, we first examine the streamwise evolution of the velocity deficits $\Delta U_C(x)$ along the wake centreline for different yaw angles and different $\tilde {\lambda }$. Figure 7(a) shows that $\Delta U_C(x)$ with different $\tilde {\lambda }$ and different yaw angles deviates from each other, which reveals the influence of both factors. For the same $\tilde {\lambda }$, $\Delta U_C(x)$ decreases with the increase of yaw angle. For the same yaw angle, $\Delta U_C(x)$ increases when increasing $\tilde {\lambda }$ owing to the increase of $\tilde {C}_T$. More importantly, at all considered downstream locations ($1D\le x \le 10D$), the $\Delta U_C(x)$ profiles are observed to vary in a similar manner. This similarity is further confirmed in figure 7(b), which shows the $\Delta U_C(x)$ profiles overlap with each other when normalized using the characteristic velocity $U_N$. In the figure, it is also noticed that $\Delta U_C(x)$ first increases until approximately $x = 2.5D$ and then gradually decreases to far-wake locations. This flow deceleration is related to the recovery of the pressure. In the far wake, the scaled wake centre velocities $\Delta U_C(x)/U_N$ show an excellent agreement with each other, which indicates $U_N$ is the proper velocity scale for the velocity deficit.

Figure 7. The characteristic velocity deficit $\Delta U_C(x)$ for different $\tilde {\lambda }$ and yaw angles: (a) normalized by the inflow velocity $U_\infty$; and (b) normalized by the velocity scale as $U_N = U_\infty (1-\sqrt {1-\tilde {C}_T\cos ^{2} \gamma })$.

We then examined the wake half-width $R_{1/2}(x) = \sqrt {2\ln {2}}S(x)$, where $S(x)$ is obtained by fitting the velocity deficit to the Gaussian function ((3.3)). The initial wake width $R_N$, derived using the streamwise one-dimensional momentum theory (details can be found in Appendix A), was employed for the normalization. The expression of $R_N$ is shown as follows:

(3.6)\begin{equation} R_N = R \cos{\gamma}\sqrt{\frac{1+\sqrt{1-\tilde{C}_T\cos^{2} \gamma}}{2\sqrt{1-\tilde{C}_T\cos^{2} \gamma}}}, \end{equation}

which is a function of the thrust coefficient $\tilde {C}_T$, yaw angle $\gamma$ and the rotor's radius $R$. Figure 8(a) shows that there were differences in $R_{1/2}/R$ for the cases with different yaw angles and different $\tilde {\lambda }$, but the curves overlapped when normalizing $R_{1/2}$ using $R_N$, as shown in figure 8(b). That the curves overlap with each other at all considered downstream locations in the range of $1D< x<10D$ shows that the proposed length scale $R_N$, a characteristic wake width scale in the near wake, is still valid in the far wake. It is also noticed that the wake width decreased slightly in the initial near-wake region, which was attributed to the difficulties in fitting the velocity deficit profiles to the Gaussian function when the root loss phenomenon is present, as shown in figure 4. It is noted that $R_N$ is different from that in the paper by Bastankhah & Porté-Agel (Reference Bastankhah and Porté-Agel2016) for predicting the wake width at the onset of far wake (see details in Appendix B).

Figure 8. The characteristic wake width $R_{1/2}(x)$ for different $\tilde {\lambda }$ and yaw angles: (a) normalized by the rotor radius $R$; and (b) normalized by the length scale $R_N =R \cos {\gamma }\sqrt {({1+\sqrt {1-\tilde {C}_T\cos ^{2} \gamma }})/({2\sqrt {1-\tilde {C}_T\cos ^{2} \gamma }})}$.

In figure 9, the similarity of the horizontal profiles of the streamwise velocity deficit at different turbine downstream locations is examined. Figure 9(a–d) shows the magnitude of the velocity deficit and the wake width decrease when increasing the yaw angles. In contrast, the transverse profiles of the streamwise velocity deficit shown in figure 9(e–h) overlap with each other when normalized using the velocity scale $U_N$ and the length scale $R_N$, and are shifted with respect to the wake centreline $Y_C(x)$. It is seen from figure 9(e–h) that the velocity deficit profiles are symmetrical to the wake centreline for different yaw angles including $\gamma = 30^{\circ }$. The characteristic wake width scale $R_N$ and velocity deficit scale $U_N$ are proper for describing the similarity of velocity deficits for different yaw angles. It is noticed that at $3D$ turbine downstream, the velocity deficits $\Delta U/ U_N$ on the plateau are approximately equal to 1, which indicates that the present scaling factor $U_N$ is a good estimation of the velocity deficit in the near wake. Furthermore, the overlapping of the velocity deficit profiles at further turbine downstream locations (figure 9f–h) shows that $U_N$ successfully describes the similarity of the velocity deficits for different yaw angles, although it is derived at the imminent near wake.

Figure 9. The time-averaged velocity deficit $\Delta U = U_\infty - U$ at downstream locations of ($x \in \{3D,5D,7D,9D\}$) on the hub-height plane ($z = z_{hub}$) for different $\tilde {\lambda }$ and yaw angles: (a–d) normalized by the inflow velocity $U_\infty$; and (e–h) the wake centres are shifted with respect to $Y_C$, with the values normalized by the proposed velocity and length scale $U_N$ and $R_N$.

3.5. Transverse velocity $V$

In this subsection, we investigate the similarity of the transverse velocity $V$. Recall that in figure 4, the transverse velocity field is shown to be significantly different from the streamwise velocity and is thus not expected to have the same characteristic velocity and length scales. Consequently, the characteristic quantities scales defined using the transverse component of the thrust $T_y$, following Jiménez et al. (Reference Jiménez, Crespo and Migoya2010) and Shapiro et al. (Reference Shapiro, Gayme and Meneveau2018), will be employed for scaling the quantities related to the transverse velocity, i.e. $R^{V}_{1/2}(x)$ and $V_C(x)$. Because the transverse velocity is induced by the CVP in the wake (Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2016; Howland et al. Reference Howland, Bossuyt, Martínez-Tossas, Meyers and Meneveau2016), the transverse profiles of $V$ are fitted with a CVP shape function derived based on the Biot–Savart law instead of the Gaussian function (see Appendix E for details).

We first examine the wake half-width $R_{1/2}^{V}(x)$ defined using the transverse profiles of the transverse velocity. Figure 10 shows the streamwise variation of $R^{V}_{1/2}$. Recall that both wake expansion and yaw can induce transverse velocity near the rotor, as shown in figure 4, so using the fitting function based on the CVP to define $R_{1/2}^{V}$ works only in the far wake $(x>4D)$. Despite some fluctuations, it is found that the variation of $R_{1/2}^{V}(x)$ with yaw angles is smaller than that of $R_{1/2}$. The curves of the same $\tilde {\lambda }$ cluster together, which shows that $R_{1/2}^{V}(x)$ is less dependent on the yaw angle. The yaw-independent scaling of $R_{1/2}^{V}$ is consistent with the theory of Glauert (Reference Glauert1926) for auto-gyros, where the initial influencing zone of $V$ is found to have the same size as the rotor disc instead of the projected size of the rotor disc (Burton et al. Reference Burton, Jenkins, Sharpe and Bossanyi2011). The fluctuations of $R_{1/2}^{V}$ may be explained by the complex dynamics related with the transverse velocity, e.g. the CVP and its interaction with the wake rotation (Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2016; Howland et al. Reference Howland, Bossuyt, Martínez-Tossas, Meyers and Meneveau2016). In the far-wake locations ($x > 4D$), both $R_{1/2}^{V}(x)$ and $R_{1/2}(x)$ gradually increase with the distance from the turbine.

Figure 10. Streamwise variations of the half-width $R^{V}_{1/2}(x)$ computed by fitting the transverse profiles of the transverse velocity $V$. The grey curves replot $R_{1/2}(x)$ computed from the streamwise velocity in figure 8(a) for comparison.

We then examined the similarity of the transverse velocity using the velocity scale $V_N$, which is proportional to the initial transverse velocity given by Jiménez et al. (Reference Jiménez, Crespo and Migoya2010) and Shapiro et al. (Reference Shapiro, Gayme and Meneveau2018):

(3.7)\begin{equation} V_N = \tfrac{1}{2}\tilde{C}_T \cos^{2} \gamma \sin \gamma U_\infty. \end{equation}

Figure 11 shows the streamwise variations of the maximum transverse velocity $V_C(x)$ obtained from with the CVP fitting function. In figure 11(b), $V_C(x)/V_N$ increases from approximately $0.4$ at $x = 2D$ to $0.7$ at $x=4D$ and then gradually decreases to approximately $0.4$ at $x = 10D$. More importantly, the $V_C/V_N$ profiles for different $\gamma$ and $\tilde {\lambda }$ approximately overlap with each other especially at far downstream locations, which shows the characteristic velocity $V_N$ successfully captures the main feature of the similarity. (The standard deviations of the relative differences of the different curves from the mean of the curves shown in figure 11(a,b) are approximately 31 % and 3.5 %, respectively, in the far-wake region ($x>4D$).) However, there are slightly larger differences between these curves when comparing with $U/U_N$ (figure 7), which implies that $V_N$ does not include all the underlying physics.

Figure 11. The characteristic transverse velocity $V_C(x)$ for different yaw angles: (a) normalized by inflow velocity $U_\infty$; and (b) normalized by the proposed velocity scale $V_N = (\tilde {C}_T \cos ^{2}{\gamma } \sin {\gamma } U_\infty$)/2.

In figure 12(a–d), the $V$ profiles are normalized using the incoming wind speed $U_{\infty }$. For a given $\tilde {\lambda }$, the velocity amplitude increases when increasing the yaw angle for all the considered streamwise locations. At $x= 3D$, the $V$ profiles show complex variations, especially for the cases with a yaw angle $\gamma = 10^{\circ }$. In figure 12(e–h), the $V$ profiles are normalized by the proposed velocity scale $V_N$ and shifted with respect to the centre location $Y_C$ based on the streamwise velocity deficit. It is noted that scaling of the transverse velocity profiles is not as accurate as that for the streamwise velocity deficit. This is probably because many other factors including the wake rotation and vertical shear, which affect the transverse velocity, are not taken into account in the proposed scaling. However, the profiles of the transverse velocity are approximately symmetrical about the peak, and the overall shape can be fitted with the proposed CVP function (the grey curves) especially in the far wake. Additionally, the peak of the transverse velocity locates on one side of the wake centreline $Y_C$ at approximately $0.5D$ and varies slightly with the downstream locations, which was assumed to be approximately equal to $S$ in (3.3) in the paper by Bastankhah & Porté-Agel (Reference Bastankhah and Porté-Agel2016). This shift of the maximum of the transverse velocity from the wake centreline may explain the overestimation of $Y_C$ computed using a uniform $V$, as in the analytical wake models of Jiménez et al. (Reference Jiménez, Crespo and Migoya2010), as shown by e.g. Bastankhah & Porté-Agel (Reference Bastankhah and Porté-Agel2016). However, it should also be noticed that the shape of $V$ is different from the idealized shape of the CVP fitting function in the near wake, as shown in figure 13, as a result of complex flow dynamics in the near wake caused by the nacelle, the wake rotation (Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2016; Howland et al. Reference Howland, Bossuyt, Martínez-Tossas, Meyers and Meneveau2016) and the inflow vertical shear (Gebraad et al. Reference Gebraad, Thomas, Ning, Fleming and Dykes2017). The tower, which was not considered in the work, may also have an effect.

Figure 12. The time-averaged transverse velocity $V$ profile at downstream locations of ($x \in \{3D,5D,7D,9D\}$) on the hub-height plane ($z = z_{hub}$) for different yaw angles: (a–d) $V$ normalized by the inflow velocity $U_\infty$; (e–h) the wake centres are shifted with respect to $Y_C$, with the values normalized by the velocity scale $V_N = (\tilde {C}_T \cos ^{2} \gamma \sin \gamma U_\infty )/2$. The grey thick curves show the corresponding fitting function.

Figure 13. Wake characteristics at $x= 7D$ for cases of different yaw angles with $\tilde {\lambda } = 8$: (a–d) the streamwise velocity deficit $\Delta U$; (e–h) the turbine-added turbulence kinetic energy $\Delta k$; (i–l) the turbine-added Reynolds shear stress $\Delta \langle u'v' \rangle$; and (m–p) the turbine-added Reynolds shear stress $\Delta \langle u'w' \rangle$. The dashed circles illustrate the projected rotor.

4.1. Turbine-added turbulent kinetic energy and Reynolds shear stresses

In this subsection, we examine the similarity of the turbine-added turbulence kinetic energy $\Delta k$ and the turbine-added Reynolds shear stresses $\Delta \langle u^{\prime } v^{\prime } \rangle$ and $\Delta \langle u^{\prime } w^{\prime } \rangle$. The employed characteristic velocity scale is given as follows:

(4.1)\begin{equation} U_T = U_\infty\cos\gamma\sqrt{\tilde{C}_T/2}. \end{equation}

Figure 13 shows the streamwise velocity deficit, the turbine-added turbulence kinetic energy and Reynolds shear stresses at $x=7D$ for cases of different yaw angles with $\tilde {\lambda }=8$. The streamwise velocity deficit is plotted in figure 13(a–d), to illustrate the wake region. The wake centre is shown to be deflected to the left and the wake's cross-section is deformed to a curled-shape with top-down asymmetry when increasing the yaw angle, as described in the literature (Howland et al. Reference Howland, Bossuyt, Martínez-Tossas, Meyers and Meneveau2016). Interestingly, the curled-shape is less obvious for the turbine-added turbulence kinetic energy $\Delta k$ , which shows an approximate left-right symmetry, as shown in figure 13(e–h). The turbine-added Reynolds shear stress $\Delta \langle u'v' \rangle$ is plotted in figure 13(i-l) and $\Delta \langle u'w' \rangle$ is plotted in figure 13(m–p). Overall it is found that these Reynolds stresses become asymmetric at large yaw angles. Moreover, it is found that $\Delta \langle u'v' \rangle > \Delta \langle u'w' \rangle$, which may be related to the stronger horizontal wake meandering. In the following, the transverse profiles of $\Delta k$ and $\Delta \langle u'v' \rangle$ at the hub height and different streamwise locations will be examined in detail for different yaw angles and turbine operation conditions. The vertical profiles $\Delta \langle u'w' \rangle$, however, will not be be examined because of its irregular shape and the difficulty to define the location for comparison.

Figure 14 compares the turbine-added turbulence kinetic energy $\Delta k$ computed from cases with different $\gamma$ and $\tilde {\lambda }$ at different downstream locations. The curves in figure 14(a–d) are $\Delta k$ profiles normalized using $U^{2}_{\infty }$. Here, the $\Delta k$ profiles of different yaw angles show apparent differences in terms of the magnitude, the locations of the $\Delta k$ peaks and the width of the zone with increased turbulence. Figure 14(e–h)show the same results with the abscissa scaled by $U^{2}_T$ and the ordinate shifted with $Y_C(x)$ then scaled by $R_N$. With this scaling, the $\Delta k$ profiles of different yaw angles overlap with each other in the far wake, which shows that the characteristic velocity $U_T$ is still the proper velocity scale for the wake of yawed wind turbines. By further probing into the $\Delta k$ profiles in the near wake, it is found that the profile of $\Delta k$ contains two peaks with the distance between the two peaks scaled well by $R_N$, as the location of $\Delta k$ peaks appear at $Y_C \pm S$, which is consistent with results in the literature (Schottler et al. Reference Schottler, Bartl, Mühle, Sætran, Peinke and Hölling2018). It is also observed that $\Delta k$ in the near wake ($x = 3D$, figure 14a,e) is asymmetric when $\gamma \neq 0^{\circ }$ (especially for $\gamma = 30^{\circ }$) and $\Delta k$ behind the trailing half of the rotor ($y-y_C<0$) is larger than that behind the leading half ($y-y_C>0$), which is similar to the flow passed an inclined circular disc (Calvert Reference Calvert1967; Gao et al. Reference Gao, Tao, Tian and Yang2018). This asymmetry increases with $\gamma$ and cannot be captured by the proposed velocity scale $U_T$. However, it only manifests in the near wake and may be considered as immaterial in utility-scale wind farms, where the turbine spacing is often larger than $5D$ (Méchali et al. Reference Méchali, Barthelmie, Frandsen, Jensen and Réthoré2006). At far-wake locations, this asymmetry disappears with the $\Delta k/U_T^{2}$ profiles overlapping well with each other on the $(y-Y_C)/R_N$ coordinate.

Figure 14. The turbine-added turbulence kinematic energy $\Delta k$ at downstream locations of ($x \in \{3D,5D,7D,9D\}$) on the hub-height plane ($z = z_{hub}$) for different yaw angles and $\tilde {\lambda }$ : (a–d) normalized by the square of the inflow velocity $U_\infty ^{2}$; (e–h) the wake centres are shifted with respect to $Y_C$ and normalized by the length scale $R_N$, with the values normalized by the square of the velocity scale $U^{2}_T$.

Figure 15 plots the turbine-added Reynolds shear stress $\Delta \langle u^{\prime } v^{\prime } \rangle$ at the hub-height plane ($z = z_{hub}$) presented in the same manner as figure 14. In figure 15(a), it is observed that $\Delta \langle u^{\prime } v^{\prime } \rangle$ profiles are featured by complex variations near the wake centre at $x=3D$, which is caused by the turbine nacelle and root loss near the hub. When $\Delta \langle u^{\prime } v^{\prime } \rangle$ are plotted on the coordinate $(y-Y_C)/R_N$, the turbine-added Reynolds shear stress profiles normalized using $U_T^{2}$ overlap with each other in the far wake $(x \ge 5D)$ for different yaw angles and $\tilde {\lambda }$, which confirms $U_T$ as the proper velocity scale for the turbine-added Reynolds stresses.

Figure 15. The turbine-added Reynolds shear stress $\Delta \langle u'v' \rangle$ on the hub-height plane ($z = z_{hub}$): (a–d) normalized by the square of the inflow velocity $U_\infty ^{2}$; (e–h) the wake centres are shifted with respect to $Y_C$ and normalized by the length scale $R_N$, with the values normalized by the square of the velocity scale $U^{2}_T$. For the legend, see figure 14.

4.2. Statistics of instantaneous wake centre and width

This subsection analyses the influence of the yaw misalignment on the dynamics of the horizontal large-scale wake motion (wake meandering), which has a significant impact on the wake expansion, recovery and the fatigue loads on downstream wind turbines (Ainslie Reference Ainslie1988; Högström et al. Reference Högström, Asimakopoulos, Kambezidis, Helmis and Smedman1988; Larsen et al. Reference Larsen2007). In the same way as for the time-averaged wake centre, the instantaneous wake centre is defined as the centre of the Gaussian fit of the instantaneous velocity deficit, which is obtained by subtracting the velocity field computed from the case without wind turbines at exactly the same instant. Figure 16 shows the velocity deficit behind a turbine with $\gamma = 0^{\circ }$ and a turbine with $\gamma = 30^{\circ }$ at the same instant. To obtain the wake centre position and the width along the downstream locations, the velocity deficit $\Delta U$ profile is first spatially filtered with a width of $0.5D$ in the streamwise direction ($x$) and then is fitted to the Gaussian function at each downstream location with (3.3). The curve_fit function in the scipy.optimize package is employed for the fitting process. At each time step, the fitting starts from the near wake with the fitting results serving as the initial guess for the next downstream location to help to converge the fitting process. The filtered velocity deficit and the fitted Gaussian curve are illustrated in figure 16(b) for $\gamma = 0^{\circ }$ and figure 16(d) for $\gamma = 30^{\circ }$ both with $\tilde {\lambda } = 8$. Despite some fluctuations in the instantaneous velocity profiles (blue lines), the fitted Gaussian curves (red dashed lines) capture the wake centre position $y_c$ and width $r_{1/2}$ satisfactorily. In figure 16(a,c), the fitted wake centre lines $y_c(x)$ and the characteristic wake half-width (defined by $r_{1/2}(x)$) are plotted with red dotted lines and red dashed lines, respectively. Generally, they provide a good estimation of the overall trend of the wake in the entire region except for $9D< x<10D$, where the velocity deficit is small and irregular. It is found that the wake behind a yawed wind turbine is generally narrower and that the velocity deficit is smaller than the wake behind a non-yawed turbine. This observation is in accordance with the time-averaged wake quantities. In the following, we examine whether the proposed scales for the time-averaged wake are still appropriate for the statistics of the instantaneous wakes. The following analyses in this section are based on the data collected on the horizontal plane located at the turbine hub height with inflow 1.

Figure 16. Instantaneous flow field behind the wind turbines at hub-height plane $z = z_{hub}$. The contour of the velocity deficit $\Delta U$ at the same simulation time for yaw angles $\gamma = 0^{\circ }$ (a) and $\gamma = 30^{\circ }$ (c). The black solid lines illustrate the wind turbines. Red dash–dotted lines denote the wake centreline $y_c(x)$ and the red dashed lines denote the wake width defined by $r_{1/2}(x)$ obtained from the Gaussian fit of the instantaneous streamwise velocity. The velocity deficit at $x = 7D$ and the corresponding Gaussian fit for yaw angles $\gamma = 0^{\circ }$ (b) and $\gamma = 30^{\circ }$ (d).

We first show the probability density function (p.d.f.) of the wake centre location at $x \in \{5D,7D,9D\}$ for different yaw angles in figure 17. In this figure, panels in the same rows are at the same downstream location and the panels in the same columns are at the same yaw angle. For brevity, only the results of $\tilde {\lambda } = 8$ are plotted as the results of all considered $\tilde {\lambda }$ were similar. The p.d.f. is plotted as grey bars and is fitted to normal distribution curves. In each panel, two vertical lines are plotted to denote the wake centreline location obtained from the time-averaged velocity field ($Y_C$) and the mean value of the instantaneous centre location ($\overline {y_c}$), respectively. These two lines overlap. When increasing the yaw angle, the p.d.f. moves in the $-y$ direction owing to the wake deflection. In the last column, the p.d.f. profiles for the different yaw angles are shifted with respect to $Y_C$ and plotted in the same figure, and are shown to overlap with each other. This indicates that the transverse distribution of the instantaneous wake centres is independent of the yaw angle at different downstream locations for the present cases, which further indicates that the amplitude of wake meandering, defined as the standard deviation of the wake centre location, does not change with yaw angle. By comparing the p.d.f.s at different locations, it is found that the wake centres distribute in a wider range as the wake travels downstream, which is in accordance with the wake expansion phenomenon. It is also noted that the p.d.f. distribution is not exactly Gaussian. To examine this, we compute the skewness of the distribution in table 3. The skewness for cases with non-zero yaw angles generally has a positive value, which shows that the distribution has a longer tail on the right (Lovric Reference Lovric2011), i.e. in the opposite direction of the mean wake deflection. However, because these skewness is rather small (in the range of $-0.5$ to 0.5), the distribution is thus approximately symmetrical and can be approximated by the normal distribution (Cramer Reference Cramer2002). Similarly, the p.d.f.s of the instantaneous wake width $r_{1/2}$ are compared in figure 18. The peak of the p.d.f. is found to move towards the left when increasing the yaw angle, which indicates that the wake width decreases with yaw angle. Travelling downstream, the peak moves slightly to the right owing to the wake expansion. For the cases with yaw angles $\gamma = 0^{\circ }$ and $10^{\circ }$, the distributions of $r_{1/2}$ are well represented by the normal distribution, while for the cases with $\gamma =20^{\circ }$ and $30^{\circ }$, the p.d.f.s are slightly skewed with a longer right tail especially at $x = 9D$. On the other hand, the standard deviation of $r_{1/2}$ is not greatly affected by the yaw angles, as observed in the last column.

Figure 17. The probability density function of the instantaneous wake centre at $x \in \{5D,7D,9D\}$ for different yaw angles at $\tilde {\lambda } = 8$. The histogram is fitted to a Gaussian distribution function in each panel, and the fitted curves are compared in panels (e,j,o). The orange dashed vertical lines denote the mean value of the instantaneous wake centre locations and the red dotted vertical lines are plotted at the wake centre obtained from the time-averaged velocity shown in figure 6.

Figure 18. The probability density function of the instantaneous wake width at $x \in \{5D,7D,9D\}$ for different yaw angles at TSR $\tilde {\lambda } = 8$. The histogram is fitted to a Gaussian distribution in each panel, and the fitted curves are compared in panels (e,j,o).

Table 3. The skewness of the instantaneous wake centre position $y_c$ for different yaw angles. The skewness is defined as $g_1 = {m_3}/{\sigma _{y_c}^{3}}$ with $m_3 = ({1}/{N}) \sum _{i=1}^{N} (y_c -\overline {y_c})^{3}$ and $N$ is the number of samples.

We further analyse the streamwise evolution of instantaneous wake characteristics at different $\gamma$ and $\tilde {\lambda }$. In figure 19(a), the streamwise variations of $\overline {y_c}$ computed from different yaw angles and $\tilde {\lambda }$ are compared. The mean of the instantaneous wake centre locations is shown to scale well using the proposed length scale $Y_N$ and agrees with the wake centre obtained from the time-averaged field $Y_C$ (in grey). Figure 19(b) plots the streamwise evolution of the standard variation of the wake centre position $\sigma _{y_c}$, which indicates the amplitude of the meandering motion of the wake centre around its mean position. In figure 19(b), the rotor diameter $D$, instead of $Y_N$, is used for the normalization. Here, the $\sigma _{y_c}$ for the different cases overlap with each other. This observation further confirms the finding from the p.d.f. of the wake centre in figure 17, where the wake meandering amplitudes are independent of yaw angle and TSR. Furthermore, it is observed that $\sigma _{y_c}$ increases approximately linearly with the distance to the turbine, which shows that the meandering amplitude increases as the wakes travel downstream, and this conforms with Taylor's frozen eddy hypothesis.

Figure 19. Streamwise distribution of the statistics of instantaneous quantities of the wakes behind yawed wind turbines, where the values are scaled with the same factors as for the time-averaged quantities: (a) the mean value of instantaneous wake centre positions $y_c(x)$, scaled by $Y_N$; (b) the standard deviation of $y_c(x)$ normalized by the rotor diameter $D$; (c) the mean value of the characteristic instantaneous wake width $r_{1/2}(x)$ scaled by $R_N$; (d) the standard deviation of $r_{1/2}(x)$ normalized by $D$. In (a,c), the grey curves plot the corresponding normalized characteristics obtained from the time-averaged velocity field.

The instantaneous wake width $r_{1/2}$ is analysed in figure 19(c,d). From figure 19(c), it is found that the mean value of the instantaneous wake width $\overline {r_{1/2}(x)}$ for different $\gamma$ and $\tilde {\lambda }$ can be properly scaled by the proposed length scale $R_N$. When compared with the time-averaged wake width $R_{1/2}(x)$ shown in grey, the widths defined by both approaches agree well in the near wake ($x=1D$), then $R_{1/2}$ becomes larger at further downstream locations. This is because $R_{1/2}$, defined using the time-averaged velocity field, combines both the expansion of the instantaneous wake width $r_{1/2}$ and the wake meandering, and the latter smears the instantaneous velocity deficit over a wider region; in contrast, the wake width follows the instantaneous wake centre $y_c$ and expands slower. The difference between $r_{1/2}$ and $R_{1/2}$ increases with turbine downstream distance in accordance with the meandering amplitude represented by $\sigma _{y_c}$. Figure 19(d) plots the standard deviation of the instantaneous wake half-width $\sigma _{r_{1/2}}$, which represents the wake deformation that expands and shrinks with time. It is noticed that the standard deviations of the wake width for different yaw angles and $\tilde {\lambda }$ overlap with each other when normalized by the rotor diameter $D$, which suggests again that the wake deformation is not greatly influence by the turbine operation condition. Furthermore, the standard deviation of $r_{1/2}$ gradually increases travelling in the streamwise direction, which indicates a larger wake deformation in the far wake.

In summary, figure 19 mainly reveals that the mean of the instantaneous wake quantities and their standard deviation are scaled differently. The proposed scaling factors derived for the time-averaged velocity field still work for the mean value of $y_c$, $r_{1/2}$, but the standard deviations of these quantities, on the other hand, are independent of the length scales defined using the turbine operation condition. The influence of inflow turbulence on the meandering amplitude and the wake deformation will be analysed in detail in § 5.

5.1. Characteristics of the time-averaged wake

Figure 20 shows the characteristics of the time-averaged wake for different yaw angles and different inflow conditions. In figure 20(a), the wake centreline deflection is compared. When comparing curves from the cases with the same inflow, the wake deflections $Y_C(x)$ for the cases with different yaw angles overlapped with each other when normalized by $Y_N$. Furthermore, a reduction of the wake deflection was observed at far-wake locations when the inflow turbulence intensity was increased, which is in agreement with previous studies (Jiménez et al. Reference Jiménez, Crespo and Migoya2010; Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2016; Bartl et al. Reference Bartl, Mühle, Schottler, Sætran, Peinke, Adaramola and Hölling2018). Furthermore, the spanwise velocity $V_C(x)$ in figure 20(b) was found to decrease with increasing inflow turbulence intensity and the normalized values in general agreed with each other for the same inflow. The influence of the inflow turbulence on the streamwise velocity deficit is shown in figure 20(c). The normalized streamwise velocity deficits $\Delta U_C(x)/U_N$ showed a remarkably good agreement between different yaw angles, but they decreased with the turbulence intensity. Moreover, all the curves started at similar $y$ intercepts at approximately 1 in the near wake, which showed that the effect of the inflow turbulence on the imminent wake is minor, where the streamwise velocity deficit can be predicted well by the one-dimensional momentum theory. In contrast, the effect of the inflow condition becomes important in the far wake, where it affects wake recovery, as reported in the literature (Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2017). The length of the near-wake region can be defined using the maximum streamwise velocity deficit. The so-computed near-wake region is shown in figure 20(d). Interestingly, it is found that increasing the yaw angle and inflow turbulent intensity both result in a shorter near-wake region, with a larger impact from the inflow turbulence. The streamwise variation of the wake width $R_{1/2}$ is shown in figure 20(e). The curves overlap well with each other in the near wake regardless of the inflows. Beyond this near-wake region, the wake expands approximately linearly. For the cases with higher turbulence intensity, the wake width is also larger as the wake expansion starts earlier and grows faster. It is also worth noting that small discrepancies of the normalized wake width for different yaw angles are observed in the far wake at $x>6D$, especially for the inflow 3 cases with the highest turbulence intensity. This indicates that the length scale $R_N$, which is derived from the one-dimensional momentum theory, may encounter difficulties to scale the far-wake width when the inflow turbulence intensity is high.

Figure 20. Effects of inflow turbulence on the statistics of time-averaged wakes: (a) the centreline deflection; (b) the characteristic spanwise velocity; (c) the characteristic streamwise velocity deficit; (d) near-wake length defined using the streamwise location where the velocity deficit reaches the maximum; and (e) the width of the wake.

The transverse profiles of the streamwise velocity deficit $\Delta U$ and the transverse velocity $V$ are presented in figure 21. In general, these results confirm that the proposed scaling factors are able to overlap the profiles of different yaw angles for the same inflow, and that the magnitudes of the turbine induced streamwise velocity deficit $\Delta U$ and the transverse velocity $V$ decrease faster for inflows with higher turbulence intensity.

Figure 21. Effects of inflow turbulence on horizontal profiles of the normalized time-averaged wake velocity: (a–d) the streamwise velocity deficit; and (e–h) the spanwise velocity. See figure 20 for the legend.

5.2. Characteristics of instantaneous wakes and the turbine-added turbulence

The statistics of the turbine-added turbulent characteristics, which include the added TKE $\Delta k$, and the added Reynolds shear stress $\Delta \langle u'v' \rangle$ are plotted in figure 22. These curves are normalized by the proposed velocity scale $U_T$ and $R_N$, as in § 4. In general, these curves have similar shapes regardless of the incoming flow. Under the same inflow condition, the curves overlap, which indicates that the proposed velocity and length scales can still describe the similarity of these characteristics for different yaw angles. Because these scales do not consider the influence of the inflow condition, discrepancies between the curves from cases with different inflows are observed, as expected. For $\Delta k$, the influence of inflow is the most obvious at $x=7D$, where the added TKE decreases when the inflow is changed from inflow 1 and inflow 2 to inflow 3, which is in agreement with the work of Crespo & Herna (Reference Crespo and Herna1996). For the added Reynolds stress $\Delta \langle u'v' \rangle$, an overall good scaling is observed for the cases under different inflows.

Figure 22. Effects of inflow turbulence on the normalized added turbulence: (a–d) the added kinetic turbulent energy; and (e–h) the added in-plane Reynolds shear stress.

The effects of inflow conditions on the instantaneous wake quantities are examined in figure 23. The same technique explained in § 4 is employed to extract the instantaneous wake quantities from the flow field for the three inflow conditions. In figure 23(a), the mean of the instantaneous wake centre location $\overline {y_c}/Y_N$ as a function of downstream location is compared with the wake centreline $Y_C/Y_N$ obtained from the time-averaged velocity field. The mean wake centrelines defined by these two approaches agree fairly well with each other even in the far wake. Figure 23(b) shows the standard deviations of the instantaneous wake position for different inflows, as in § 4. It is seen that $\sigma _{y_c}(x)$ is independent of yaw angle for the considered cases. Considering that $\sigma _{y_c}(x)$ essentially represents the meandering amplitude, the observation in figure 23(b) suggests that the wake meandering is dominated by the incoming large eddies instead of turbine operation for the present cases, which is consistent with the works of Bastankhah & Porté-Agel (Reference Bastankhah and Porté-Agel2017) and Yang & Sotiropoulos (Reference Yang and Sotiropoulos2019b). Figure 23(c) compares the mean of the instantaneous wake width $\overline {r_{1/2}}$ with the wake width obtained from the time-averaged velocity field $R_{1/2}$ (plotted in grey). The normalized instantaneous wake widths $\overline {r_{1/2}}/R_N$ overlap well with each other for all the considered inflow conditions and yaw angles, which shows that the near-wake width $R_N$ is a good scaling factor even in the far wake. Moreover, the mean of the instantaneous wake width grows very slowly at all the considered locations $(1D< x<10D)$, which suggests that the expansion of the instantaneous wake width is much weaker than that of the time-averaged wake. This should be taken into account in wake meandering models (e.g. the dynamic wake meandering model Larsen et al. Reference Larsen, Madsen, Thomsen and Larsen2008). Furthermore, it is observed that the instantaneous wake widths are not very sensitive to the inflow condition, i.e. the effects of the inflow turbulence on the expansion of the instantaneous wake are negligible for the considered cases. This further suggests that the faster recovery of the time-averaged wake for cases with higher inflow turbulence intensity is mainly caused by stronger meandering motions. Lastly, the effects of inflow on the standard deviation of the instantaneous wake width are shown in figure 23(d). Here, $\sigma _{r_{1/2}}$ is shown to have a similar behaviour to $\sigma _{y_c}$, i.e. it increases linearly with downstream distance, and the slope of the increase is nearly independent of turbine yaw angles but increases with the turbulence intensity of the inflow, which indicates stronger instantaneous wake deformation for higher inflow turbulence intensities. Moreover, this figure clearly reveals that increasing the inflow turbulence intensity (TI) leads to an increase in the growth rates of both $\sigma _{y_c}$ and $\sigma _{r_{1/2}}$. The empirical formula for the dependence of the expansion rate of the time-averaged wake width on the inflow turbulence intensity can be found in the literature e.g. $k = 0.38 \textrm {TI} + 0.004$ proposed by Niayifar & Porté-Agel (Reference Niayifar and Porté-Agel2016). However, to the best of the authors’ knowledge, similar empirical formulae for $\sigma _{y_c}$ and $\sigma _{r_{1/2}}$ do not exist yet, which need further investigation.

Figure 23. Effects of inflow turbulence on the instantaneous wake quantities: (a) the mean of the instantaneous wake centre positions $y_c$; (b) the standard deviation of the instantaneous wake centre position $\sigma _{y_c}$; (c) the mean of the characteristic instantaneous wake width $r_{1/2}$; and (d) the standard deviation of the characteristic instantaneous wake width $\sigma _{r_{1/2}}$. In (a,c), the curves in grey scale plot the corresponding terms obtained from the time-averaged wake. In (b,d), the growth rates $k$ of $\sigma _{y_c}$ and $\sigma _{r_{1/2}}$ are illustrated by a short dashed line for each inflow condition.

5.3. Relation between the time-averaged wake and the instantaneous wake

This section is dedicated to relating the time-averaged wake and the instantaneous wake to the standard deviation of the wake centre position $\sigma _{y_c}$, which represents the wake meandering amplitude.

To derive the relations, hypotheses are made as follows. At any downstream location $x$,

1. (i) the wake centre position $y_c$ follows the normal distribution $\mathcal {N}(\overline {y_c},\sigma _{y_c}^{2})$, i.e. the probability function of the wake centre is

   (5.1)\begin{equation} \mathcal{P} (y_c=\xi) = \frac{1}{\sqrt{2{\rm \pi}}\sigma_{y_c}}\exp{\left(-\frac{(\xi-\overline{y_c})^{2}}{2\sigma_{y_c}^{2}}\right)}, \end{equation}
   as shown in figure 17;

2. (ii) the transverse profiles of instantaneous velocity deficits can be approximated by a Gaussian function with $\Delta u (y,\xi ) = \Delta u_c \exp {(-({(y-\xi )^{2}}/{2 s^{2}}))}$ as shown in figure 16, where $s$ is the standard deviation of the Gaussian function related to the wake width with $r_{1/2} = \sqrt {2\ln {2}} s$ and $\xi$ is the instantaneous wake centre.

The time-averaged velocity deficit $\Delta U$ can be regarded as the ensemble-average of the instantaneous velocity deficits, so it can be calculated by multiplying the p.d.f. of the wake centre $P(y_c = \xi )$ with the instantaneous wake deficit function $\Delta u(y,\xi )$ and integrating $\xi$ from $-\infty$ to $+\infty$, as follows,

(5.2)\begin{align} \Delta U (y) & = \int_{-\infty}^{+\infty} \mathcal{P}(y_c = \xi) \Delta u(y,\xi) \,\textrm{d} \xi\nonumber\\ & = \int_{-\infty}^{+\infty} \frac{1}{\sqrt{2{\rm \pi}}\sigma_{y_c}}\exp{\left(-\frac{(\xi-\overline{y_c})^{2}}{2\sigma_{y_c}^{2}}\right)} \Delta u_c \exp{\left(-\frac{(y-\xi)^{2}}{2 s^{2}}\right)}\, \textrm{d} \xi\nonumber\\ & = \sqrt{\frac{s^{2}}{\sigma_{y_c}^{2} + s^{2}}} \Delta u_c \exp{\left(\frac{(y-\overline{y_c})^{2}}{2(\sigma_{y_c}^{2} + s^{2})}\right)}. \end{align}

Comparing (5.2) and (3.3), one obtains the following relations,

(5.3)$$\begin{gather} \Delta U_C = \sqrt{\frac{s^{2}}{\sigma_{y_c}^{2} + s^{2}}} \Delta u_c , \end{gather}$$

(5.4)$$\begin{gather}S = \sqrt{\sigma_{y_c}^{2} + s^{2}}. \end{gather}$$

Recall the definition of the wake width is $R_{1/2} = \sqrt {2 \ln 2} S$ and $\overline {r_{1/2}} = \sqrt {2 \ln 2} s$ for the time-averaged and instantaneous wakes, respectively. One obtains the relation of $R_{1/2}$, $\overline {r_{1/2}}$ and $\sigma _{y_c}$ as

(5.5)\begin{equation} R_{1/2}(x) = \sqrt{2\ln{2}~ \sigma^{2}_{y_c}(x)+\overline{r_{1/2}(x)}^{2} }. \end{equation}

This relation is validated with the simulation results for all the inflow conditions and yaw angles in figures 24 and 25. Figure 24 compares the wake width predicted by (5.5) with that from simulation. The panels (a,d,g) show clearly the difference between the instantaneous wake width (in colour) and the time-averaged wake width (in grey), whereas the panels (c,f,i) show an excellent agreement between the time-averaged wake width calculated by (5.5) (in colour) and the simulation results (in grey). Figure 25 compares the mean of the instantaneous centreline velocity deficit $\Delta u_c$ with the time-averaged centreline velocity deficit $\Delta U_C$. The panels (a,d,g) show that the discrepancies between these two definitions increase with downstream distance owing to the increasing wake meandering amplitudes in the far wake, as shown in panels (b,e,h). In panels (c,f,i), the prediction $\tilde {\Delta U}_C$ using (5.3) agrees well with the computed $\Delta U_C$ from the time-averaged simulation results. This agreement shows that the contribution of the wake meandering can be well captured by (5.3).

Figure 24. $(a{,}d{,}g)$ Mean of the instantaneous wake width $\overline {r_{1/2}}$ (in colour) compared with the time-averaged wake width $R_{1/2}$ (in grey). $(b{,}e{,}h)$ Standard deviation of the wake centre position. $(c{,}\,f{,}i)$ The predicted time-averaged wake width $\tilde {R}_{1/2}(x) = \sqrt {2\ln {2} \sigma _{y_c}^{2}(x)+\overline {r_{1/2}(x)}^{2} }$ (in colour) compared with the time-averaged simulation results ${R}_{1/2}$ (in grey).

Figure 25. $(a{,}d{,}g)$ Mean of the instantaneous velocity deficit $\overline {\Delta u_c}$ (in colour) compared with the time-averaged wake velocity deficit $\Delta U_C$ (in grey). $(b{,}e{,}h)$ Standard deviation of the wake centre position. $(c{,}\,f{,}i)$ The predicted time-averaged velocity deficit $\tilde {\Delta U}_C(x)= \sqrt {{s^{2}}/({s^{2} + \sigma _{y_c}^{2}})}$ (in colour) compared with the simulation results ${U}_{C}$ (in grey). See figure 24 for the legend.

Physically speaking, the derived relations and observations in figures 24 and 25 show the relative importance of (i) wake meandering and (ii) turbulence diffusion and pressure (on the moving frame of reference located at the instantaneous wake centre) on the recovery of the wake, i.e. the decrease of velocity deficit and the growth of wake radius. When the wake meandering amplitude is small ($\sigma _{y_c} < r_{1/2}$) in the near wake, the recovery of the time-averaged wake is dominated by the second factor mentioned above. Otherwise, when the meandering amplitude is large ($\sigma _{y_c} > r_{1/2}$), the recovery of the time-averaged wake is mainly caused by the wake meandering. Additionally, because $\sigma _{y_c}$ is independent of yaw angle, this equation explains why the far-wake width, where the inflow turbulence is dominant, is poorly scaled with $R_N$ (derived based on the turbine operation condition for the near wake). These findings and the derived relations can be employed in the development of advanced wake meandering models to take into account more accurately the instantaneous wake features instead of using models based on the time-averaged wake (Jiménez et al. Reference Jiménez, Crespo and Migoya2010; Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2016; Shapiro et al. Reference Shapiro, Gayme and Meneveau2018).