2.1. Evolving vortex sheet governing equations

As shown in Shapiro et al. (Reference Shapiro, Gayme and Meneveau2020), among others, vortices shedding from the circumference of the yawed rotor represented as an actuator disk form a tubular vortex sheet. The objective of this section is to model the shape evolution of this vortex sheet with downstream distance or, equivalently, with time. Only the streamwise component of the shedding vorticity is modelled in this work because the lateral wake deflection and the deformation of the wake cross-section are mainly due to the velocities induced by the streamwise component of shedding vorticity (Martinez-Tossas & Branlard Reference Martinez-Tossas and Branlard2020). The vortex sheet consists of semi-infinite streamwise vortex lines. In order to enable solving the governing equations analytically, following Coelho & Hunt (Reference Coelho and Hunt1989) and Rottman et al. (Reference Rottman, Simpson and Stansby1987), we assume that the vortex sheet is planar and that its constituent vortex lines are infinite instead of semi-infinite, an approximation that improves at increasing distances from the origin.

Figure 1 shows a schematic of the vortex sheet in the plane normal to the incoming flow. The $(x,y,z)$ coordinate system is defined with an origin at the rotor centre and with $x$ in the streamwise direction (i.e. parallel to the incoming flow) and $y$ and $z$ in the spanwise and vertical directions, respectively. Alongside this Cartesian system, we define a polar coordinate system of $(r,\theta )$ in the $yz$-plane (i.e. plane normal to the incoming flow). This polar coordinate system is attached to the centre of the vortex sheet, denoted by $C$. The position of $C$ in the $(x,y,z)$ coordinate system is denoted by $(x_c,y_c,z_c)$. As $r$ is the radial distance from $C$, one can write $r^{2}=(y-y_c)^{2}+(z-z_c)^{2}$. The polar angle $\theta$ is measured from the positive $y$-axis toward the positive $z$-axis such that $\tan \theta =(z-z_c)/(y-y_c)$. The shape of the vortex sheet is represented by the polar function $\xi (\theta,t)$ that measures the distance of the vortex sheet from the centre, where $t$ is time. Our main objective is to describe the evolution of the vortex sheet as a function of time $t$ in a frame moving downstream with the convection velocity $U_{con}$, i.e. determine $\xi (\theta,t)$. This is equivalent to determining the downstream spatial evolution of the vortex sheet, with $x=U_{con} t$. The sheet location $\xi (\theta,t)$ obeys

(2.1)\begin{equation} \xi=\xi_0 + \int u_r(\theta,t)\,\textrm{d}t, \end{equation}

where $\xi _0=\xi (\theta,0)$ and $u_r(\theta,t)$ denotes the radial velocity of the vortex sheet, which is affected by the strength of the evolving vortex sheet, whose evolution is treated next.

Figure 1. Schematic of the vortex sheet and different velocity terms on the right-hand side of (2.4). (a) Self-induced vortex sheet velocity, $\boldsymbol {u}_I$. (b) Vortex sheet velocity induced by the point vortex at the vortex sheet centre, $\boldsymbol {u}_{I\!I}$. (c) Velocity of the vortex sheet centre, $\boldsymbol {u}_{c}$.

Let us denote the strength of the vortex sheet by $\gamma =\gamma (\theta,t)$, where the vortex strength $\gamma$ is defined as the amount of circulation per unit length. In addition to the vortex sheet, there is a point vortex at $C$ with a circulation of $\varGamma _r$ (see figure 1b) to model the rotor root vortex, which is elaborated later in § 2.3. We will show in § 2.3 that the initial condition for the strength can be written as

(2.2)$$\begin{gather} \gamma_0=\gamma(\theta,0)=\gamma_r+\gamma_b\sin\theta \quad \textrm{at}\ r=\xi_0, \end{gather}$$

(2.3)$$\begin{gather}\varGamma_r={-}2{\rm \pi}\xi\gamma_r \quad\textrm{at } r=0, \end{gather}$$

where $\gamma _r$ and $\gamma _b$ are constants depending on turbine operating conditions. Later in § 2.3, we show that $\gamma _b$ is related to the vorticity generated due to turbine yaw, while $\gamma _r$ and $\varGamma _r$ are related to the vorticity generated by turbine rotating blades. Our focus now is to predict the deformation of the vortex sheet provided that the initial conditions are given by (2.2)–(2.3).

The velocity of the vortex sheet $\boldsymbol{u}(\theta,t)$ with respect to a coordinate system attached to the centre $C$ is given by

(2.4)\begin{equation} \underbrace{\boldsymbol{u}(\theta,t)}_{\substack{\text{vortex sheet velocity} \\ \text{with respect to C}}}=\underbrace{\boldsymbol{u}_{I}}_{\substack{\text{self-induced} \\ \text{vortex sheet velocity}}}+\underbrace{\boldsymbol{u}_{II}}_{\substack{\text{vortex sheet velocity induced} \\ \text{by point vortex at C}}}-\underbrace{\boldsymbol{u}_{c}}_{\substack{\text{velocity} \\ \text{of C}}}. \end{equation}

In (2.4) and hereafter, bold letters denote vectors. Next, we employ the Biot–Savart law to determine the three velocity terms on the right-hand side of (2.4), starting with the self-induced velocity $\boldsymbol {u}_{I}=u_{{I},r}\boldsymbol {e}_r+u_{{I},\theta }\boldsymbol {e}_{\theta }$, where $\boldsymbol {e}_r$ and $\boldsymbol {e}_{\theta }$ are unit vectors in the radial and tangential directions, respectively. The radial $u_{{I},r}$ and tangential $u_{{I},\theta }$ components of the self-induced velocity at a given polar angle of $\theta$ are respectively given by

(2.5)$$\begin{gather} u_{I,r}(\theta,t)=\int_0^{2{\rm \pi}}\frac{\gamma(\theta',t)\sin{\alpha}}{2{\rm \pi} l}\xi'\,\textrm{d}\theta', \end{gather}$$

(2.6)$$\begin{gather}u_{I,\theta}(\theta,t)=\int_0^{2{\rm \pi}}\frac{\gamma(\theta',t)\cos{\alpha}}{2{\rm \pi} l}\xi'\,\textrm{d}\theta', \end{gather}$$

where $l$ and $\alpha$ are defined in figure 1(a), $\theta '$ is a dummy integration variable, and $\xi '=\xi (\theta ',t)$. According to the law of cosines, $l^{2}=\xi ^{2}+\xi ^{'2}-2\xi \xi '\cos {(\theta '-\theta )}.$ The angle $\alpha$ shown in the figure 1(a) can be related to $l$ based on the law of sines for the drawn triangle, which results in $\sin {\alpha }=(\xi '/l)\sin {(\theta '-\theta )}.$ For small values of time $t$ and yaw angle $\beta$, we can assume that the vortex sheet is approximately circular and, thus, $\xi \approx \xi '$. After some trigonometric manipulations, (2.5) and (2.6) can be simplified to

(2.7)$$\begin{gather} u_{I,r} (\theta,t) = {p.v.}\frac{1}{4{\rm \pi}} \int_{0}^{2{\rm \pi}} \frac{\gamma(\theta',t)}{\tan{\left[\left(\theta'-\theta\right)/2\right]}}\textrm{d}\theta', \end{gather}$$

(2.8)$$\begin{gather}u_{{I}, \theta} (\theta,t) = {p.v.}\frac{1}{4{\rm \pi}}\int_{0}^{2{\rm \pi}} \frac{\gamma(\theta',t)\sin{\left[\left(\theta'-\theta\right)/2\right]}}{\sin{\left[\left(\theta'-\theta\right)/2\right]}}\textrm{d}\theta'= \frac{1}{4{\rm \pi}}\int_{0}^{2{\rm \pi}}\gamma(\theta')\,\textrm{d}\theta'. \end{gather}$$

Both integrals in (2.7) and (2.8) have singularities at $\theta '=\theta$ and, thus, we use the Cauchy principal values of these two integrals. While the principal value (p.v.) of the latter can be simply obtained by removing $\sin [(\theta '-\theta )/2]$ from the numerator and the denominator, the p.v. of the former needs to be determined for a given $\gamma (\theta,t)$.

Next, we determine the velocity of the vortex sheet induced by the point vortex at C as shown in figure 1(b). We obtain

(2.9)$$\begin{gather} u_{{II},r} (\theta,t) = 0, \end{gather}$$

(2.10)$$\begin{gather}u_{{II},\theta} (\theta,t) = \frac{\varGamma_r}{2{\rm \pi} \xi}={-}\gamma_r. \end{gather}$$

Finally, we determine $\boldsymbol {u}_c$, which is the velocity of $C$ induced by the vortex sheet, shown in figure 1(c). It can be readily shown that $\boldsymbol {u}_c$ is given by

(2.11)$$\begin{gather} u_{c,r}(\theta,t)=\frac{1}{2{\rm \pi}}\int_0^{2{\rm \pi}}\gamma(\theta',t) \sin\left({\theta'-\theta}\right)\textrm{d}\theta', \end{gather}$$

(2.12)$$\begin{gather}u_{c,\theta}(\theta,t)={-}\frac{1}{2{\rm \pi}}\int_0^{2{\rm \pi}}\gamma(\theta',t) \cos\left({\theta'-\theta}\right)\textrm{d}\theta'. \end{gather}$$

If we neglect streamwise ($x$-direction) straining, vorticity is a conserved quantity, so the vorticity transport equation for the vortex sheet strength $\gamma (\theta,t)$ provides the additional required evolution equation (Moore Reference Moore1978)

(2.13)\begin{equation} \frac{\partial \gamma}{\partial t}+ \frac{\partial\left(\gamma u_s\right)}{\partial s}=0, \end{equation}

where $s$ is the arclength along the vortex sheet. At small values of time $t$ and yaw angle $\beta$, the vortex sheet remains approximately circular, so $u_s$ and $\partial s$ can be respectively replaced with $u_{\theta }$ and $\xi \partial \theta$. Therefore, (2.13) is simplified to

(2.14)\begin{equation} \frac{\partial \gamma}{\partial t}+ \frac{1}{\xi}\frac{\partial (\gamma u_\theta)}{\partial \theta}\approx 0. \end{equation}

Next, we use $\gamma _b$ and $\xi _0$ to non-dimensionalise variables in (2.4) and (2.14). This leads to the following set of dimensionless equations:

(2.15)$$\begin{gather} \hat{u}_{r}(\theta,\hat{t})\approx {p.v.}\frac{1}{4{\rm \pi}} \int_{0}^{2{\rm \pi}} \frac{\hat{\gamma}(\theta',\hat{t})}{\tan{\left[\left(\theta'-\theta\right)/2\right]}}\textrm{d}\theta' -\frac{1}{2{\rm \pi}}\int_0^{2{\rm \pi}}\hat{\gamma}(\theta',\hat{t})\sin\left({\theta'-\theta}\right)\textrm{d}\theta', \end{gather}$$

(2.16)$$\begin{gather}\hat{u}_{\theta}(\theta,\hat{t})\approx\frac{1}{4{\rm \pi}} \int_{0}^{2{\rm \pi}}\hat{\gamma}(\theta')\,\textrm{d}\theta'+\frac{1}{2{\rm \pi}} \int_0^{2{\rm \pi}}\hat{\gamma}(\theta',\hat{t})\cos\left({\theta'-\theta}\right)\textrm{d}\theta'-\chi, \end{gather}$$

(2.17)$$\begin{gather}\frac{\partial \hat{\gamma}}{\partial \hat{t}} \approx{-}\frac{1}{\hat{\xi}} \frac{\partial (\hat{\gamma} \hat{u}_\theta)}{\partial \theta}. \end{gather}$$

Here $\hat {t}=t\gamma _b/\xi _0$, $\hat {u}=u/\gamma _b$, $\hat {\gamma }=\gamma /\gamma _b$, $\hat {\xi }=\xi /\xi _0$ and $\chi =\gamma _r/\gamma _b$. Note that the dimensionless time $\hat {t}$ becomes negative for negative values of $\gamma _b$. In the following, we solve (2.15)–(2.17) using the power series method.

2.2. Analytical solution using power series approximation

We write $\hat {\gamma }(\theta,t)$, $\hat {u}_r(\theta,t)$ and $\hat {u}_{\theta }(\theta,t)$ as power series in the form of

(2.18a–c)\begin{equation} \hat{\gamma}(\theta,\hat{t}) = \sum _{n=0}^{\infty} \hat{\gamma}_{n}(\theta)\hat{t}^{n}; \quad \hat{u}_r(\theta,\hat{t}) = \sum _{n=0}^{\infty} \hat{u}_{r n}(\theta)\hat{t}^{n}, \quad \hat{u}_{\theta}(\theta,\hat{t}) =\sum _{n=0}^{\infty} \hat{u}_{\theta n}(\theta)\hat{t}^{n}. \end{equation}

According to (2.1), for $\hat {\xi }$ and the factor ${1}/{\hat {\xi }}$ in (2.17), we have

(2.19)$$\begin{gather} \hat{\xi}=1+\int \hat{u}_r\,\textrm{d}\hat{t}=1+\sum _{n=0}^{\infty} \frac{1}{n+1} \hat{u}_{r n} \hat{t}^{n+1}, \end{gather}$$

(2.20)$$\begin{gather}\frac{1}{\hat{\xi}} = \sum _{n=0}^{\infty} f_{n}(\theta)\hat{t}^{n}, \end{gather}$$

where $f_n$s are Taylor series expansion coefficients of ${1}/{\hat {\xi }}$. For example, the first three coefficients, which are used in the final solution of this paper, can be shown to be $f_{0}=1$, $f_{1}=-\hat {u}_{r 0}$ and $f_{2} =\hat {u}^{2}_{r 0} -\frac {1}{2} \hat {u}_{r 1}$. We insert the power series (2.18a–c) and (2.20) into (2.15)–(2.17) and equating coefficients we obtain

(2.21)$$\begin{gather} \hat{u}_{rn}(\theta)\approx {p.v.}\frac{1}{4{\rm \pi}} \int_{0}^{2{\rm \pi}} \frac{\hat{\gamma}_n(\theta')}{\tan{\left[\left(\theta'-\theta\right)/2\right]}}\textrm{d}\theta' -\frac{1}{2{\rm \pi}}\int_0^{2{\rm \pi}}\hat{\gamma}_n(\theta')\sin\left({\theta'-\theta}\right)\textrm{d}\theta', \end{gather}$$

(2.22)$$\begin{gather}\hat{u}_{\theta n}(\theta)\approx\frac{1}{4{\rm \pi}} \int_{0}^{2{\rm \pi}}\hat{\gamma}_n(\theta')\,\textrm{d}\theta'+\frac{1}{2{\rm \pi}} \int_0^{2{\rm \pi}}\hat{\gamma}_n(\theta')\cos\left({\theta'-\theta}\right)\textrm{d}\theta'- \begin{cases} \chi, & \text{if } n=0, \\ 0, & \text{if } n>0, \end{cases} \end{gather}$$

(2.23)$$\begin{gather}\hat{\gamma}_{n+1} \approx{-}\frac{1}{(n+1)} \left(\sum_{j=0}^{n}f_{j} \sum_{i=0}^{n-j} \frac{\partial (\hat{\gamma}_i \hat{u}_{\theta (n-j-i)})}{\partial \theta}\right). \end{gather}$$

The first term on the right-hand side of (2.21) is a Cauchy p.v. of an improper integral. The following identities are useful to solve this integral:

(2.24)$$\begin{gather} {p.v.} \int_{0}^{2{\rm \pi}} \frac{\sin{n x}}{\tan{\left[\left(x-b\right)/2\right]}}\textrm{d}\kern0.06em x=2{\rm \pi}\cos{nb}, \end{gather}$$

(2.25)$$\begin{gather}{p.v.} \int_{0}^{2{\rm \pi}} \frac{\cos mx}{\tan{\left[\left(x-b\right)/2\right]}}\textrm{d}\kern0.06em x={-}2{\rm \pi}\sin{mb}. \end{gather}$$

Here $n\in \{1,2,\ldots \}$, $m\in \{0,1,2,\ldots \}$ and $b\in [0,2{\rm \pi} ]$. The complete derivation of these integrals can be found in the Appendix A.

From (2.2), $\hat {\gamma }_0(\theta )=\sin \theta +\chi$. One can insert $\hat {\gamma }_0$ into (2.21) and (2.22) to respectively find $\hat {u}_{r0}(\theta )$ and $\hat {u}_{\theta 0}(\theta )$. Values of $\hat {u}_{r0}(\theta )$, $\hat {u}_{\theta 0}(\theta )$ and $f_0$ can be then inserted into (2.23) to find $\hat {\gamma }_1(\theta )$. This recursive process is repeated until reaching the desired order of evaluation for the power series of (2.18a–c). After $\hat {u}_{r}(\theta,\hat {t})$ is obtained using the developed recursive relations, the dimensionless shape of the vortex sheet $\hat {\xi }(\theta,\hat {t})$ is evaluated from (2.19). The solutions for $\hat {\gamma }$, $\hat{u}_r$ and $\hat{u}_{\theta}$ up to ${O}(\hat {t}^{3})$ and for $\hat {\xi }$ up to ${O}(\hat {t}^{4})$ are written below,

(2.26)\begin{align} \hat{\gamma}(\theta,\hat{t}) &= \sin (\theta ) + \chi -\frac{1}{2} \hat{t} \sin (2 \theta ) +\hat{t}^{2} \left(-\frac{1}{4} \chi \cos (2 \theta )+\frac{3}{16} \sin (3 \theta )-\frac{\sin (\theta )}{16}\right) \nonumber\\ &\quad +\hat{t}^{3} \left(\frac{1}{12} \chi ^{2} \sin (2 \theta )-\frac{1}{48} \chi \cos (\theta )+\frac{5}{32} \chi \cos (3 \theta ) \right.\nonumber\\ &\quad \left.+\frac{5}{96} \sin (2\theta )-\frac{7}{96} \sin (4 \theta )\right), \end{align}

(2.27)\begin{align} \hat{u}_r(\theta,\hat{t}) &={-}\frac{1}{4} \hat{t} \cos (2 \theta )+ \hat{t}^{2} \left(\frac{1}{8} \chi \sin (2 \theta )+\frac{3}{32} \cos (3 \theta )\right) \nonumber\\ &\quad + \hat{t}^{3} \left(\frac{1}{24} \chi ^{2} \cos (2 \theta )-\frac{5}{64} \chi \sin (3 \theta )+\frac{5}{192} \cos (2 \theta )-\frac{7}{192} \cos (4 \theta )\right), \end{align}

(2.28)\begin{align} &\qquad\ \hat{u}_{\theta}(\theta,\hat{t}) = \frac{1}{2}\sin (\theta ) - \frac{1 }{2}\chi -\frac{1}{32} \hat{t}^{2} \sin (\theta ) -\frac{1}{96} \hat{t}^{3} \chi \cos (\theta ), \end{align}

(2.29)\begin{align} \hat{\xi}(\theta,\hat{t}) &= 1 -\frac{1}{8} \hat{t}^{2} \cos (2 \theta )+\hat{t}^{3} \left(\frac{1}{24} \chi \sin (2 \theta )+\frac{1}{32} \cos (3 \theta )\right)\nonumber\\ &\quad +\hat{t}^{4} \left(\frac{1}{96} \chi ^{2} \cos (2 \theta)-\frac{5}{256} \chi \sin (3 \theta )+\frac{5}{768} \cos (2 \theta)-\frac{7}{768} \cos (4 \theta )\right). \end{align}

One can compute higher-order terms of (2.26)–(2.29), which may become relevant at increasing values of $\hat {t}$. However, since the above solution is developed based on the assumption that the vortex sheet remains approximately circular, increasing deformation of the vortex sheet makes the solution inaccurate at large values of $\hat {t}$. For practical applications at large $\hat {t}$, in Appendix B we propose an empirical formula that merges smoothly with the theoretical expression at small $\hat {t}$ (i.e. $|\hat {t}|\leq$2), while it has desired reasonable properties at large times (i.e. $|\hat {t}|>$2). In the next section we prove the validity of the initial conditions in (2.2) and (2.3). Moreover, values of $\xi _0$, $\gamma _b$ and $\gamma _r$ are determined as functions of turbine operating conditions.

2.3. Setting vortex sheet initial conditions at yawed turbine location

In this section we determine the vorticity shedding from a yawed rotor disk (i.e. $\gamma _0(\theta )=\gamma (\theta,0)$). According to the Kutta–Joukowsky theorem, lift force is proportional to the amount of circulation around a lifting airfoil. This means that an airfoil can be replaced with a bound vortex. Also, for any airfoil with finite span, free vortices must trail downstream from both sides of the bound vortex to infinity, forming a horseshoe vortex (Anderson Reference Anderson2011). Turbine blades rotate and produce power due to their generated lift force, and vorticity is shed from the root and the tip of rotor blades. In addition, the whole yawed rotor can be assumed as a big finite-span airfoil with the lateral component of the thrust force regarded as the lift force. Therefore, in order to find the total shedding vorticity at the rotor disk, we need to determine those due to both yaw offset and rotating blades. In the following, we assume that the yawed rotor can be modelled as a rotating actuator disk.

2.3.1. Vorticity shedding due to turbine yaw

Prior studies have suggested two different approaches to model the distribution of circulation at a yawed disk. By modelling a yawed disk as a lifting line, the circulation is concentrated on a vertical line at the centre of the rotor with an elliptical distribution spanning from the bottom tip to the top tip of the rotor (Shapiro et al. Reference Shapiro, Gayme and Meneveau2018). Alternatively, vorticity due to yaw offset can be assumed to shed from the circumference of the rotor (Zong & Porté-Agel Reference Zong and Porté-Agel2020). Martinez-Tossas & Branlard (Reference Martinez-Tossas and Branlard2020) used vortex cylinder theory to state the equivalency of these two methods. Shapiro et al. (Reference Shapiro, Gayme and Meneveau2020) proved that both vorticity distributions yield the same induced velocity inside the radius of the rotor disk. To determine the reference circulation density $\gamma _b$ needed in (2.2) and to provide more physical insight, we build upon the literature to show that the equivalency of these two vorticity distributions can be also verified simply by rearranging the position of horseshoe vortices over a yawed disk.

Figure 2(a) shows a schematic of a yawed actuator disk modelled as a lifting surface with a constant vortex strength of $\gamma _b$ in the $\bar {z}$-direction. The coordinate system $(\bar {x},\bar {y},\bar {z})$ is defined based on the rotor plane as shown in figure 2(c), and its respective polar coordinate system ($\bar {r},\bar {\theta }$) is defined such that $\bar {r}=\sqrt {\bar {y}^{2}+\bar {z}^{2}}$ and $\tan \bar {\theta }=\bar {z}/\bar {y}$. The lifting surface shown in figure 2(a) can be envisaged as a surface with an infinite number of horseshoe vortices uniformly distributed across the yawed disk. The bound circulation at a given vertical position $z$ is given by

(2.30)\begin{equation} \varGamma(\bar{z})=\int_{-\sqrt{R^{2}-\bar{z}^{2}}}^{\sqrt{R^{2}-\bar{z}^{2}}} \gamma_b\,\mathrm{d}\bar{y}=2\gamma_b\sqrt{R^{2}-\bar{z}^{2}}, \end{equation}

where $R$ is the rotor radius. Note that according to (2.30), the vertical distribution of the bound circulation for the lifting surface with a constant vortex strength $\gamma _b$ is elliptical. This means that if we concentrate all these horseshoe vortices on a vertical line at the centre of the disk, the lifting surface is transformed to a lifting line with an elliptical distribution of circulation, like the one used by Shapiro et al. (Reference Shapiro, Gayme and Meneveau2018), as shown in figure 2(b). In this case, trailing vortices shed from all along the lifting line, because it consists of horseshoe vortices that vary in size. From (2.30), the maximum value of bound circulation for the lifting line, denoted by $\varGamma _b$, occurs at $\bar {z}=0$, and its value is equal to $2R\gamma _b$. From the lifting line theory, we know that $\varGamma _b=-U_{h} C_T R \cos ^{2}{\beta }\sin {\beta }$ (Shapiro et al. Reference Shapiro, Gayme and Meneveau2018), where $C_T$ is the turbine thrust coefficient. The value of $C_T$ is given by

(2.31)\begin{equation} C_T= \frac{2T}{\rho {\rm \pi}R^{2} U_{h}^{2} \cos^{2}\beta}, \end{equation}

where $T$ is the total magnitude of the turbine thrust force, $\rho$ is the air density and $U_{h}$ is the inflow velocity at the hub height. Note that this definition is the same as the one used in Shapiro et al. (Reference Shapiro, Gayme and Meneveau2018), but it is different from the one used in some other prior studies (e.g.  Burton et al. Reference Burton, Sharpe, Jenkins and Bossanyi1995; Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2016). Since $\gamma _b=\varGamma _b/2R$, the value of $\gamma _b$ is given by

(2.32)\begin{equation} \gamma_b={-}\tfrac{1}{2}U_{h} C_T \cos^{2}{\beta}\sin{\beta}. \end{equation}

Figure 2. Vorticity shedding from a yawed actuator disk, modelled (a) as a lifting surface, and (b) as a lifting line. (c) A schematic of different coordinate systems used in this paper.

Next, we determine the distribution of trailing vortices shedding from the circumference of the lifting surface. The circulation of the vortex shedding from an infinitesimal circumferential element $\mathrm {d}s$, where $\mathrm {d}\bar {s}=R\,\mathrm {d}\bar {\theta }$, is $\mathrm {d}\varGamma _{shed,yaw}=\gamma _b \,\mathrm {d}\bar {y}$. Given that $\mathrm {d}\bar {y}=\mathrm {d}\bar {s}\sin {\bar {\theta }}$, we obtain

(2.33)\begin{equation} \mathrm{d}\varGamma_{shed,yaw}=\gamma_b\sin{\bar{\theta}}\,\mathrm{d}\bar{s}. \end{equation}

For the lifting line, on the other hand, the magnitude of circulation of trailing vortex over the segment $\mathrm {d}\bar {z}$ is equal to -$(\mathrm {d}\varGamma /\mathrm {d}\bar {z})\,\mathrm {d}\bar {z}$ (Anderson Reference Anderson2011). From (2.30) and $\varGamma _b=2R\gamma _b$, we obtain

(2.34)\begin{equation} \textrm{d}\varGamma_{shed,yaw}=\frac{\varGamma_b \bar{z}\,\mathrm{d}\bar{z}}{R\sqrt{R^{2}-\bar{z}^{2}}}. \end{equation}

Using the variable change $\bar {z}=R\sin {\bar {\theta }}$, one can easily show that $\mathrm {d}\varGamma _{shed}$ for the lifting surface at any $\bar {\theta }$ (2.33) is half of that of the lifting line (2.34) at the respective height $\bar {z}$. Note that for the lifting surface, at a given height, trailing vortices shed at both angles of $\bar {\theta }$ and $({\rm \pi} - \bar {\theta })$ with the same magnitude of circulation. Therefore, trailing vortices shedding from the lifting line and the lifting surface vary with height in a similar manner. It is also worth mentioning that the results presented here are in agreement with those obtained from the skewed vortex cylinder theory (Coleman et al. Reference Coleman, Feingold and Stempin1945; Branlard & Gaunaa Reference Branlard and Gaunaa2016). By modelling a yawed turbine wake as a skewed vortex cylinder, Martinez-Tossas & Branlard (Reference Martinez-Tossas and Branlard2020) stated that the dominant vorticity shedding from the rotor is the tangential vorticity vector, which lies in the rotor plane. The streamwise projection of this tangential vorticity is equal to the one found in the present work (2.33) (c.f. (9) in Martinez-Tossas & Branlard Reference Martinez-Tossas and Branlard2020).

2.3.2. Inclusion of wake angular momentum effects

In this section we determine the value of $\gamma _r$ required in (2.2) and (2.3). Based on the method of Joukowsky that models a turbine blade as one single horseshoe vortex with constant bound circulation (see Okulov & Van Kuik (Reference Okulov and Van Kuik2012) for historical background), two free trailing vortices with the same magnitude of circulation are shed from both root and tip ends of each turbine blade. Under the assumption of a large number of blades, this creates a vortex system consisting of a bound vortex disk, an axial root vortex and a tubular vortex sheet as shown in figure 3(a). Let us denote the circulation of the root trailing vortex with $\varGamma _r$. As the amount of circulation along any horseshoe vortice remains constant, the bound circulation on the rotor disk at any radial position should be the same as $\varGamma _r$. According to the Kutta–Joukowsky theorem, the bound circulation over an annular ring at a radial position $\bar {r}$ and thickness of $\textrm {d}\bar {r}$ on the rotor disk generates a lift force $\textrm {d}\boldsymbol {L}$, which amounts to (Okulov & Sørensen Reference Okulov and Sørensen2010)

(2.35)\begin{equation} \textrm{d}\boldsymbol{L}=\rho \boldsymbol{V_0}\times \boldsymbol{e_{\bar{r}}}\varGamma_r\,\textrm{d}\bar{r}, \end{equation}

where $\boldsymbol {V_0}$ is the resultant relative wind velocity experienced by the blade element as shown in figure 3(b). In this figure, $\phi$ denotes the angle between $\boldsymbol {V_0}$ and $\bar {\theta }$-direction, and $U_d$ is the component of $\boldsymbol {V_0}$ in the $\bar {x}$-direction. The tangential component of $\textrm {d}\boldsymbol {L}$ produces power $P$, which is given by

(2.36)\begin{equation} \textrm{d}P=\varOmega\,\textrm{d}Q=\rho\varOmega\varGamma_r U_d \bar{r}\,\textrm{d}\bar{r}, \end{equation}

where $\varOmega$ is the turbine rotational velocity and $\textrm {d}Q$ is the torque generated by the given annular ring. From the axial momentum theory, the power generated by the annular ring can be written as the product of $U_d$ and $\textrm {d}T$, where $\textrm {d}T$ is the thrust force exerted on the annular ring. Therefore, we obtain an additional equation for $\textrm {d}P$ as

(2.37)\begin{equation} \textrm{d}P=U_d\textrm{d}T=\rho U_d U_{\infty}^{2}C_T\cos^{2}\beta({\rm \pi}\bar{r})\,\textrm{d}\bar{r}. \end{equation}

Note that to derive (2.37), the local thrust coefficient for a given annular ring is assumed to be the same as its value for the whole rotor defined in (2.31). This is a correct assumption for the Joukowsky vortex model (Van Kuik, Sørensen & Okulov Reference Van Kuik, Sørensen and Okulov2015). Equating (2.36) and (2.37) leads to

(2.38)\begin{equation} \varGamma_r=\frac{{\rm \pi} R}{\lambda} U_{h}C_T\cos^{2}\beta, \end{equation}

where $\lambda$ is the tip-speed ratio and defined as $\varOmega R/U_{h}$.

Figure 3. (a) Modelling a turbine rotor as a rotating actuator disk. (b) The velocity triangle for a rotor blade element.

As seen in figure 3(a), the trailing vorticity sheds over the circumference of the rotor disk. It is evident that the value of circulation for the vorticity shedding over the circumferential element $\mathrm {d}\bar {s}$ is given by

(2.39)\begin{equation} \textrm{d}\varGamma_{shed,rot}=\gamma_r\,\textrm{d}\bar{s}, \end{equation}

where $\gamma _r=-\varGamma _r/(2{\rm \pi} R)$. The variable $\gamma _r$ denotes the strength of the shedding tubular vortex sheet, and from (2.38),

(2.40)\begin{equation} \gamma_r={-}\frac{1}{2\lambda} U_{h}C_T\cos^{2}\beta. \end{equation}

It is important to note that as discussed by Branlard & Gaunaa (Reference Branlard and Gaunaa2016), the above-mentioned shedding vortices are actually in the direction of the wake centreline axis that forms an angle with the streamwise coordinate $x$. Prior studies (e.g.  Coleman et al. Reference Coleman, Feingold and Stempin1945) however showed that the angle between the wake centreline and the $x$-coordinate is expected to be much smaller than the turbine yaw angle. Therefore, we assume that both axial root and tip shedding vortices are in the $x$-direction for simplicity.

2.3.3. Total vorticity shedding from a yawed turbine

From findings of §§ 2.3.1 and 2.3.2, we can determine the total vorticity shedding from a yawed rotor due to both yaw offset and rotating blades as a function of turbine operating conditions. From (2.33) and (2.39), the value of the dimensionless initial vortex strength $\hat {\gamma }_0$ in the $(r,\theta )$ polar coordinate system is given by

(2.41)\begin{equation} \hat{\gamma}_0(\theta)=\frac{\gamma_0}{\gamma_b}=\sin\theta+\chi \quad \textrm{at}\ r=\xi_0, \end{equation}

where from (2.32) and (2.40),

(2.42)\begin{equation} \chi=\frac{\gamma_r}{\gamma_b}=\frac{1}{\lambda\sin\beta}. \end{equation}

The variable $\chi$, called the rotation rate, is the ratio of the strength of vortex generation due to rotating blades to the one generated due to yaw offset. For the limiting cases of $\lambda =\pm \infty$, $\chi$ is equal to zero, and the shedding vorticity is only due to the yaw offset. As the tip-speed ratio $\lambda$ goes to infinity, the amount of torque generated by the turbine goes to zero. Therefore, according to the conservation of angular momentum, there should be no wake rotation downwind of the actuator disk in the limiting cases of $\lambda =\pm \infty$. Hereafter, the term non-rotating wake refers to the wake of an actuator disk with an infinite tip-speed ratio $\lambda$, while rotating wake refers to the wake of an actuator disk with a finite value of $\lambda$.

2.3.4. Initial shape of the vortex sheet

The vortex sheet sheds from the circumference of the rotor, so it initially has a shape similar to the projected frontal area of the yawed disk, which is an ellipse with a semi-major axis of $R$ in the $z$-direction, and a semi-minor axis of $R\cos \beta$ in the $y$-direction. The disk-averaged velocity normal to the rotor $U_d$ is equal to $U_{\infty }\cos \beta (1-a)$, where $a$ is the turbine induction factor, and it is given by $0.5(1-\sqrt {1-C_T\cos ^{2}\beta })$ (Burton et al. Reference Burton, Sharpe, Jenkins and Bossanyi1995). Behind the turbine, the rotor streamtube area expands further as pressure recovers to the background value (Manwell, McGowan & Rogers Reference Manwell, McGowan and Rogers2010). At this location, the streamwise velocity is given by $U_{\infty }(1-2a)$ (Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2016; Shapiro et al. Reference Shapiro, Gayme and Meneveau2018). From continuity, $A_*$, the ratio of the expanded streamtube area to the projected frontal area of the rotor is therefore given by

(2.43)\begin{equation} A_*=\frac{\left(1-a\right)\cos\beta}{1-2a}\frac{1}{\cos\beta}= \frac{1+\sqrt{1-C_T\cos^{2}\beta}}{2\sqrt{1-C_T\cos^{2}\beta}}. \end{equation}

Neglecting the distance between the rotor and the end of the streamtube expansion, we set the initial wake area enclosed by the vortex sheet at $t=0$ to be the projected frontal area of the rotor times $A_*$. Therefore, $\xi _0(\theta )$ has an elliptical shape expressed by

(2.44)\begin{equation} \xi_{0}(\theta)=R\sqrt{A_*}\frac{|\cos{\beta}|}{\sqrt{1-\sin^{2}{\beta}\sin^{2}{\theta}}}. \end{equation}

For a small yaw angle $\beta$, the vortex sheet initially has an approximately circular shape. Therefore, one can approximate $\xi _0$ with $\tilde {\xi }_0$ given by

(2.45)\begin{equation} \tilde{\xi}_0\approx R\sqrt{A_*}. \end{equation}

2.4. Vortex sheet lateral deflection

The analytical solutions of the vortex sheet shape developed earlier are represented in the $(r,\theta )$ polar coordinate system, which is attached to the vortex sheet centre $C$. Therefore, in order to fully determine the locus of the vortex sheet with respect to a stationary coordinate system, we also need to compute how $y_c$ and $z_c$ vary with time or downstream distance (i.e. wake deflection). From figure 1(c), the Kutta–Joukowsky theorem can be used to obtain

(2.46)$$\begin{gather} \hat{y}_c=\frac{1}{2{\rm \pi}}\int_0^{\hat{t}} \int_0^{2{\rm \pi}}\hat{\gamma}(\theta',\hat{t})\sin\theta'\,\textrm{d}\theta'\,\textrm{d}\hat{t}, \end{gather}$$

(2.47)$$\begin{gather}\hat{z}_c=\frac{-1}{2{\rm \pi}}\int_0^{\hat{t}}\int_0^{2{\rm \pi}} \hat{\gamma}(\theta',\hat{t})\cos\theta'\,\textrm{d}\theta'\,\textrm{d}\hat{t}, \end{gather}$$

where $\hat {y}_c=y_c/\xi _0$ and $\hat {z}_c=z_c/\xi _0$. Inserting $\hat {\gamma }(\theta,\hat {t})$ from (2.26) into (2.46) and (2.47) and performing the integration lead to

(2.48)$$\begin{gather} \hat{y}_c=\frac{\hat{t}}{2}-\frac{\hat{t}^{3}}{96}, \end{gather}$$

(2.49)$$\begin{gather}\hat{z}_c=\frac{\chi\hat{t}^{4}}{384}. \end{gather}$$

From (2.49), the vertical displacement of $C$ is zero when $\chi =0$ (i.e. actuator disks with non-rotating wake) as expected from symmetry and consistent with prior experimental and numerical works (e.g.  Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2016; Howland et al. Reference Howland, Bossuyt, Martínez-Tossas, Meyers and Meneveau2016; Bartl et al. Reference Bartl, Mühle, Schottler, Sætran, Peinke, Adaramola and Hölling2018). Comparison of (2.48) and (2.49) for non-zero values of $\chi$ shows that $\hat {z}_c$ is non-zero but still considerably smaller than $\hat {y}_c$ for small values of $\hat {t}$. Therefore, we neglect $\hat {z}_c$ in this work for simplicity.

It is worth remembering that to derive the analytical solution for the deformation of the vortex sheet, we assumed that the shape of the vortex sheet does not largely deviate from a circle. Although this is an acceptable assumption for small values of yaw angle and time, it is less accurate for large values of time, when the vortex sheet rolls up and forms a CVP. As shown in figure 4(a), based on (2.48), the value of $\hat {y}_c$ may even decrease with an increase of $\hat {t}$, which is clearly unphysical. Since we expect that at large times (or downstream distances) a CVP can more realistically represent the vorticity shedding from the yawed rotor, we enhance our model for $\hat {y}_c$ so that at large distances it tends to the situation of a CVP instead of using the truncated series vortex sheet solution.

Figure 4. (a) Lateral deflection of the vortex sheet centre $y_c$ based on modelling the shed vorticity either as an approximately circular vortex sheet (2.48) or a CVP (2.53). The empirical relation (2.55) provides predictions similar to the former approach at small $\hat {t}$, while it tends to the latter solution at large times. (b) Schematic of modelling the CVP shedding from a yawed rotor.

Figure 4(b) shows a schematic of a CVP. The CVP has the circulation of $\varGamma _b$ (Shapiro et al. Reference Shapiro, Gayme and Meneveau2018) as shown in figure 4(b). The lateral position of the CVP is denoted by $y_{cvp}$, and the lateral distance between the wake centre $C$ and the CVP is denoted by $\delta _c$ in the figure. The vertical spacing between counter-rotating vortices is equal to $2\xi _0$ and it is assumed to remain constant. Initial values of $y_c$, $y_{cvp}$ and $\delta _c$ are zero. In this analysis, only the vorticity shed due to the yaw offset is considered as the effect of shed vorticity due to rotating blades on the wake deflection is expected to be small. Our objective is to find the variation of $y_c$ with time (or downstream distance), but let us first determine how $\delta _c$ varies with time. According to the Biot–Savart law,

(2.50)\begin{equation} v_c=\frac{2\varGamma_b}{2{\rm \pi} L}\cos\alpha=\frac{\varGamma_b \xi_0}{{\rm \pi} \left(\xi_0^{2}+\delta_c^{2}\right)}, \end{equation}

where $v_c$ is the lateral velocity of $C$ induced by the CVP, and $L$ and $\alpha$ are defined in figure 4(b). The CVP also moves with a lateral velocity of $v_{cvp}=\varGamma _b/(4{\rm \pi} \xi _0)$ due to its self-induced velocity. Therefore, one can write

(2.51)\begin{equation} \frac{\textrm{d}\delta_c}{\textrm{d}t}=v_c-v_{cvp}= \frac{\varGamma_b \xi_0}{{\rm \pi} \left(\xi_0^{2}+\delta_c^{2}\right)}- \frac{\varGamma_b}{4{\rm \pi}\xi_0}=\frac{\varGamma_b}{4{\rm \pi} \xi_0} \left(\frac{3\xi_0^{2}-\delta_c^{2}}{\xi_0^{2}+\delta_c^{2}}\right). \end{equation}

It is interesting to note that, according to (2.51), $\delta _c$ increases until $\delta _c$ approaches $\sqrt {3}\xi _0$. At this time $v_c=v_{cvp}$ and, therefore, the relative position of C with respect to the CVP does not change anymore, and $\delta _c$ remains equal to $\sqrt {3}\xi _0$ afterwards.

Next, we approximate $\varGamma _b \approx 2\xi _0\gamma _b$ in (2.51). We then integrate (2.51) (with separation of variables, $\delta _c$ and $t$) and write the solution in the dimensionless form. This yields an implicit expression for $\hat {\delta }_c(\hat {t})$ (it will be later expressed explicitly using an empirical formula),

(2.52)\begin{equation} \frac{\hat{t}}{2{\rm \pi}}=\frac{-2}{\sqrt{3}}\ln{\left(\frac{\sqrt{3}-\hat{\delta}_c}{\sqrt{3}+ \hat{\delta}_c}\right)}-\hat{\delta}_c, \end{equation}

where $\hat {\delta }_c=\delta _c/\xi _0$. Note that to derive (2.52), we assume that $\hat {\delta }_c < \sqrt {3}$. Given that $\hat {y}_c=\hat {\delta }_c+\hat {y}_{cvp}$, we have

(2.53)\begin{equation} \hat{y}_c=\hat{\delta}_c +\frac{\hat{t}}{2{\rm \pi}}. \end{equation}

Predictions of $\hat {y}_c$ based on modelling the vortex sheet as a CVP using (2.52) and (2.53) are shown in figure 4(a). As discussed earlier, the solution for $\hat {y}_c$ based on the CVP is expected to provide acceptable predictions at large values of $\hat {t}$, while the solution based on an approximately circular vortex sheet (2.48) works better at short values of $\hat {t}$. An empirical formula that merges these two behaviours to leading order in both limits can be written as

(2.54)\begin{equation} \hat{y}_c=\frac{c_1|\hat{t}|^{3}+c_2\hat{t}^{2}+c_3|\hat{t}|}{c_4\hat{t}^{2}+c_5|\hat{t}|+c_6}\textrm{sgn}(\hat{t}), \end{equation}

where $\textrm {sgn}(x)$ is the sign function of $x$, and $c_1,\ldots,c_6$ are polynomial coefficients, which need to be determined. Note that, similar to the vortex sheet and CVP solutions, the empirical relation is an odd function, so the wake deflection is opposite for turbines with opposite yaw angles. To find suitable values of polynomial coefficients ($c_1,\ldots,c_6$), we match the series expansion of (2.54) at $\hat {t}\to 0$ and $\hat {t}\to \infty$ with $\hat {t}/2-\hat {t}^{3}/96$ and $\hat {t}/(2{\rm \pi} )+\sqrt {3}$, respectively. This leads to a system of equations that needs to be solved. So we obtain

(2.55)\begin{equation} \hat{y}_c=\frac{({\rm \pi}-1)|\hat{t}|^{3}+2\sqrt{3}{\rm \pi}^{2}\hat{t}^{2}+48({\rm \pi}-1)^{2} |\hat{t}|}{2{\rm \pi}({\rm \pi}-1)\hat{t}^{2}+4\sqrt{3}{\rm \pi}^{2}|\hat{t}|+96({\rm \pi}-1)^{2}}\textrm{sgn}(\hat{t}). \end{equation}

Equation (2.55) provides predictions similar to (2.48) at small values of $\hat {t}$ and approaches (2.53) at large values of $\hat {t}$ as shown in figure 4(a).

2.5. Comparison with numerical simulations

In this section we compare predictions of the vortex sheet (i.e. wake edge) shape $\xi (\theta,t)$ based on the new proposed model with numerical simulation data. For simulations, the pseudo-spectral LES code LESGO is applied. The LESGO code has been used in prior works (Calaf, Meneveau & Meyers Reference Calaf, Meneveau and Meyers2010; Stevens, Graham & Meneveau Reference Stevens, Graham and Meneveau2014; VerHulst & Meneveau Reference VerHulst and Meneveau2015; Martínez-Tossas et al. Reference Martínez-Tossas, Churchfield, Yilmaz, Sarlak, Johnson, Sørensen, Meyers and Meneveau2018; Shapiro et al. Reference Shapiro, Gayme and Meneveau2018; Stevens, Martınez-Tossas & Meneveau Reference Stevens, Martınez-Tossas and Meneveau2018) to simulate flow past wind turbines and wind farms. It has been validated by detailed comparisons with several other LES codes (Martínez-Tossas et al. Reference Martínez-Tossas, Churchfield, Yilmaz, Sarlak, Johnson, Sørensen, Meyers and Meneveau2018). Turbines are simulated using the actuator disk model with rotation (ADM-R). See the Appendix C for more information about the LESGO code and the LES set-up of this study. Under uniform inflow conditions the role of turbulence is minimal, and the code runs mostly as an inviscid solver with regularization, as it was also used in Shapiro et al. (Reference Shapiro, Gayme and Meneveau2018). Simulations are performed for a range of local thrust coefficients $C_T'= 0.8, 1.0$ and $1.33$, yaw angles $\beta = 10^{\circ }$, $20^{\circ }$ and $30^{\circ }$, and rotation rates $\chi = 0$, 0.25 and 0.5. Note that according to (2.32) this means that $\gamma _b$ and $\hat {t}$ are negative and the curling is expected to be in the opposite direction of that shown in the sketch in figure 1, i.e. in the LES the wake is being deflected in the negative $y$-direction. Also, it is worth remembering that the rotation rate $\chi$ depends on both yaw angle $\beta$ and tip-speed ratio $\lambda$. According to (2.42), for a utility-scale wind turbine with a tip-speed ratio $\lambda =8$ and yaw angle $\beta =15^{\circ } - 30^{\circ }$, rotation rate $\chi$ varies between 0.25–0.5. The non-rotating case commonly used in the LES corresponds to an infinite tip-speed ratio and reduces to the standard actuator disk model (ADM) without rotation. The local thrust coefficient $C_T'$ is related to the thrust coefficient $C_T$ through $C_T=16C_T'/(4+C_T'\cos ^{2}\beta )^{2}$ (Shapiro et al. Reference Shapiro, Gayme and Meneveau2018). A fringe forcing region is used to force the flow back to laminar inflow when using periodic boundary conditions in the $x$-direction. Excluding this fringe region the effective domain has sides that are $L_x = 15.12D$, $L_y = 5.76D$ and $L_z = 5.76D$ long. A uniform grid with $N_x = 384$ effective grid points in the streamwise direction and $N_y = N_z = 192$ grid points in the spanwise and vertical directions are used. The centre of the actuator disk is placed 3.6D from the inlet of the domain.

In order to determine the shape of the wake edge based on the developed model, we need to first compute the value of $\hat {t}=\gamma _bt/\xi _0$, where $\xi _0$ can be approximated with $\tilde {\xi }_0$ (2.45) and $t=x/U_{con}$. Although the convection velocity $U_{con}$ in turbine wakes changes with the streamwise distance, it is approximated with a constant value in this study, as done in prior studies (e.g.  Shapiro et al. Reference Shapiro, Gayme and Meneveau2020). For cases with no incoming turbulence, the turbine wake does not significantly interact with the surrounding flow, and it experiences a slow recovery. In this case, the streamwise velocity profile in the central part of the wake can be modelled as a top-hat core (i.e. potential core) with a constant velocity $U_0$ equal to $U_{{in}}\sqrt {1-C_T\cos ^{2}\beta }$ (Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2016; Shapiro et al. Reference Shapiro, Gayme and Meneveau2018), where $U_{{in}}$ is the incoming velocity. The top-hat core is surrounded by a shear layer in which the velocity changes from $U_0$ to $U_{{in}}$. Therefore, we approximate the convection velocity with $U_{con}=0.5(U_0+U_{{in}})$. For instance, based on this definition, $|\hat {t}|=2$ (i.e. limiting value for using the analytical model) corresponds to a streamwise distance in the range of $12R-29R$ for a turbine with $C_T'=1.33$ (i.e. $C_T\approx 0.75$ for $\beta =0$), and $\beta =30^{\circ } - 10^{\circ }$.

The analysis presented in § 2 suggests that the wake shape, non-dimensionalised by $\xi _0(\theta )$, only depends on the dimensionless time $\hat {t}$, and the rotation rate $\chi =1/\lambda \sin \beta$. As a first test of the model, we compare the dimensionless model predictions with LES results normalised such that they can be presented as function of $\hat {t}$ and $\chi$. For the LES data, the edges of the wake are identified by tracking the edge of the streamtube that passes through the face of the actuator disk which is appropriate in this case due to the lack of turbulent mixing. Results are shown in figure 5, which shows that the LES data for turbines with different operating conditions approximately collapse onto the same wake profile curve for given values of $\hat {t}$ and $\chi$, in agreement with the proposed theory. The figure also shows that the proposed analytical model is able to capture the scaled wake shape. At larger time magnitudes and large rotation rates some discrepancies appear, especially in the bottom right quadrant. The governing equations are developed by assuming small deviations of the vortex sheet from its initial shape. In addition, a severely truncated series expansion is used to solve governing equations. Therefore, the model cannot fully capture the vortex roll-up and transition of the vortex sheet to a CVP at large times, and some discrepancies are noticeable. For non-zero values of rotation rate, an additional level of vorticity $\gamma _r$ sheds from the rotor circumference due to blade rotation. Given the sinusoidal nature of vorticity due to yaw offset (2.33), the additional shedding vorticity due to rotation increases the vorticity magnitude on either bottom or top halves of the wake (e.g.  bottom half for the data shown in figure 5), which in turn accelerates the vortex roll-up. At larger values of rotation rate, discrepancy is thus expected to be higher due to the earlier occurrence of vortex roll-up. This is confirmed in figure 5 by comparing model predictions at the same dimensionless time (e.g.  $\hat {t}=-1.6$), but different values of rotation rate.

Figure 5. Dimensionless shape of the wake of yawed wind turbines in uniform inflow for $C_T' = 0.8$ ($\circ$), $C_T' = 1$ ($\square$) and $C_T' = 1.33$ ($\triangle$) and yaw angles $\beta =10^{\circ }$ (red), $\beta =20^{\circ }$ (blue), $\beta =30^{\circ }$ (green) at various evolution times $\hat {t}$ and rotation rates $\chi$. The analytical model (black solid line) is shown for comparison. Note that results of $\beta =10^{\circ }$ are not shown for $\hat {t}=-1.6$ and $-2.0$ because for this yaw angle they correspond to downwind distances that exceed the computational domain.

Next, we expand the comparison and plot the results in terms of physical parameters that are more directly related to the flow configuration: downstream distance $x/R$, local thrust coefficient $C_T^{\prime }$, yaw angle $\beta$ and dimensionless rotation rate $\chi$. Figure 6 shows wake edge predictions of the analytical model together with the LES data for different values of $C_T'$, $\beta$ and $\chi$, at several downwind locations. The figure shows that the degree of wake curling increases with yaw angle, thrust coefficient and streamwise distance, as expected from the analysis presented in § 2. Moreover, wake rotation breaks the vertical symmetry of the wake. The results presented in figure 6 show that the wake shape depends strongly on all of the varied parameters: thrust coefficient, yaw angle and rotation rate. The analytical model developed is seen to agree well with the LES results, up to intermediate levels of wake curling. The analytical model successfully predicts the shape of the wake for various operating conditions. As the wake deformation grows further downstream or at increasing $C_T^{\prime }$, differences appear, as mentioned before due to the limitations of the model that is based on a severely truncated series expansion. It is clear that there is reduced agreement in the lower half of the wake for cases with large values of yaw angle and rotation rate. As discussed earlier, for these cases, the lower half of the wake cross-section is subject to a strong vortex roll-up caused by the cumulative vorticity due to both yaw offset and rotation. Still, the model is able to predict many qualitative features of the wake shape, including its vertical asymmetry for cases with rotation. In addition, the sideways displacement of the entire wake is also captured quite well in all cases.

Figure 6. Wake of yawed wind turbines in uniform inflow for yaw angles $\beta =10^{\circ }$ (red), $\beta =20^{\circ }$ (blue), $\beta =30^{\circ }$ (green) at various downstream locations $x/R$ and rotation rates $\chi$. Large-eddy simulation measurements are shown with symbols and modelled wake locations are shown with solid lines.

3.1. Vortex sheet evolution in turbulent atmospheric boundary layer

For the ideal flow case, we assumed that the vortex sheet does not decay as it moves downstream, and the strength of the streamwise vorticity only evolves along the vortex sheet in time following idealised vortex dynamics. Although this may be an acceptable assumption for turbines with uniform non-turbulent inflows, it is not expected to be valid for turbines immersed in turbulent environments such as the atmospheric boundary layer (ABL). The vertically varying mean inflow velocity in the ABL can be approximated as $U_\textrm {in}(z) = (u_*/\kappa )\ln (z/z_0)$, where $u_*$ is the friction velocity, $\kappa$ is the von-Kàrmàn constant and $z_0$ is the roughness height. As discussed in Shapiro et al. (Reference Shapiro, Gayme and Meneveau2020), the vorticity shedding from a yawed rotor decays due to the turbulent diffusion; i.e. $\gamma _b(t)$ and $\gamma _r(t)$ are functions of time. Instead of solving the full governing equations for a diffusing vortex sheet, we approximate the effects of diffusion by scaling $\gamma$ and the velocities by the circulation $\gamma _b(t)$ that is decaying according to a previously obtained analytical solution for the decay of CVP in a turbulent boundary layer (Shapiro et al. Reference Shapiro, Gayme and Meneveau2020). Specifically, we define new scaled variables $\hat {\gamma }$, $\hat {u}_r$ and $\hat {u}_\theta$ such that

(3.1a–c)\begin{equation} \gamma(\theta,t) = \gamma_b(t) \hat{\gamma}(\theta,\hat{t}),\quad {u}_r(\theta,t)= \gamma_b(t)\hat{u}_r(\theta,\hat{t}),\quad {u}_{\theta}(\theta,t)=\gamma_b(t) \hat{u}_{\theta}(\theta,\hat{t}). \end{equation}

If one defines scaled displacement and time according to

(3.2)$$\begin{gather} \hat{\xi}(\theta,\hat{t})=\frac{\xi(\theta,t)}{\xi_0}, \end{gather}$$

(3.3)$$\begin{gather}\hat{t}(t)=\frac{1}{\xi_0}\int_0^{t}\gamma_b(t')\,\textrm{d}t', \end{gather}$$

one recovers the original governing equations (2.15), (2.16) and (2.55) for the new definition of variables $\hat {\gamma }, \hat {u}_r, \hat {u}_\theta, \hat {\xi }$ and $\hat {t}$. An additional term is however introduced in (2.17), so

(3.4)\begin{equation} \frac{\partial \hat{\gamma}}{\partial \hat{t}}+\frac{1}{\hat{\xi}} \frac{\partial (\hat{\gamma} \hat{u}_\theta)}{\partial \theta}\approx{-} \frac{\xi_0\hat{\gamma}}{\gamma^{2}_b}\frac{\textrm{d} \gamma_b}{\textrm{d} t}. \end{equation}

It will be shown later in (3.8) that $\int (\gamma _b/\gamma _{b0})\,\textrm {d} t \propto c(1-\textrm {exp}(k_vt/c))$, where $\gamma _{b0}=\gamma _b(t=0)$, $c$ is a constant and $k_v$ is the expansion rate of the turbulent diffusive scale, and it is modelled as $u_*/U_{in}(z)$ (Shapiro et al. Reference Shapiro, Gayme and Meneveau2020). So the Taylor expansion of $\gamma _b(t)$ is given by

(3.5)\begin{equation} \gamma_{b}(t)\approx \gamma_{b0}(1+{O}(k_v)t+{O}(k^{2}_v)t^{2}). \end{equation}

For atmospheric flows, the value of $k_v$ at the hub height is equal to $u_*/U_h<\!<1$ (Shapiro et al. Reference Shapiro, Gayme and Meneveau2020). Therefore, $\gamma _b$ is a slow varying reference quantity, and the additional term on the right-hand side of (3.4) is neglected for simplicity. This means that the solution already developed in § 2.2 will still be used as a model for the decaying vortex sheet once rescaled by the prescribed $\gamma _b(t)$ evolution and using the modified time $\hat {t}$. Note that for a constant $\gamma _b$, (3.1a–c)–(3.3) become the same as dimensionless variables defined earlier for a non-decaying vortex sheet.

Next, in order to evaluate the integral of (3.3) and derive a relationship for $\hat {t}$ under turbulent inflow conditions we need to specify a convection velocity under turbulent inflow conditions. Due to atmospheric turbulence, the wake mixes and recovers more quickly than in the laminar inflow case, and, thus, the mean velocity at the wake edge is comparable to the incoming velocity. We therefore assume that in this case the vortex sheet at a height $z$ is convected downstream with the incoming velocity at that height $U_{{in}}(z)$; thus, $U_{con}=U_{{in}}(z)$ and $t \approx x/U_{{in}}(z)$.

To evaluate (3.3), we must specify the decay of vorticity, i.e. of $\gamma _b$ with streamwise distance. In Shapiro et al. (Reference Shapiro, Gayme and Meneveau2020) the decay of the total vortex circulation $\varGamma _b(x)$ was studied, its relationship with the density $\gamma _b$ being $\gamma _b=\varGamma _b/2R$. The resulting derived model for the total vortex circulation $\varGamma _b(x)$ as a function of downstream distance $x$ is given by

(3.6)\begin{equation} \frac{\varGamma_b(x)}{\varGamma_{b0}}=\frac{\sqrt{\rm \pi}}{4}\frac{R}{\eta(x)} \exp\left(-\frac{R^{2}}{8\eta^{2}(x)}\right)\left[I_0\left(\frac{R^{2}}{8\eta^{2}(x)}\right)+I_1 \left(\frac{R^{2}}{8\eta^{2}(x)}\right)\right], \end{equation}

where $\varGamma _{b0}=\varGamma _{b}(x=0)$, $I_n$ is the modified Bessel function of the first kind with order $n$, $\eta (x)=k_{\nu }(x-x_0)/24^{1/4}$ is the turbulent diffusive scale, and $x_0$ is the virtual origin assumed to be zero in the current work for simplicity.

One can use $\gamma _b=\varGamma _b/2R$, $t\approx x/U_\textrm {in}(z)$ and $\eta (x)\approx k_{\nu }x/24^{1/4}$ to rewrite (3.3) as

(3.7)\begin{equation} \hat{t}=\frac{24^{1/4}}{2k_{\nu}U_{{in}}(z)\xi_0 R }\int_0^{\eta}\varGamma_b(\eta')\,\textrm{d}\eta'. \end{equation}

Numerical integration of $\varGamma _b$, expressed by (3.6), yields results that can be conveniently approximated by the following (fitted) expression:

(3.8)\begin{equation} \frac{1}{\varGamma_{b0} R}\int_0^{\eta}\varGamma_b(\eta')\,\textrm{d}\eta' \approx 1.3\left[1-\textrm{exp}\left(-\frac{\eta(x)}{1.3R}\right)\right]. \end{equation}

We then use (2.32) to express $\varGamma _{b0}$ as a function of operating conditions, approximate $\xi _0(\theta )$ with $\tilde {\xi }_0$ given by (2.45), and insert (3.8) into (3.7) to obtain

(3.9)\begin{equation} \hat{t}(x,z)\approx{-}1.44 \frac{U_h}{u_*}\frac{R}{\tilde{\xi}_0} C_T \cos^{2}\beta \sin\beta \left[1-\textrm{exp}\left({-}0.35\frac{u_*}{U_\textrm{in}(z)}\frac{x}{R}\right)\right]. \end{equation}

Due to vorticity decay, $\hat {t}$ for turbulent inflow cases increases at a slower rate than the one for laminar inflow cases. For instance, according to (3.9), for a turbine with $C_T'=1.33$ subject to an ABL with $k_\nu = 0.05$, the streamwise position associated with $|\hat {t}|=2$ varies between $17R$ and $61R$ for $\beta =30^{\circ } - 10^{\circ }$. As mentioned in § 2.5, for the same turbine with $\beta =30^{\circ } - 10^{\circ }$ subject to a laminar flow, $|\hat {t}|=2$ at $x=12R-29R$. It is also worth noting that the above definition of $\hat {t}$, (3.9), is reduced to the one used for non-decaying vortex sheets (i.e. $\hat {t}=\gamma _{b0}t/\tilde {\xi }_0$) as $u_*$ tends to zero.

The effect of the ground on the wake deflection was not modelled in the uniform inflow cases. To model the effect of the ground, we use an image technique to modify the wake centre location $y_c$, as shown in figure 7. Modelling the vortex sheet as a CVP, the image CVP induces a lateral velocity in the opposite direction, termed as $v_g$ given by

(3.10)\begin{equation} v_g=\frac{\varGamma_b}{2{\rm \pi}} \left[\frac{1}{z+z_h-\xi_0} -\frac{1}{z+z_h+\xi_0}\right] = \frac{\varGamma_b \xi_0}{{\rm \pi}\left[(z+z_h)^{2}-\xi_0^{2}\right]}. \end{equation}

Therefore, the lateral wake deflection caused by the ground is given by

(3.11)\begin{equation} y_g(z) = \int_0^{t} v_g(z,t') \,\textrm{d} t' = \frac{\xi_0 \displaystyle\int_0^{t} \varGamma_b(t') \,\textrm{d} t}{{\rm \pi} \left[(z+z_h)^{2}-\xi_0^{2}\right]}. \end{equation}

Approximating $\varGamma _b \approx 2\xi _0\gamma _b$ and using $\hat {t} = \int _0^{t} \gamma _b(t') / \xi _0 \,\textrm {d} t'$, we find that

(3.12)\begin{equation} \hat{y}_g = \frac{2}{\rm \pi} \frac{\hat{t}}{\left[ ({z+z_h)}/{\xi_0}\right]^{2}-1}, \end{equation}

where $\hat {y}_g=y_g/\xi _0$. Substituting $\xi _0$ with $\tilde {\xi }_0$ in (3.12) for simplicity and subtracting this result from (2.55) yields

(3.13)\begin{equation} \hat{y}_c = \frac{({\rm \pi}-1)|\hat{t}|^{3}+2\sqrt{3}{\rm \pi}^{2}\hat{t}^{2}+48 ({\rm \pi}-1)^{2}|\hat{t}|}{2{\rm \pi}({\rm \pi}-1)\hat{t}^{2}+4\sqrt{3}{\rm \pi}^{2}|\hat{t}|+96({\rm \pi}-1)^{2}} \textrm{sgn}(\hat{t}) - \frac{2}{\rm \pi}\frac{\hat{t}}{[({z+z_h})/{\tilde{\xi}_0}]^{2}-1}. \end{equation}

It is worth mentioning that for $z_h\to \infty$, the second term on the right-hand side of (3.13) vanishes, and, thus, the equation is reduced to (2.55).

Figure 7. Schematic of modelling the effect of ground using an image technique.

3.2. Analytical model for mean velocity distribution

The shape of the wake edge predicted earlier can now be used to model the spatial distribution of the velocity deficit in the curled wake at each streamwise position. For non-curled (i.e. non-yawed) wakes, a number of wake profiles have already been proposed in the literature, including top hat (Katić, Højstrup & Jensen Reference Katić, Højstrup and Jensen1986; Frandsen et al. Reference Frandsen, Barthelmie, Pryor, Rathmann, Larsen, Højstrup and Thøgersen2006), Gaussian (Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2014; Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2016), double Gaussian (Schreiber, Balbaa & Bottasso Reference Schreiber, Balbaa and Bottasso2020) and super-Gaussian (Shapiro et al. Reference Shapiro, Starke, Meneveau and Gayme2019; Blondel & Cathelain Reference Blondel and Cathelain2020). Most of these conserve flux of momentum deficit only in its linearised version valid far downstream; see discussion in Bastankhah & Porté-Agel (Reference Bastankhah and Porté-Agel2014). Here we demonstrate the use of the shape deformation for an analytical model in the context of the Gaussian wake model (Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2014; Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2016), but it could be implemented in other wake models as well.

The modelled streamwise mean velocity,

(3.14)\begin{equation} U(\boldsymbol{x}) = U_{{in}}(z) - \Delta U(\boldsymbol{x}), \end{equation}

is defined based on the incoming velocity field $U_{{in}}(z)$ and the modelled velocity deficit in the wake, $\Delta U(\boldsymbol {x})$. In the Gaussian model, the velocity deficit profile is modelled as

(3.15)\begin{equation} \frac{\Delta U}{U_h} = C(x) \exp \left[ - \frac{(y-y_c)^{2} + (z-z_h)^{2}}{2\sigma^{2}} \right], \end{equation}

where $C(x)$ is the normalised maximum velocity deficit at each streamwise location and $\sigma$ is the characteristic wake width. In (3.15) the wake centre location is $y_c=\hat {y}_c\xi _0$, where $\hat {y}_c$ is obtained from (3.13) and $\xi _0$ is approximated with $\tilde {\xi }_0$ given in (2.45). In prior versions of the Gaussian wake model, it is assumed that the characteristic wake width depends only on downstream distance. Moreover, it is assumed that it grows linearly downstream at a rate $k$. The linear growth of the wake arises from the similarity solution when eddy viscosity is assumed to scale with a constant velocity, friction velocity $u^{*}$ and the wake scale $\sigma$ itself (Shapiro et al. Reference Shapiro, Starke, Meneveau and Gayme2019), i.e.

(3.16)\begin{equation} \sigma(x) = kx+\sigma_0, \end{equation}

where $k$ is the wake expansion rate and $\sigma _0$ the initial wake size. We now propose to include wake curling and deformation by making $\sigma$ dependent also on the angle $\theta$ according to

(3.17)\begin{equation} \sigma(x,\theta) = kx+0.4 \xi(\theta,x), \end{equation}

where $\xi (x,\theta )=\xi _0(\theta )\hat {\xi }(\theta,\hat {t})$, and $\xi _0(\theta )$ is given by (2.44), and the dimensionless wake shape $\hat {\xi }(\theta,\hat {t})$ is given by either the analytical relation of (2.29) for $|\hat {t}|\leq 2$ or the empirical relation (B1) for given values of polar angle $\theta$ and dimensionless time $\hat {t}$. The polar angle is determined at each position from $\tan \theta =(z-z_h)/(y-y_c)$, and $\hat {t}$ is given by (3.9). In (3.17) the first term on the right-hand side of the equation expands the wake in all radial directions due to turbulent mixing, while the second term deforms the wake cross-section according to the vortex sheet solution derived earlier. According to (3.17), for an unyawed turbine, the initial characteristic wake width is reduced to $0.4R\sqrt {A_*}$, which is the same as the one suggested by Bastankhah & Porté-Agel (Reference Bastankhah and Porté-Agel2014). For the wake expansion rate, we assume that $k = \alpha u_*/U_{{in}}(z)$ (Shapiro et al. Reference Shapiro, Starke, Meneveau and Gayme2019), where $\alpha$ is an empirical constant. Alternatively, $k$ can be estimated based on the turbulence intensity $I$ of the incoming boundary layer flow (i.e. $k=\alpha ' I$, where $\alpha '$ is an empirical constant) as suggested in prior studies (see Niayifar & Porté-Agel (Reference Niayifar and Porté-Agel2016), Carbajo Fuertes, Markfort & Porté-Agel (Reference Carbajo Fuertes, Markfort and Porté-Agel2018), Zhan, Letizia & Iungo (Reference Zhan, Letizia and Iungo2020), among others). At each height $z$, turbulence intensity $I(z)$ is defined as $\sqrt {\overline {u^{'2}}}/U_{in}(z)$, where $u'$ is the turbulent fluctuation of streamwise velocity and the overbar denotes time averaging. Note that by invoking the logarithmic law for the fluctuating velocity variance in high-Reynolds-number turbulent boundary layers (Marusic et al. Reference Marusic, Monty, Hultmark and Smits2013; Meneveau & Marusic Reference Meneveau and Marusic2013), one can show that $u_*/U_{in}$ and $I$ are related to each other by $I=(u_*/U_{in}) [B_1-A_1\ln (z/\delta )]^{1/2}$, where $\delta$ is the boundary layer thickness, and $A_1$ and $B_1$ are constants.

The maximum velocity deficit $C(x)$ in (3.15) is obtained by enforcing the conservation of streamwise momentum deficit flux $\rho \int \Delta U (U_h - \Delta U) \,\textrm {d} A \approx T \cos \beta$. To simplify the integration and avoid dependence on $\theta$, we approximate $\sigma ^{2}(x,\theta )$ with $\tilde {\sigma }^{2}(x)$ where the latter is given by

(3.18)\begin{equation} \tilde{\sigma}^{2}(x) = (kx + 0.4\tilde{\xi}_0)(kx + 0.4\tilde{\xi}_0\cos\beta). \end{equation}

In this expression, $\tilde {\xi }_0$ is given by (2.45) and $k=\alpha u_*/U_h$ is the wake expansion rate at $z=z_h$. This yields

(3.19)\begin{equation} C(x) = 1 - \sqrt{1 - \frac{C_T \cos^{3} \beta}{2 \tilde{\sigma}^{2}(x)/R^{2}}}. \end{equation}

Note that in stating conservation of flux of streamwise momentum deficit, pressure and turbulent and viscous shear stress effects are assumed to be negligible. This may be a questionable assumption in the near wake region as well as far wake of turbines deep inside a wind farm (Bastankhah et al. Reference Bastankhah, Welch, Martínez-Tossas, King and Fleming2021).

For the sake of completeness, a summary of the steps required to implement the proposed model and predict wake velocity deficit distributions at a given downwind location $\boldsymbol {x}=(x,y,z)$ is provided below.

1. (1) Compute the approximate form of the initial wake shape $\tilde {\xi }_0$ (2.45).

2. (2) Determine the dimensionless time $\hat {t}$ from (3.9).

3. (3) Find the wake centre location $y_c\approx \hat {y}_c\tilde {\xi }_0$, where $\hat {y}_c$ is given by (3.13).

4. (4) Find the polar angle $\theta$, which is measured from the positive $y$-axis toward the positive $z$-axis such that $\tan \theta =(z-z_h)/(y-y_c)$.

5. (5) Evaluate the initial wake shape $\xi _0(\theta )$ (2.44).

6. (6) Calculate the wake shape function $\xi (\theta,x)=\xi _0(\theta )\hat {\xi }(\theta,\hat {t})$, where $\hat {\xi }(\theta,\hat {t})$ can be estimated either from the analytical solution (2.29) (if $|\hat {t}|<2$) or the empirical one (B1), and $\chi$ is given by (2.42).

7. (7) Find the wake width $\sigma (x,\theta )$ based on (3.17).

8. (8) Evaluate the maximum velocity deficit $C(x)$ from (3.19), where $\tilde {\sigma }^{2}(x)$ is given by (3.18).

9. (9) Determine the wake velocity deficit $\Delta U$ according to (3.15), where $\sigma$ obtained in step (7) is used.

3.3. Comparison with LES

In the following, the streamwise mean velocity distribution based on the proposed model is compared with the LES data for turbulent ABL inflow cases. Simulations of yawed wind turbines represented as rotating actuator disks (ADM-R) are performed using $C_T' = 1.33$ and a local tip-speed ratio of $\lambda '=10.67$ at yaw angles of $\beta = 15^{\circ }$, $20^{\circ }$, $25^{\circ }$ and $30^{\circ }$, where the local tip-speed ratio is defined as $\lambda '=\varOmega R/U_d$. In unyawed conditions the selection of local thrust coefficient and tip-speed ratio corresponds to $C_T = 0.75$ and $\lambda =8$, which are realistic of modern utility-scale turbines. The actuator disk with diameter $D = 100$ m and a hub height of $z_h = 100$ m is placed 500 m from the inlet of a domain with an effective size of $L_x = 3.75$, $L_y = 3$ and $L_z = 1$ km. The domain is divided into $N_x = 360$, $N_y = 288$ and $N_z = 432$ grid points. The velocity field is averaged for a time $\mathcal {T} u_*/L_z \approx 8$, where $u_* = 0.45$ m s$^{-1}$. A $0.49L_z$ shift is used to reduce streamwise streaks in the time-averaged velocity field. The roughness height is $z_0 = 0.1$ m.

Figures 8 show contour maps of streamwise velocity on representative planes across the LES domain. It shows the streaks in the turbulent ABL at the turbine hub height and the generated wake behind the yawed wind turbine, which is shown as a black circle. At downstream locations, the effect of the curled yawed wind turbine wake on the cross-plane velocity field is shown at downstream locations of $x/R = 8$, $24$ and $40$.

Figure 8. Contour plots of instantaneous streamwise velocity including a wind turbine with turbulent boundary layer inflow from LES. Turbine operating parameters are $C_T' = 1.33$ and yaw angle $\beta =25^{\circ }$. Contours are shown through the turbine centre at $z=z_h$, at the back of the domain at $y=L_y$ and $x=L_x$ and at cross-planes of $x/R=8$, $24$ and $40$. The swept area of the rotor is denoted as a black circle. A zoomed in flow field around the turbine (red box) is also shown and white arrows highlight the sense of rotation of the induced CVP.

Figures 9 and 10 show the wake mean velocity distributions based on the LES results and the analytical model for $\beta = 15^{\circ }$ and $25^{\circ }$ at various downstream locations. Top panels show the normalised velocity deficit $\Delta U/U_h$ and bottom panels show the normalised streamwise velocity distribution.

Figure 9. (a) Contour plots of normalised wake velocity deficit behind a wind turbine in turbulent inflow with a thrust coefficient of $C_T' = 1.33$ and local tip-speed ratio $\lambda '=10.67$ at a yaw angle of $\beta =15^{\circ }$. White circles indicate the frontal area of wind turbines. (b) Contour plots of normalised streamwise velocity behind the same turbine.

Figure 10. Same as figure 9 but for a yaw angle of $\beta =25^{\circ }$.

As in past studies (Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2016; Shapiro et al. Reference Shapiro, Gayme and Meneveau2018), the wake recovery rate $k$ in the analytical model is calculated by fitting a Gaussian profile to the downstream wake profiles at $z=z_h$. This gives the resulting wake expansion rate as $k=0.6 u_*/U_{in}$. As seen in figures 9 and 10, the model captures the curling and deflection of the wake as well as some variation in the wake deflection as a function of vertical distance due to ground effects. The LES results display slightly less curling than the model as well as more noticeable wake deflection towards the ground that is opposite to the deflection of the bulk of the wake. Also, the analytical model does not predict small vertical wake deflections observed in the LES data. According to the vortex sheet analysis performed in the current study (see (2.49)), this vertical deflection is due to rotation effects (i.e. non-zero values of rotation rate $\chi$). As mentioned in § 2.4, for simplicity, we did not include the vertical wake deflection in the final version of the model. Despite these small differences, figures 9 and 10 show that overall model predictions are in acceptable agreement with the LES data.

Wake flow results are also used to compute the power of a downwind turbine based on both the proposed analytical model and the LES data. In both cases, a virtual wind turbine is placed in the flow field at different distances downstream of the yawed turbine. The idea is that as yawing increases, the downstream turbine will overlap less and less with the wake due both to the sideways displacement of the wake and to the curled crescent shaped wake that creates lower velocity deficit at the centre of the hypothetical downstream turbine.

In order to evaluate the power generated by the downwind turbine, we require the disk-averaged streamwise velocity defined according to

(3.20)\begin{equation} U_d = (1-a)\frac{1}{{\rm \pi} R^{2}} \iint _{disk} U(x_T,y,z) \,\textrm{d} z\,{\textrm{d}y}, \end{equation}

where $U(x_T,y,z)$ is the mean velocity in the flow at the turbine location $x_T$ computed from the model or from the LES and the integration covers the turbine disk area. For the latter case, we evaluate the mean velocity by time averaging $U(x_T,y,z)=\langle \tilde {u}_1(x_T,y,z)\rangle$. Since the turbine is not included in the simulation, the turbine disk velocity that would occur there includes the $(1-a)$ prefactor, where $a$ is the turbine's assumed induction factor. The power is subsequently calculated as

(3.21)\begin{equation} P = \tfrac{1}{2} \rho {\rm \pi}R^{2} C_T' U_d^{3}, \end{equation}

and it is normalised by the power that an unyawed free-standing turbine would generate under similar conditions, $P_0 = (1/2) \rho {\rm \pi}R^{2} C_T'(1-a)^{3} U_h^{3}$ (thus making the result independent of assumed $C_T'$, etc.).

Figure 11 shows the normalised power of this hypothetical turbine as a function of streamwise spacing for different yaw angles of the upwind turbine. The figure demonstrates good agreement between the analytical model and the LES results for a broad range of streamwise spacings and yaw angles.

Figure 11. (a) Sketch of a hypothetical turbine placed at various locations downstream of the yawed turbine in the LES field. (b) Normalised power of the hypothetical turbine operating at $C_T' = 1.33$. Different yaw angles in LES are $\beta =15^{\circ }$ ($\circ$, red), $\beta =20^{\circ }$ ($\square$, green), $\beta =25^{\circ }$ ($\triangle$, blue) and $\beta =30^{\circ }$ (, purple). The predictions based on the analytical curled wake model are shown as solid lines.

3.4. Comparison with experimental data

Model predictions are also compared with wind-tunnel experiments by Bastankhah & Porté-Agel (Reference Bastankhah and Porté-Agel2016). Flow measurements were performed to quantify the wake of a yawed wind turbine with a diameter of $15$ cm and hub height of $12.5$ cm, and the turbine is subject to a turbulent boundary layer, naturally developed over the smooth surface of the wind-tunnel floor. Additional information on turbine properties ($C_T$, $\lambda$, etc.) and inflow conditions ($U_h$, $I$, etc.) may be found in Bastankhah & Porté-Agel (Reference Bastankhah and Porté-Agel2016). The wake recovery rate $k$ for the analytical model is estimated based on $k=0.35I$ (Carbajo Fuertes et al. Reference Carbajo Fuertes, Markfort and Porté-Agel2018). Figure 12 shows contours of normalised velocity deficit in $yz$-planes at different downwind locations and different yaw angles based on both experiments and model predictions. Overall the figure shows that the proposed model is able to successfully predict the complex curled shape of the wake and its lateral deflection. The opposite wake deflection close to the ground is also well captured by the model. While the overall trends and qualitative features of the velocity defect distribution of the curled wake are reproduced by the analytical model, some differences between measurements and model can still be discerned. The figure shows that the vertical extent of the lower half of the wake is underestimated by the analytical model, which is mainly due to neglecting the wake of the turbine tower in the analytical model. Tower wake effects are however expected to be less significant for utility-scale wind turbines which tend to have less bulky towers (with respect to the rotor diameter).

Figure 12. Contours of normalised velocity deficit in $yz$-planes at different downwind locations and different yaw angles based on: wind-tunnel experiments (Bastankhah & Porté-Agel Reference Bastankhah and Porté-Agel2016) and the new proposed analytical model. White circles indicate the frontal area of wind turbine.