2.1. Modelling framework and assumptions

The proposed modelling framework is presented schematically in figure 1. The schematic distinguishes between data sets, denoted by numerical indicators (1–8), and mathematical models, represented by alphabetical labels (a to d). The inputs, main body of the model and the outputs are explained as follows.

1. (i) Input. To evaluate the wake response to a specific AWC approach, the two inputs of the model are (1) the AWC approach to be evaluated and (2) a non-controlled baseline wake.

2. (ii) Main body. The main body of the model starts with (a) a forcing model for transferring the AWC into (3) equivalent periodic rotor forcing. Then an iterative loop starts, involving three interconnected components, including: (b) a resolvent model for predicting (4) the coherent structures in the wake ( $\hat {\boldsymbol{u}}$ ); (c) a RS model that predicts (5) the RS generated by both coherent structures ( $\boldsymbol{\tau }^c$ ) and small-scale turbulence ( $\boldsymbol{\tau }^s$ ); (d) a mean flow model to compute (6) the mean flow correction, i.e. the enhanced wake recovery ( $\delta \boldsymbol{u}$ ), which is added back to (2) the non-controlled time-averaged wake to produce the baseflow for the resolvent analysis in the next iteration.

3. (iii) Output. Upon achieving convergence of the iteration loop, the model yields self-consistent predictions of (7) the AWC-induced coherent structures and (8) the time-averaged wake.

The main assumptions for the modelling include: (i) the periodicity of forcing and response, i.e. both are simple harmonic oscillations in time, required by the resolvent model; (ii) the axisymmetry of the baseflow, i.e. no shear or veer of incoming wind are considered. The former assumption facilitates the solution in frequency domain, and the latter reduces the wake to a two-dimensional problem for enhancing the model’s efficiency.

In the present work, an LES with a wind turbine modelled with an actuator surface method (see § 3.1) is conducted to obtain the non-controlled baseline wake and the time-averaged aerodynamic forces on the rotor.

2.2. Governing equations

In this work, we seek to predict the mean and the large-scale coherent fluctuations of the wake. To this end, we employ the triple decomposition (Hussain & Reynolds Reference Hussain and Reynolds1970; Yim et al. Reference Yim, Meliga and Gallaire2019) to separate the different components of the flow variables as follows:

(2.1) \begin{align} \boldsymbol{u} & = \overline {\boldsymbol{u}} + \boldsymbol{u}^{c}(t) + \boldsymbol{u}^{s}(t) = \widetilde {\boldsymbol{u}}(t) + \boldsymbol{u}^{s}(t), \\[-12pt] \nonumber \end{align}

(2.2) \begin{align} p & = \overline {p} + p^{c}(t) + p^{s}(t) = \widetilde {p}(t) + p^{s}(t), \\[-12pt] \nonumber \end{align}

(2.3) \begin{align} \boldsymbol{f} & = \overline {\boldsymbol{f}} + \boldsymbol{f}^{c}(t) + \boldsymbol{f}^{s}(t) = \widetilde {\boldsymbol{f}}(t) + \boldsymbol{f}^{s}(t), \end{align}

where $\boldsymbol{u}$ , $p$ , $\boldsymbol{f}$ represent velocity, pressure and rotor forcing, respectively. The symbol $\overline {\{\cdot \}}$ represents the time-averaged part, $\{\cdot \}^{c}$ denotes the coherent fluctuating part, $\{\cdot \}^{s}$ denotes the small-scale turbulent part. Here $\widetilde {{\{\cdot \}}}$ denotes the filtered field containing both mean and coherent fluctuation, while the small-scale turbulent component is excluded. As the AWC focuses on the mean and the large-scale coherent fluctuation of the wake, we focus on the filtered fields and assume that they verify the incompressible Navier–Stokes equation as follows:

(2.4) \begin{align} \frac {\partial \widetilde {\boldsymbol{u}}}{\partial t}+ (\widetilde {\boldsymbol{u}} \cdot \nabla) \widetilde {\boldsymbol{u}} + \widetilde {p} - \nu _{\textit {eff}} \nabla ^2 \widetilde {\boldsymbol{u}} & = \widetilde {\boldsymbol{f}}, \end{align}

(2.5) \begin{align} \nabla \cdot \widetilde {\boldsymbol{u}} &= 0. \end{align}

Here, an effective viscosity $\nu _{\textit {eff}}$ is introduced to account for the influence of small-scale turbulence (Symon et al. Reference Symon, Madhusudanan, Illingworth and Marusic2023). The closure will be provided in $\S$ 2.5 through an eddy viscosity model.

By assuming that both the coherent parts of forcing ( $\boldsymbol{f}^c$ ) and response ( $\boldsymbol{u}^c$ , ${p}^c$ ) are simple harmonic, they are written as

(2.6) \begin{equation} [ \boldsymbol{u}^c(t), p^c(t), \boldsymbol{f}^c(t)] = \frac {1}{2}[\hat {\boldsymbol{u}}, \hat {p}, \kern3pt \hat {\boldsymbol{\kern-3pt f}}] \text{e}^{-\mathrm{i} \omega t} + \frac {1}{2} [\hat {\boldsymbol{u}}^H, \hat {p}^H, \kern3pt \hat {\boldsymbol{\kern-3pt f}}^H] \text{e}^{\mathrm{i} \omega t} \end{equation}

where $\mathrm{i} = \sqrt {-1}$ is the imaginary unit, $\omega$ is the angular frequency and $\{\cdot \}^H$ denotes complex conjugate. The spatial modes are represented by $\hat {\{\cdot \}}$ , while $\text{e}^{\pm \mathrm{i} \omega t}$ describes the temporal oscillation.

Substituting this ansatz into (2.5) and separating coherent fluctuating and time-averaged components yields two coupled systems. The dynamics of the coherent fluctuation at forcing frequency $\omega$ is described by the harmonic linearized Navier–Stokes equation,

(2.7) \begin{align} -\mathrm{i} \omega \hat {\boldsymbol{u}} + (\overline {\boldsymbol{u}} \cdot \nabla) \hat {\boldsymbol{u}} + (\hat {\boldsymbol{u}} \cdot \nabla) \overline {\boldsymbol{u}} + \nabla \hat {p} - \nu _{\textit {eff}} \nabla ^2 \hat {\boldsymbol{u}} & = \kern3pt \hat {\boldsymbol{\kern-3pt f}}, \\[-12pt] \nonumber \end{align}

(2.8) \begin{align} \nabla \cdot \hat {\boldsymbol{u}} & = 0. \end{align}

These equations govern the resolvent model (Li & Yang Reference Li and Yang2024), which will be detailed in § 2.4.

For the time-averaged field, the governing equations are

(2.9) \begin{align} \left (\overline {\boldsymbol{u}} \cdot \nabla \right) \overline {\boldsymbol{u}} + \nabla \overline {p} & = \overline {\boldsymbol{f}} + \nabla \cdot \boldsymbol{\tau ^c} + \nu _{\textit {eff}} \nabla ^2 \overline {\boldsymbol{u}}, \\[-12pt] \nonumber \end{align}

(2.10) \begin{align} \nabla \cdot \overline {\boldsymbol{u}} &= 0. \end{align}

Here, $\boldsymbol{\tau ^c} = \tau ^c_{ij} = -\overline {u^c_{i}u^c_{j}}$ represents the RS tensor arising from AWC-induced coherent fluctuation ( $\boldsymbol{u}^c$ ), where indices $i,j \in \{1,2,3\}$ denote velocity components in three-dimensional space. Once coherent fluctuations are excited by AWC, $\nabla \cdot \boldsymbol{\tau ^c}$ drives mean flow modification. Here $\nu _{\textit {eff}} \nabla ^2 \overline {\boldsymbol{u}}$ represents the effect of small-scale turbulence on the mean wake. These two terms are included in the RS model and will be explained in § 2.6.

The above two equation systems are interconnected: the coherent fluctuation of the wake is computed with (2.7) and (2.8) based on the mean wake $\overline {\boldsymbol{u}}$ ; and the mean wake is also influenced by the coherent fluctuation through RS ( $\boldsymbol{\tau }^c$ ) in (2.9).

2.3. Rotor forcing model

In the present framework, the rotor forcing model serves to transfer the different AWC approaches to coherently oscillatory forcing $\kern3pt \hat {\boldsymbol{\kern-3pt f}}$ , which will be employed as an input of the resolvent model when solving (2.7). Here, we develop analytical forcing models for three different AWC approaches that induce helical-type meandering, namely the DIPC, in-plane circular translation and precession of the rotor, as shown in figure 3. The first approach introduces wake perturbation through DIPC, which is proposed by Frederik et al. (Reference Frederik, Doekemeijer, Mulders and van Wingerden2020) and is known as the helix approach. The latter two approaches perturb the wake via rotor displacements, which exist in FOWT due to platform motion (Tran & Kim Reference Tran and Kim2015). Following the axisymmetry assumption of the time-averaged wake, we model only the global effects at the rotor disk level rather than individual blade effects, and assume small perturbation amplitudes to enable linearization of rotor aerodynamics about the base state.

Figure 3.

Schematic of AWC approaches: (a) DIPC; (b) in-plane circular translation; (c) precession of the rotor. Shaded areas in (b) and (c) indicate the original rotor swept area.

2.3.1. Dynamic individual pitch control

As shown in figure 3(a), the DIPC approach periodically adjusts the pitch angle $\beta$ of each blade at a common frequency $\omega _p$ , with well-designed phase differences between blades, as follows:

(2.11) \begin{equation} \beta _k(t) = A_\beta \sin \left(\omega _p t + \frac {2\pi }{3}(k-1) \right), \end{equation}

where $k \in \{1,2,3\}$ is the blade index and $A_{\beta }$ is the angular amplitude of pitch control. The angular velocity of the pitch variation is $\omega _p = \Omega \pm \omega$ with $\Omega$ the angular velocity of the rotor. This control strategy leads to periodic variation of the angle of attack (AoA) and the aerodynamic forces. The first blade points upwards (in the $z$ direction) at $t = 0$ , and its azimuth angle is $\psi _1 = \Omega t$ . According to the multiblade coordinate transformation (Bir Reference Bir2008; Van der Hoek et al. Reference van der Hoek, den Abbeele, Simao Ferreira and van Wingerden2024), DIPC introduces an equivalent rotor level rotating radial force (Korb et al. Reference Korb, Asmuth and Ivanell2023). A linear relationship between the force and the pitch angle variations has been proposed (Korb et al. Reference Korb, Asmuth and Ivanell2023), but the coefficient preceding the linear relation is undetermined. Here we propose an approach to determine this coefficient by using simply the blade aerodynamic data obtained from the non-controlled baseline case.

The AoA of a blade under DIPC varies as a function of time $(t)$ and radial location $(r)$ :

(2.12) \begin{equation} \alpha _k (t,r) = \alpha ^0(r) - \beta _k(t). \end{equation}

Here $\alpha ^0 (r)$ is the AoA of the baseline case without AWC and $\beta _k(t)$ is the pitch variation defined in (2.11). Following Korb et al. (Reference Korb, Asmuth and Ivanell2023), we consider only the contribution from the lift, whose variation due to DIPC can be expressed as

(2.13) \begin{equation} \begin{aligned} \frac {\delta L_k(t,r)}{\overline {L^0}(r)} & = \frac { C_L(\alpha _k(t,r)) - C_L(\alpha ^0(r))}{C_L(\alpha ^0(r))} \\ & \approx - \frac {\beta _k(t)}{C_L(\alpha ^0(r))} \left.\frac {\partial C_L}{\partial \alpha }\right |_{\alpha ^0(r)}. \\ \end{aligned} \end{equation}

Here, $\overline {L^0}$ denotes the baseline time-averaged lift force and $\delta L_k$ denotes the lift variation. This expression employs a linearization around the baseline state using Taylor expansion, assuming a small pitch variation. By using the geometric similarity of force components, this relationship can be projected into the streamwise $(x)$ and azimuthal $(\theta)$ directions to obtain the thrust and tangential force variations:

(2.14) \begin{equation} \frac {\delta F_{x,k}(t,r)}{\overline {F^0_{x}}(r)} = \frac {\delta F_{\theta,k}(t,r)}{\overline {F^0_{\theta }}(r)} = \frac {\delta L_k(t,r)}{\overline {L^0}(r)} =- \frac {\beta _k(t)}{C_L(\alpha ^0(r))} \left.\frac {\partial C_L}{\partial \alpha }\right |_{\alpha ^0(r)}. \end{equation}

Once the baseline AoA ( $\alpha ^0(r)$ ) and the baseline blade forces ( $\overline {F^0_x}(r)$ , $\overline {F^0_\theta }(r)$ ) are obtained, the force variation can be quantitatively computed with the derivative of the lift coefficient $ (({\partial \overline { C_L}}/{\partial \alpha }) |_{\alpha ^0(r)})$ from the airfoil database and the pitch variation ( $\beta _k(t)$ ).

Next, the aerodynamic force variation is transformed from the blade-following reference frame to a ground-fixed one, to model the effect of DIPC on a stationary actuator disk exerting equivalent unsteady normal and tangential forces. To achieve this, we assume that the rotor’s rotational frequency is significantly larger than the forcing frequency ( $|\Omega | \gg |\omega |$ ). This assumption allows us to neglect the time interval between blade passages and to average the forces on the rotor swept area. At any point $(r, \theta)$ on the rotor swept area, the DIPC-controlled rotor, as described in (2.11), generates a streamwise blade force at three blade passages per revolution ( $t_1$ , $t_2$ , $t_3$ ) as follows:

(2.15) \begin{equation} \delta F_{x,k}(t_k,r) = - \frac {A_\beta \overline {F^0_x}(r)}{C_L(\alpha _0(r))} \left.\frac {\partial C_L}{\partial \alpha }\right |_{\alpha _0(r)} \sin { \left ( \left (\Omega \pm \omega \right) t_k + \frac {2\pi }{3}(k-1) \right) } \end{equation}

with $k \in {1,2,3}$ the blade index. The blade passage instant is expressed as

(2.16) \begin{equation} t_k = \frac {\theta - \dfrac {\pi }{2}- \dfrac {2\pi }{3}(k-1) }{\Omega }. \end{equation}

Combining these equations, yields

(2.17) \begin{equation} \delta F_{x,k}(t_k,r) = \frac {A_\beta \overline {F^0_x}(r)}{C_L(\alpha _0(r))} \left.\frac {\partial C_L}{\partial \alpha }\right |_{\alpha _0(r)} \cos { \left (\theta \pm \omega t_k \right)}. \end{equation}

By assuming $t_k \approx t$ and averaging for a period of rotor evolution, the streamwise coherent force variation per unit area exerted by the equivalent actuator disk on the flow can be expressed as

(2.18) \begin{equation} f^c_x = -A_\beta \cos (\theta \pm \omega t) \frac {1}{2} U_\infty ^2 \frac {\overline {C^0_x}(r)}{C_L(\alpha _0(r))} \left.\frac {\partial C_L}{\partial \alpha }\right |_{\alpha _0(r)}, \end{equation}

where $\overline {C^0_x}(r)$ denotes the local thrust coefficient at radial position $r$ , defined as

(2.19) \begin{equation} \overline {C^0_x}(r) = \frac {3 \cdot \overline {F^0_x}(r) \text{d} r}{ 2\pi r \text{d} r \cdot \frac {1}{2} U_\infty ^2} = \frac {3 \overline {F^0_x}(r) }{ \pi r U_\infty ^2 }. \end{equation}

Note that the sign of the force on the blade and that on the flow are opposite.

Similarly, the coherent tangential force variation per unit area exerted by the actuator disk can be calculated as

(2.20) \begin{equation} f^c_\theta = -A_\beta \cos (\theta \pm \omega t) \frac {1}{2} U_\infty ^2 \frac {\overline {C^0_\theta }(r)}{C_L(\alpha _0(r))} \left.\frac {\partial C_L}{\partial \alpha }\right |_{\alpha _0(r)}, \end{equation}

with $\overline {C^0_\theta } (r)$ the time-averaged tangential force coefficient of the baseline case. By assuming radial force variation is negligible, the rotor force variation can be expressed as

(2.21) \begin{equation} \boldsymbol{f^c}(t) = -\frac {A_\beta \cos (\theta \pm \omega t) U_\infty ^2}{2 C_L(\alpha _0(r))} \left.\frac {\partial C_L}{\partial \alpha }\right |_{\alpha _0(r)} \left [\overline {C^0_x}(r), 0, \overline {C^0_\theta }(r) \right]. \end{equation}

The frequency-domain transformation of these forces reveals their spatial distribution, as shown in figure 4(a–c). The streamwise component dominates and exhibits antisymmetry about the rotor diameter, generating an out-of-plane moment in the wake flow (Van der Hoek et al. Reference van der Hoek, den Abbeele, Simao Ferreira and van Wingerden2024). It rotates around the rotor’s axis over time, providing initial perturbation to the wake.

Figure 4.

The spatial modes of unsteady rotor forcing on the actuator disk for DIPC (a–c), translation (d–f),and precession (g–i). Panels (a,d,g), (b,e,h) and (c,f,i) display the streamwise, radial and azimuthal components, respectively. The forcings are normalized by $A U^2_\infty$ , where $A$ is the amplitude of each control strategy. Yellow arrows indicate the rotating direction of the rotor.

2.3.2. Circular translation

In-plane circular translation involves circular movement of the rotor within its rotation plane (figure 3 b). The hub centre traces a circular path of radius $A_t$ in the $yz$ plane relative to its original position ( $O_0$ ). This motion is equivalent to a combination of synchronized sway and heave movements, with horizontal and vertical displacements expressed as

(2.22) \begin{align} \delta y (t) &= A_{t}\sin (\omega t), \\[-12pt] \nonumber \end{align}

(2.23) \begin{align} \delta z (t) &= \pm A_{t}\cos (\omega t). \\[6pt] \nonumber \end{align}

Here, the sign of the vertical displacement $\delta z (t)$ will define the helix direction of the wake. The associated unsteady rotor forcing can be determined straightforwardly by using the motion-to-forcing model proposed in Li & Yang (Reference Li and Yang2024), which computes the aerodynamic force variation due to harmonic rotor position oscillations. The resulting thrust and tangential force oscillations are

(2.24) \begin{equation} f^c_x = \frac {U_\infty ^2 }{2} \frac {\partial {\overline {C^0_x}}(r)}{\partial r} \left ( \delta y(t) \cos {\theta } + \delta z(t) \sin {\theta } \right) = A_t \frac {U_\infty ^2}{2} \frac {\partial {\overline {C^0_x}}(r)}{\partial r} \sin (\pm \theta + \omega t) \end{equation}

and

(2.25) \begin{equation} f^c_\theta = \frac {U_\infty ^2 }{2} \frac {\partial {\overline {C^0_\theta }}(r)}{\partial r} \left ( \delta y(t) \cos {\theta } + \delta z(t) \sin {\theta } \right) = A_t \frac {U_\infty ^2}{2} \frac {\partial {\overline {C^0_\theta }}(r)}{\partial r} \sin (\pm \theta + \omega t). \end{equation}

The sign $\pm$ preceding $\theta$ determines the azimuthal direction of the circular in-plane translation. The radial forcing is neglected.

The spatial modes of the forcing are shown in figure 4(d–f). As seen, the motion generates an antisymmetric streamwise forcing concentrated near the rotor’s edge and centre, corresponding to the location where the radial gradient of the thrust is highest. The tangential forcing appears mainly near the wake centre.

2.3.3. Rotor precession

Precession refers to the circular motion of the turbine rotor’s axis, which traces out a cone-like path, as shown in figure 3(c). This movement results in a misalignment of the rotor axis out of the wind direction, generating a side force to the wake. The normal direction of the rotor is defined as

(2.26) \begin{align} n_x & = -\sqrt {1 - A_\gamma ^2} , \\[-12pt] \nonumber \end{align}

(2.27) \begin{align} n_y(t) & = A_\gamma \sin (\omega t), \\[-12pt] \nonumber \end{align}

(2.28) \begin{align} n_z(t) & = \pm A_\gamma \cos (\omega t), \end{align}

where $A_\gamma$ is the misalignment between the rotor axis and the oncoming wind. The above expression indicates that the precession can be considered as a combined horizontal and vertical side-to-side oscillatory rotor inclination. The precession direction and the wake helix direction are both determined by the sign of $n_z (t)$ in (2.28).

Analogously to the circular translation, the forcing generated by the circular rotation can also be obtained straightforwardly by the motion-to-forcing model (Li & Yang Reference Li and Yang2024) as combine yaw and tilt motions. By incorporating (2.26)–(2.28) into the motion-to-forcing model, the radial and tangential components of the rotor forcing are expressed as follows:

(2.29) \begin{align} f^c_r & = -\frac {A_\gamma }{2} U_\infty ^2 {\overline {C^0_x}}(r) \sin \left ( \pm \theta - \omega t \right), \\[-12pt] \nonumber \end{align}

(2.30) \begin{align} f^c_\theta & = -\frac {A_\gamma }{2} U_\infty ^2 {\overline {C^0_x}}(r) \cos \left ( \pm \theta - \omega t \right). \end{align}

Here, $A_\gamma$ represents the angle of the rotor’s misalignment. The sign preceding $\theta$ is dependent on the movement azimuthal direction. For this movement, the streamwise forcing is one order smaller than the radial and tangential components (Li & Yang Reference Li and Yang2024), and thus it is neglected here. The spatial modes are shown in figure 4(g–i).

2.4. Resolvent model

The resolvent model is employed to predict the coherent wake fluctuation induced by AWC. As shown in figure 1, the model requires two inputs: (i) a rotor forcing representing the AWC-induced perturbation, computed by the rotor-forcing model $\S$ 2.3; and (ii) a time-averaged baseflow $\overline {\boldsymbol{u}}$ . The model outputs the spatial modes of coherent wake fluctuation ( $\hat {\boldsymbol{u}}$ ). While the inputs and outputs match the previously proposed forcing-to-wake model exactly, our implementation requires updating the resolvent operator in each iteration using the modified time-averaged flow field. In the following, the key principles of the resolvent model will be presented. For a detailed understanding of the fundamental principles, readers may refer to Li & Yang (Reference Li and Yang2024) for the forcing-to-wake model and Rolandi et al. (Reference Rolandi, Ribeiro, Yeh and Taira2024) for broader context on resolvent analysis.

The resolvent operator is derived from the harmonic linearized Navier–Stokes equations that govern the harmonic fluctuating part of the flow, i.e. (2.7) and (2.8). Writing these equations in matrix form, yields

(2.31) \begin{equation} \left ( -i\omega \begin{bmatrix} \boldsymbol{I} & \\ & 0 \end{bmatrix} - \boldsymbol{L}(\overline {\boldsymbol{u}}) \right) \begin{bmatrix} \hat {\boldsymbol{u}} \\ \hat {p} \end{bmatrix} = \begin{bmatrix} \kern3pt \hat {\boldsymbol{\kern-3pt f}} \\ 0 \end{bmatrix}. \end{equation}

Here, $\boldsymbol{I}$ represents a $3 \times 3$ identity matrix arising from the temporal derivatives of $\boldsymbol{u}^c$ , while $\boldsymbol{L}(\overline {\boldsymbol{u}})$ encompasses all remaining linear operations acting on the spatial modes of the fluctuating velocity $\hat {\boldsymbol{u}}$ and pressure $\hat {p}$ in (2.7) and (2.8). The resolvent operator $\boldsymbol{R}(\overline {\boldsymbol{u}})$ is expressed as

(2.32) \begin{equation} \boldsymbol{R}(\overline {\boldsymbol{u}}) = \left ( -i\omega \begin{bmatrix} \boldsymbol{I} & \\ & 0 \end{bmatrix} - \boldsymbol{L}(\overline {\boldsymbol{u}}) \right)^{-1}. \end{equation}

Using this operator, the coherent fluctuations $[\hat {\boldsymbol{u}}, \hat {p}]$ are directly linked to the AWC-induced rotor forcing $\kern3pt \hat {\boldsymbol{\kern-3pt f}}$ , as follows:

(2.33) \begin{equation} \begin{bmatrix} \hat {\boldsymbol{u}} \\ \hat {p} \end{bmatrix} = \boldsymbol{R}(\overline {\boldsymbol{u}}) \begin{bmatrix} \kern3pt \hat {\boldsymbol{\kern-3pt f}} \\ 0 \end{bmatrix}. \end{equation}

To improve computational efficiency, we assume an axisymmetric mean wake and apply Fourier transform in the azimuthal direction, reducing the problem to the two inhomogeneous directions $(x,r)$ for spatial discretization,

(2.34) \begin{equation} \begin{bmatrix} \boldsymbol{u}^c (x,r,\theta,t) \\ p^c (x,r,\theta,t) \\ \boldsymbol{f}^c (x,r,\theta,t) \end{bmatrix} = \frac {1}{2} \begin{bmatrix} \hat {\boldsymbol{u}}(x,r) \\ \hat {p}(x,r) \\ \kern3pt \hat {\boldsymbol{\kern-3pt f}}(x,r) \end{bmatrix} \text{e}^{\mathrm{i} (m \theta -\omega t)} + \frac {1}{2} \begin{bmatrix} \hat {\boldsymbol{u}}^H(x,r) \\ \hat {p}^H(x,r) \\ \kern3pt \hat {\boldsymbol{\kern-3pt f}}^H(x,r) \end{bmatrix} \text{e}^{-\mathrm{i} (m \theta - \omega t)}, \end{equation}

where $m$ is the azimuthal wavenumber. Substituting this ansatz into the linearized Navier–Stokes equations (2.7) and (2.8) yields

(2.35) \begin{align} -\textrm {i}\omega \hat {u}_x & = -\Omega \hat {u}_x -\hat {u}_r \frac {\partial \overline {u}_x}{\partial r}-\hat {u}_x \frac {\partial \overline {u}_x}{\partial x} - \frac {\partial \hat {p}}{\partial x}+ \nu _{\textit {eff}} (\unicode{x1D6E5} + \frac {1}{r^2}) \hat {u}_x +\hat {f}_x, \end{align}

(2.36) \begin{align} -\textrm {i}\omega \hat {u}_r & = -\Omega \hat {u}_r -\hat {u}_r \frac {\partial \overline {u}_r}{\partial r}-\hat {u}_x \frac {\partial \overline {u}_r}{\partial x} + \frac {2\overline {u}_\theta }{r}\hat {u}_\theta -\frac {\partial \hat {p}}{\partial r}+ \nu _{\textit {eff}} \left [\Delta \hat {u}_r-\frac {2\textrm {i}m}{r^2} \hat {u}_\theta \right] + \hat {f}_r, \end{align}

(2.37) \begin{align} -\textrm {i}\omega \hat {u}_\theta & = - \Omega \hat {u}_\theta - \hat {u}_r \frac {\partial \overline {u}_\theta }{\partial r} - \frac {\overline {u}_r}{r} \hat {u}_\theta -\frac {\overline {u}_\theta }{r} \hat {u}_r - \frac { \textrm {i}m}{r} \hat {p} + \nu _{\textit {eff}} \left [\Delta \hat {u}_\theta +\frac {2\textrm {i}m}{r^2} \hat {u}_r\right] + \hat {f}_\theta, \end{align}

(2.38) \begin{align} 0 & = \frac {1}{r} \frac {\partial (r \hat {u}_r)}{\partial r}+\frac {\textrm {i}m}{r} \hat {u}_\theta +\frac {\partial \hat {u}_x}{\partial x}, \end{align}

with the following shorthand operators:

(2.39) \begin{align} \Omega & =\overline {u}_r \frac {\partial }{\partial r}+\frac {\overline {u}_\theta }{r} \textrm {i} m +\overline {u}_x \frac {\partial }{\partial x}, \end{align}

(2.40) \begin{align} \unicode{x1D6E5} & =\frac {1}{r} \frac {\partial }{\partial r}\left (r \frac {\partial }{\partial r}\right)-\frac {m^2+1}{r^2}+\frac {\partial ^2}{\partial x^2}. \end{align}

We utilize the finite difference method to discretize the governing equations. A second-order central differencing scheme is applied to compute all spatial derivatives of the variables in (2.35)–(2.38). Numerical tests demonstrate that this scheme is much more efficient than the previously employed Chebyshev collocation method (Gupta & Wan Reference Gupta and Wan2019; Li & Yang Reference Li and Yang2024), while maintaining the accuracy for wind turbine wake problems.

Figure 5.

The DIPC-induced coherent structures predicted by the linear resolvent model. Helix direction: corotating.

Figure 6.

Same as figure 5, shown at the hub-height for all three velocity components: streamwise (a–c), lateral (d–f) and vertical (g–i) components, respectively.

We employ the resolvent model to predict the linear response of the wake induced by DIPC. To this end, the resolvent operator is constructed from the time-averaged wake field of the baseline case ( $\overline {\boldsymbol{u^0}}$ ) and the rotor forcing model is employed to compute the perturbation source $\kern3pt \hat {\boldsymbol{\kern-3pt f}}$ . Note that this is a linear prediction aiming to demonstrate qualitatively correct behaviour rather than precise quantitative results. Results including mean flow modification effects will be presented and validated with LES in § 3.3.

Figure 5 illustrates three-dimensional coherent structures at different frequencies ( $St = 0.10$ , $St=0.25$ , $St=0.50$ ) with a corotating helix. The visualization shows isosurfaces of streamwise velocity fluctuation $\hat {u}_x = \pm\, 0.05 U_\infty$ in red and green, while the streamwise mean flow velocity is depicted by light blue shadow. The resolvent model, combined with the rotor-forcing model, successfully captures the helical structures and frequency-dependent wake responses. At low frequency ( $St = 0.10$ ), a relatively weak double-helix structure emerges in the wake shear layer, extending from $x = 2D$ to the far wake. The wake response intensifies significantly at $St = 0.25$ , characterized by earlier onset of fluctuation in the near wake and larger helical structures in the far wake, albeit with reduced wavelength. At high frequency ( $St = 0.50$ ), the coherent structures exhibit further reduced wavelength, and the double-helix structure ( $\hat {u}_x = \pm 0.05 U_\infty$ ) terminates before $x = 8D$ , indicating decreased far-wake fluctuation energy. Additionally, a secondary structure appears behind the hub at all frequencies, but remains confined to the near wake region.

These wake responses are further investigated in the hub-height plane for all three velocity components in figure 6. On this plane, the three-dimensional coherent structures are reflected by antisymmetric distribution of $\hat {u}_x$ and symmetric distribution of $\hat {u}_y$ and $\hat {u}_z$ , indicating that the helical meandering wake motion is correctly captured by the resolvent model. Again, the wake fluctuation magnitude is the weakest for $St = 0.10$ and becomes stronger for $St=0.25$ and $St=0.50$ . The reduction of streamwise wavelength is also observed. While the integrated rotor-forcing and linear resolvent models demonstrate qualitatively accurate wake response predictions, subsequent LES comparisons in § 3.3 reveal amplitude overestimation. This discrepancy, also observed in motion-to-wake modelling studies (Li & Yang Reference Li and Yang2024), arises from unaccounted mean wake modifications. To address this limitation, later sections will introduce iterative mean flow correction to improve the predictive accuracy of both mean flow and coherent structures.

2.5. Reynolds stress model

Reynolds stress plays a pivotal role in the momentum exchange between the low-speed wake and the high-speed free stream (Gambuzza & Ganapathisubramani Reference Gambuzza and Ganapathisubramani2023; Van der Hoek et al. Reference van der Hoek, den Abbeele, Simao Ferreira and van Wingerden2024), serving as the crucial link between wake fluctuation and mean recovery. Traditional RS modelling based on mean flow gradients proves insufficient for free shear flows containing large-scale coherent structures (Wu Reference Wu2019), because these structures, significantly contributing to RS, are strongly influenced by upstream disturbances and cannot be effectively described by local mean flow strain rate and effective viscosity alone. This additional contribution from coherent fluctuation manifests mathematically as the additional source term ( $\nabla \cdot \boldsymbol{\tau }^c$ ) in the RANS equations for a mean velocity field (2.9).

Figure 7.

Coherent Reynolds shear stress $\tau _{xy}^c$ on the hub-height plane for DIPC at different frequencies computed using the linear resolvent model’s prediction in figure 6.

In this section, we develop a model to predict the RS generated by both coherent wake fluctuation and small-scale turbulence, where the total RS is decomposed into two components:

(2.41) \begin{equation} \boldsymbol{\tau } = \boldsymbol{\tau ^{c}} + \boldsymbol{\tau ^{s}}. \end{equation}

Here, $\boldsymbol{\tau ^{c}}$ represents the contribution from coherent wake fluctuations and $\boldsymbol{\tau ^{s}}$ accounts for incoherent small-scale turbulence. The coherent component $\boldsymbol{\tau ^{c}}$ is directly computed from the resolvent-predicted wake response $\hat {\boldsymbol{u}}$ through the time averaging of velocity fluctuation products,

(2.42) \begin{equation} \begin{aligned} \boldsymbol{\tau } ^{c}_{ij} & = - \overline {u^c_i u^c_j} \\ & = -\overline { \frac {1}{2}\left (\hat {u}_i e^{-\text{i}\omega t} + \hat {u}_i^H e^{\text{i}\omega t} \right) \cdot \frac {1}{2}\left (\hat {u}_j e^{-\text{i}\omega t} + \hat {u}_j^H e^{\text{i}\omega t} \right)} \\ & = -\frac {1}{4} \left (\hat {u}^{H}_i \hat {u}_j + \hat {u}^{H}_j \hat {u}_i \right) \\ & = -\frac {1}{2} \mathcal{R} \left ( \hat {u}^{H}_i \hat {u}_j \right), \end{aligned} \end{equation}

where $\mathcal{R}$ denotes the real part of a complex number and $\{\cdot \}^H$ denotes the complex conjugate operator. Figure 7 shows the coherent component of the $\tau _{xy}^c$ on the hub-height plane using coherent structures shown in figure 6. As seen, the $\tau _{xy}^c$ for different Strouhal numbers shows distinct spatial patterns and strength, indicating that the model predicts different momentum flux contributed by the coherent structures at different Strouhal numbers. At $St = 0.10$ , the RS remains relatively weak and spatially confined to the centreline of the near wake as a result of the weak coherent structures shown in figures 6(a) and 6(d). As the Strouhal number increases to $St = 0.25$ , the coherent structures generate pronounced $\tau _{xy}^c$ extending to the far wake. The region of elevated $\tau _{xy}^c$ shifts towards the wake edge rather than remaining at the wake centre, in contrast to the $St=0.10$ case. At $St = 0.50$ , both the magnitude and spatial extent of the $\tau _{xy}^c$ diminish compared with $St = 0.25$ , characterized by a notable decrease beyond the intermediate wake region $(x\gt 5D)$ . This reduction implies an attenuation of coherent velocity fluctuations in the far field at $St=0.50$ , which is also observed in figures 6(c) and 6(f). Again, it is worth noting that the result presented in figure 7 is only a qualitative prediction and may deviate from reality, since no mean wake correction is taken into account.

The second component of $\tau$ , i.e. the contribution from incoherent small-scale turbulence, is modelled using the Boussinesq assumption, which introduces an effective eddy viscosity $\nu _{\textit {eff}}$ :

(2.43) \begin{equation} \tau ^{s}_{ij} = \nu _{\textit {eff}} \frac {\partial \overline {u}_i}{\partial x_j}. \end{equation}

It is then connected to the mean flow (2.9) by assuming $\nu _{\textit {eff}} \nabla ^2 \overline {\boldsymbol{u}} \approx \boldsymbol{\nabla \cdot \tau ^s}$ . To determine $\nu _{\textit {eff}}$ , we utilize data from the uncontrolled baseline wake, where the RS arises solely from small-scale turbulence, with the following equation:

(2.44) \begin{equation} - \overline {u_x^{\prime} u_y^{\prime}} = \nu _{\textit {eff}} \frac {\partial \overline {u_x}}{\partial y}. \end{equation}

Here, $-\overline {u_x^{\prime} u_y^{\prime}}$ denotes the hub-height in-plane primary shear stress and ${\partial \overline {u_x}}/{\partial y}$ represents the mean flow strain rate, both obtained from the baseline wake. The effective viscosity $\nu _{\textit {eff}}$ is then computed using a least-squares method at each downstream location. For prediction purposes, this effective viscosity, once determined, remains fixed and is applied to all controlled wake cases. The validity of this assumption will be assessed using large-eddy simulation results in $\S$ 3.2.

2.6. Mean flow model

The mean flow model proposed here predicts the enhanced wake recovery induced by AWC, serving as a critical component in the nonlinear iterative framework (figure 1). Although (2.9) and (2.10) allow direct solution of the time-averaged wake from RS, we adopt a more efficient approach by leveraging the known uncontrolled baseline wake. Specifically, using the RS predictions from § 2.5 as input, the model computes a correction to the baseline flow. This correction-based formulation expresses the modified time-averaged wake as

(2.45) \begin{align} \overline {\boldsymbol{u}} & = {\overline {\boldsymbol{u}^0}} + {\delta \boldsymbol{u}}, \end{align}

(2.46) \begin{align} \overline {p} & = \overline {p^0} + {\delta p}, \end{align}

(2.47) \begin{align} \overline {\boldsymbol{f}} & = {\overline {\boldsymbol{f}^0}} + {\delta \boldsymbol{f}} \approx \overline {\boldsymbol{f}^0}, \end{align}

where $[\overline {\boldsymbol{u}}, \overline {p}, \overline {\boldsymbol{f}}]$ represent the time-averaged velocity, pressure and rotor forcing of the corrected case, $\{\cdot \}^0$ denotes the result of baseline wake and $\delta \{\cdot \}$ denotes the modifications. We assume the time-averaged rotor forcing modification is negligible ( $\delta \boldsymbol{f} \approx 0$ ), which is consistent with the rotor-forcing model developed based on linear assumption.

Assuming that the uncontrolled wake satisfies the following governing equations:

(2.48) \begin{align} \big ({\overline {\boldsymbol{u^{0}}}} \cdot \nabla \big) \overline {\boldsymbol{u^{0}}} + \nabla \overline {p^0} - \nu _{\textit {eff}} \nabla ^2 \overline {\boldsymbol{u^{0}}} & = \overline {\boldsymbol{f^{0}}}, \end{align}

(2.49) \begin{align} \nabla \cdot \overline {\boldsymbol{u^{0}}} & = 0. \end{align}

Subtracting them from the mean flow governing (2.9) and (2.10) results in the governing equations for the correction variables,

(2.50) \begin{align} \big (\overline {\boldsymbol{u^0}} \cdot \nabla \big) \delta \boldsymbol{u} + \left ({\delta \boldsymbol{u}} \cdot \nabla \right) \overline {\boldsymbol{u^0}} + \left ({\delta \boldsymbol{u}} \cdot \nabla \right) \delta \boldsymbol{u} + \nabla ({\delta p}) - \nu _{\textit {eff}} \nabla ^2 \delta \boldsymbol{u}& = \nabla \cdot \boldsymbol{\tau ^c}, \end{align}

(2.51) \begin{align} \nabla \cdot \delta \boldsymbol{u} & = 0, \end{align}

with $\delta \boldsymbol{u}$ and $\delta p$ the unknowns to be solved and $\nabla \cdot \boldsymbol{\tau ^c}$ the source term driving the wake modification. Equations (2.50)–(2.51) are also solved in a two-dimensional $(x,r)$ domain with azimuthal uniformity, given the statistically axisymmetric nature of the forcing and responses. The equations are discretized using a second-order central differencing scheme. The pressure-velocity coupling is solved via the projection method (Chorin Reference Chorin1967). The nonlinear term $ ({\delta \boldsymbol{u}} \cdot \nabla) \delta \boldsymbol{u}$ is treated through inner iterations.

2.7. Iterative solution for self-consistent prediction

The submodels introduced above are integrated and solved through an iterative procedure, as illustrated in figure 1 and detailed in table 1. Beginning with a user-defined AWC and an uncorrected baseline wake, the model sequentially computes the rotor forcing ( $\kern3pt \hat {\boldsymbol{\kern-3pt f}}$ ), the coherent structures ( $\hat {\boldsymbol{u}}$ ), the RS ( $\boldsymbol{\tau }$ ) and the mean flow correction ( $\delta \boldsymbol{u}$ ). After obtaining $\delta \boldsymbol{u}$ , the time-averaged wake $\overline {\boldsymbol{u}}$ is updated through (2.45), serving as the baseflow for the subsequent iteration’s resolvent analysis.

Table 1.

The iterative procedure for computing self-consistent wake response under AWC.

The solution naturally converges through an inherent self-correcting mechanism. For instance, an overestimated velocity fluctuation $\hat {\boldsymbol{u}}$ generates a larger $\delta \boldsymbol{u}$ , which weakens the wake deficit in the subsequent iteration, thereby reducing $\hat {\boldsymbol{u}}$ ; the converse also holds true. To ensure numerical stability, the AWC amplitude is gradually increased from an initial value $A_0$ to the target amplitude $A_f$ over $N = 5$ iterations, after which it remains constant at $A_f$ . Convergence is monitored through the relative variation of mean flow modification $\delta \boldsymbol{u}$ between consecutive steps,

(2.52) \begin{equation} \frac {|\delta \boldsymbol{u}_n - \delta \boldsymbol{u}_{n-1} |}{|\delta \boldsymbol{u}_n|} \lt \epsilon, \end{equation}

where the tolerance $\epsilon = 10^{-2}$ is selected to balance between computational efficiency and solution accuracy, which is discussed in Appendix A.1.

Figure 8.

Convergence characteristics and wake modification predicted by the proposed model: (a) convergence history showing the relative step difference of the mean flow modification; (b–d) hub-height contours of time-averaged streamwise velocity at iterations $n=0$ (initial state), $n=5$ (intermediate) and $n=13$ (converged solution), respectively; (e) streamwise evolution of time- and rotor-area-averaged streamwise velocity, comparing present predictions with large-eddy simulation results. The AWC approach is DIPC with $St=0.25$ and corotating helix.

The nonlinear solution procedure and convergence characteristics are demonstrated in figure 8 for the corotating DIPC strategy and forcing frequency $St=0.25$ . The convergence history, monitored through the relative step difference defined by (2.52), exhibits nearly monotonic decay from unity to below $10^{-2}$ over 14 iterations (figure 8 a). This robust convergence is interrupted only once at $n=5$ , corresponding to the transition from the initial to the target forcing amplitude $A_f$ . Beyond this point, the relative step difference decreases exponentially until reaching the threshold $\epsilon = 10^{-2}$ at iteration $n = 13$ .

The evolution of the hub-height streamwise velocity field $\overline {u}_x$ is illustrated in figure 8(b–d). Starting from the unmodified initial state ( $n=0$ ), the wake undergoes substantial expansion by iteration $n=5$ , with subsequent iterations refining the solution until convergence at $n=13$ . This progression demonstrates the model’s capability to capture the AWC-induced wake recovery through the nonlinear iterative procedure.

Figure 8(e) plots the streamwise evolution of time- and rotor-area-averaged streamwise velocity ( $\widetilde {U}_x$ ) and provides quantitative validation against LES results obtained for the same scenario. In the near-wake ( $0 \lt x \lt 2D$ ), all curves collapse, indicating a decrease of $\widetilde {U}_x$ to $0.4U_\infty$ that matches well with LES predictions. The intermediate wake region ( $2D \lt x \lt 6D$ ) shows the onset of improved wake recovery by LES, which is successfully predicted by the nonlinear solutions at $n=5$ and $n=13$ . The amount of velocity increase reaches approximately $0.2U_\infty$ at $x=6D$ . In the far-wake ( $x \gt 6 D$ ), this velocity difference between the controlled and uncontrolled cases ceases to grow. At $x = 10D$ , the controlled wake reaches approximately 90 % of free stream wind speed. The converged solution ( $n=13$ ) shows the optimal agreement with LES, while the intermediate solution ( $n=5$ ) slightly overpredicts wake velocity. Overall, the good agreement between the present method and LES results suggests that the proposed model successfully captures the essential physics of wake development.

3.1. The LES method and configuration

We simulate the fully nonlinear evolution of wind turbine wakes with the LES code VFS-Wind (Yang et al. Reference Yang, Sotiropoulos, Conzemius, Wachtler and Strong2015), which has been validated and applied to study the fluid mechanisms for wind turbine wake evolution (Yang et al. Reference Yang, Hong, Barone and Sotiropoulos2016; Foti, Yang & Sotiropoulos Reference Foti, Yang and Sotiropoulos2018; Li & Yang Reference Li and Yang2021; Dong et al. Reference Dong, Qin, Li and Yang2023; Li & Yang Reference Li and Yang2024).

The code implements the filtered incompressible Navier–Stokes equations to solve the three-dimensional turbulent airflow. The Coriolis force is not considered to simplify the problem. The dynamic Smagorinsky model (Smagorinsky Reference Smagorinsky1963; Germano et al. Reference Germano, Piomelli, Moin and Cabot1991) is employed for subgrid-scale stress modelling. A well-validated actuator surface method (Yang & Sotiropoulos Reference Yang and Sotiropoulos2018) is employed to model the rotor and nacelle, which are both represented by equivalent aerodynamic forces. The smoothed discrete delta function (Yang et al. Reference Yang, Zhang, Li and He2009) is employed to transfer the aerodynamic forces to the computational grid to avoid singularity issues. The present work employs a structured Cartesian grid for discretizing the computational domain, a second-order central differencing scheme for spatial discretization and a second-order fractional step method for temporal integration. A detailed description of the flow solver can be found in Ge & Sotiropoulos (Reference Ge and Sotiropoulos2007).

3.1.2. Parameter space of LES

Table 2 presents the parameter space for the LES. In total 25 LES are carried out, including one reference baseline case without AWC and 24 cases with three AWC techniques, four Strouhal numbers and two helix directions. The control amplitude is kept small to avoid disturbing wind turbine operation (Wei & Dabiri Reference Wei and Dabiri2023), yet sufficient to generate noticeable wake response and recovery within sensitive frequency ranges (Cheung et al. Reference Cheung, Brown, Houck and Develder2024). The translation and precession amplitudes are equivalent in terms of blade tip displacement: $A_\gamma R = 0.01D = A_t$ . The selected Strouhal numbers span the wake’s characteristic frequency range. Among them, $St=0.25$ and $St=0.35$ are close to the optimal forcing frequency for DIPC, according to the LES of Frederik et al. (Reference Frederik, Doekemeijer, Mulders and van Wingerden2020). The other two frequencies, $St = 0.10$ and $St=0.50$ , are below and above the wake’s most sensitive frequency, as determined by linear resolvent analysis in a previous study on the same wind turbine’s wake (Li & Yang Reference Li and Yang2024).

Table 2.

Parameter space of LES investigation.

3.1.3. Numerical configurations

The computational domain is rectangular, as depicted in figure 9, extending $14D$ in the streamwise direction ( $x$ ) and $7D$ in both transverse ( $y$ ) and vertical ( $z$ ) directions. The wind turbine is positioned on the domain centreline, $3.5D$ downstream of the inlet boundary, with the hub centre defining the origin of the coordinate system. The slip boundary condition is applied to the four lateral boundaries. At the inlet, a uniform inflow of $U_\infty = 10\,\text{ms}^{-1}$ is imposed, neglecting the effects of wind shear, veer, turbulence and the ground. At the outlet, the Neumann boundary condition is applied to the velocity components. The grid interval is $\Delta x = D/20$ and $\Delta y = \Delta z = D/40$ in the wake region and the time step is controlled so that the blade does not pass through more than one grid per step.

Figure 9.

Schematic of the computational domain and the grid for LES.

3.2. The LES results

This section analyses the impact of AWC on coherent structures and their role in mean wake recovery by analysing LES data, to obtain physical insight and to verify important assumptions of the wake modelling.

We begin by examining the instantaneous helical wakes produced by the three different rotor actuation approaches in § 3.2.1, followed by an exploration of the coherent structures in § 3.2.2 and their effect on enhanced wake recovery in § 3.2.3. These analyses highlight the similarities and differences caused by various AWC configurations. Special emphasis is placed on the RS in § 3.2.4, which is an essential element to relate AWC-induced coherent structures and recovery for achieving self-consistent wake prediction in the proposed model. Only the results at $St = 0.10$ , 0.25 and 0.50 are presented in detail, because the results at $St = 0.35$ closely resemble those of $St = 0.25$ and are reserved only for validating the proposed wake model in § 3.3 to avoid redundancy.

3.2.1. Instantaneous wake

Figure 2 provides a three-dimensional visualization of instantaneous wind turbine wakes under various AWC strategies at a forcing frequency of $St=0.25$ . Despite the differences in AWC methods, the wake structures demonstrate notable similarity across cases, characterized by the downstream amplification of perturbations due to convective instability (Gupta & Wan Reference Gupta and Wan2019). Near the rotor, the wake closely resembles the baseline case; however, farther downstream, the helical structure amplifies significantly, causing the wake to deflect from the centreline and form a pronounced meandering pattern. These characteristics align with those observed in DIPC-controlled wakes in earlier studies (Frederik et al. Reference Frederik, Doekemeijer, Mulders and van Wingerden2020; Korb et al. Reference Korb, Asmuth and Ivanell2023), regardless of the specific AWC techniques applied.

Figure 10.

Instantaneous wind turbine wake forced by DIPC (a,b,c), translation (d,e,f) and precession (g,h,i) at different frequencies, St = 0.10 (a,d,g); St = 0.25 (b,e,h), St = 0.50 (c,f,i). Helix direction: corotating.

The effect of excitation frequency on the wake structures under different AWC techniques is illustrated in figure 10. The instantaneous wakes exhibit consistent frequency-dependent trends across all AWC approaches. At the lowest frequency ( $St=0.10$ ), the impact of AWC is minimal, with the wake showing little deformation (figure 10 a,d,g). As the frequency increases to $St=0.25$ , significant amplification of the wake response occurs, producing pronounced large-scale helical deformations (figure 10 b,e,h). However, further increasing the frequency to $St=0.50$ suppresses these large-scale deformations in the far wake, leading instead to the emergence of small-scale wavy deformations (figure 10 c,f,i).

These observations indicate that while initial disturbances differ among AWC techniques, the frequency-dependent behaviour of downstream wake evolution remains qualitatively similar. The reason is that the sensitivity of the wake to specific excitation frequencies is governed by the wake’s intrinsic properties, unaffected by the choice of AWC. The primary effect of different AWCs is to induce different rotor forcing so that the magnitude of the wake response can vary quantitatively.

3.2.2. Coherent structures

To understand in depth the difference made by different AWC, the coherent structures in the wake are analysed. We employ the fast Fourier transform to transform the velocity field from the time domain to the frequency domain. Then, the spatial mode at the forcing frequency on the hub-height horizontal plane ( $xOy$ ) is extracted, since previous studies have shown that the wake dynamics is locked to the excitation frequency (Messmer et al. Reference Messmer, Hölling and Peinke2024). This method facilitates the identification of distinct spatial patterns associated with various AWC techniques and excitation frequencies, enabling detailed examination of the similarities and differences among AWC approaches.

Figure 11.

Coherent structures derived from Fourier transform for various AWC approaches at different excitation frequencies. Streamwise velocity (a), horizontal velocity (b) and vertical velocity components (c). Subpanels (i–iii), (iv–vi) and (vii–ix) represent the same AWC approach, while subpanels (i,iv,vii), (ii,v,viii) and (iii,vi,ix) correspond to the same forcing frequency. Helix direction: corotating.

In figure 11(a) the structure of the streamwise velocity fluctuation is found to be quasiantisymmetric for all considered AWC approaches. This pattern indicates that acceleration on one side of the wake is always accompanied by deceleration on the other side, reflecting the phenomenon that the entire instantaneous wake meanders laterally (Korb et al. Reference Korb, Asmuth and Ivanell2023). In figures 11(b) and 11(c), both horizontal and vertical velocities show a symmetric pattern extending laterally in the wake region, indicating that the entire wake is deflected in the directions perpendicular to the streamwise direction. For all velocity components, the fluctuation alternates in sign as travelling downstream, demonstrating the wake’s meandering behaviour. The wavelength of the coherent structures is inversely related to the actuation frequency. The unstable nature of the wake (Iungo et al. Reference Iungo, Viola, Camarri, Porté-Agel and Gallaire2013; Gupta & Wan Reference Gupta and Wan2019) is revealed by the downstream growth of velocity fluctuation magnitude. The frequency-dependent amplification observed in figure 10 is clearly illustrated by the magnitude of velocity fluctuations. The streamwise velocity fluctuation is strongest at $St=0.25$ and weaker at $St=0.50$ and $St=0.10$ across all AWC approaches. For the lateral velocity components, the fluctuations at $St=0.50$ are comparable in magnitude to those at $St=0.25$ but are segmented into shorter structures, explaining the smaller deformation seen in the instantaneous wake in figure 10.

Qualitative differences across AWC approaches are observed primarily in the vicinity behind the rotor and the near wake. These differences diminish in the far wake, where the coherent structures exhibit qualitative similarity across all AWC approaches. This similarity in the far wake implies a filter effect of the wake which selectively amplifies certain modes determined by its inherent instability properties. Nevertheless, distinct AWC approaches lead to varying fluctuation amplitudes, indicating quantitative differences in their wake modulation effects. The precession approach is found to trigger the strongest fluctuations at $St = 0.10$ , whereas the translation approach results in slightly stronger fluctuations at $St=0.50$ . At $St = 0.25$ , all three AWC approaches yield comparable results in the far wake. Another notable distinction observed in figure 11 is the systematic phase variation across AWC techniques, which is caused by the specific implementation of each AWC technique, as detailed in § 2.3.

3.2.3. Wake recovery

In this section, we investigate the recovery of the time-averaged wake velocity, which is the leading concern for wake mitigation controls (Meyers et al. Reference Meyers, Bottasso, Dykes, Fleming, Gebraad, Giebel, Göçmen and van Wingerden2022).

Figure 12.

Contour of the time-averaged wind speed for different AWC approaches, forcing frequencies. Helix direction: corotation. The Ref case shows the baseline wake without AWC.

Figure 12 presents time-averaged velocity field contours at the hub-height plane to reveal the influence of AWC approaches on wake recovery. The wake recovery exhibits strong frequency dependence for all three AWC techniques, with maximum recovery at $St = 0.25$ , followed by progressively diminishing recovery effects at $St = 0.50$ and $St = 0.10$ . For cases with strong wake recovery, the AWC fundamentally alters the spatial distribution of the wake velocity, characterized by a broadened wake width and a dual-minima lateral profile. According to Korb et al. (Reference Korb, Asmuth and Ivanell2023), the non-Gaussian profile can be explained by the helical deflection of the instantaneous wake, which makes the peak of the probability density function of the instantaneous wake centre moves off the wake centreline. Such dual-peak distributions have also been observed in side-to-side meandering wakes (Li et al. Reference Li, Dong and Yang2022).

Figure 13.

Comparison of time and rotor-area averaged wind speed $\widetilde {U_x}$ for different AWC approaches, forcing frequencies and helix directions.

Figure 13 shows the time- and rotor-area-averaged streamwise velocity ( $\widetilde {U}_x$ ) as a function of downstream distance $(x)$ for various AWC approaches, compared with the baseline reference case without AWC. At $St=0.10$ , all cases show minimal deviation from the reference, indicating negligible wake recovery. For $St=0.25$ and $St=0.50$ , wake recovery is consistently enhanced regardless of AWC approach or helix direction. Specifically, AWC at $St=0.25$ accelerates wind speed more rapidly in the far wake, while $St=0.50$ achieves earlier recovery in the near wake $(x \lt 4D)$ . This behaviour aligns with the strength of velocity fluctuations, which are more pronounced in the far wake at $St=0.25$ but are stronger in the near wake at $St=0.50$ , as shown in figure 11. Additionally, the wake recovery enhancement varies with helix rotation direction, attributable to the non-zero rotational velocity in the baseline wake, as demonstrated by Frederik et al. (Reference Frederik, Doekemeijer, Mulders and van Wingerden2020). Finally, the observed variations among AWC approaches in the far wake reflect the quantitative impact of different rotor actuations on the wake recovery enhancement.

3.2.4. Reynolds stress

In $\S$ 2, we have followed the previous studies of shear flows to reconstruct the RS from contributions of the coherent structures ( $\boldsymbol{\tau ^c}$ ) and the incoherent small-scale turbulence ( $\boldsymbol{\tau ^s}$ ) separately. This section aims to verify this decomposition strategy using LES results presented in figures 14 and 15.

Here, the coherent RS component is approximated by the coherent structures extracted from the LES result using Fourier transform $\boldsymbol{\tau ^c}$ as follows:

(3.1) \begin{equation} \tau ^{c}_{ij,{LES}} = -\frac {1}{2} \mathcal{R} \left (\hat {u}^{H}_{i,{LES}} \times \hat {u}_{j,{LES}} \right), \end{equation}

where $\hat {\{\cdot \}}_{{LES}}$ denotes the Fourier modes at the rotor forcing frequency ( $\omega$ ). The factor ${1}/{2}$ arises from the employed one-sided Fourier transform of a real time-domain signal, which is decomposed into the two-sided complex frequency-domain including both negative and positive frequencies ( $\pm \omega$ ), see (2.6). The incoherent small-scale RS component is approximated using an eddy viscosity and the mean flow strain rate, as follows:

(3.2) \begin{equation} \tau ^s_{ij,{LES}} = \nu _{\textit {eff}} \frac {\partial \overline {u}_{i,{LES}}}{\partial x_j}, \end{equation}

where $\nu _{\textit {eff}}$ is fitted from the baseline wake, according to (2.44). While all components of the RS tensor are computed, we focus on the in-plane Reynolds shear stress due to its dominant role in momentum exchange and recovery of the time-averaged wake (Gambuzza & Ganapathisubramani Reference Gambuzza and Ganapathisubramani2023).

Figure 14.

Contours of in-plane Reynolds shear stress on the hub-height plane for the reference case (a) and cases with corotating DIPC (b-e). Panel (b) shows the total shear stress ( $\tau _{xy}$ ), (c) displays the coherent component ( $\tau _{xy}^c$ ), (d) shows the incoherent component ( $\tau _{xy}^s$ ), and (e) presents the sum of both components ( $\tau _{xy}^c + \tau _{xy}^s$ ).

Figure 15.

Same data with figure 14 at $x=6D$ downstream of the rotor.

In figure 14, we employ the in-plane Reynolds shear stress at the hub-height $(\tau _{xy})$ obtained from LES to analyse the contribution from coherent and incoherent small-scale parts. Figure 14(a) depicts $\tau _{xy}$ of the reference baseline case, i.e. $\tau _{xy}$ without AWC-induced coherent contribution. Figure 14(b–e) show the $\tau _{xy}$ subject to AWC. Here, DIPC is employed with corotating helix. The results plotted are the total $\tau _{xy}$ obtained directly from LES ( $\tau _{xy} = -\overline {u^{\prime }_x u^{\prime }_y}$ ), the coherent component ( $\tau _{xy}^c$ ), the incoherent component ( $\tau _{xy}^s$ ) and the summation of both ( $\tau _{xy}^c+\tau _{xy}^s$ ). At $St=0.10$ , the total $\tau _{xy}$ exhibits patterns remarkably similar to the reference case. The coherent component ( $\tau _{xy}^c$ ) remains negligible throughout the wake region. In contrast, the incoherent component ( $\tau _{xy}^s$ ) dominates their summation ( $\tau ^c_{xy} + \tau ^s_{xy}$ ). This behaviour can be attributed to the weak coherent structures induced by the low-frequency forcing, as evidenced in figure 11. Additionally, the minimal modification of the mean wake at this frequency (figure 12) also supports the strong similarity between the $\tau _{xy}^s$ and $\tau _{xy}$ of the reference case. At $St=0.25$ , $\tau _{xy}$ increases in magnitude and expands in the far wake. Here $\tau _{xy}^c$ is greatly amplified and becomes the dominant part, due to the large-scale helical wake meandering. Conversely, $\tau _{xy}^s$ diminishes in magnitude, as the wake recovery reduces the mean velocity gradient. The sum of both components also closely approximates the directly computed $\tau _{xy}$ . At $St=0.50$ , both extent and magnitude of $\tau _{xy}$ are notably reduced compared with $St=0.25$ . While $\tau _{xy}^c$ remains the primary contributor, its dominance diminishes. Here $\tau _{xy}^s$ shows a slight increase, attributed to the thinner wake shear layer and the increased velocity gradient. The superposition of both components also effectively approximates the directly computed $\tau _{xy}$ .

To further confirm this decomposition strategy, figure 15 presents a quantitative comparison of the same data $6D$ downstream of the rotor. At $St=0.10$ , the $\tau _{xy}$ is contributed almost exclusively by its incoherent component ( $\tau _{xy}^s$ ). At $St=0.25$ , the data demonstrates the significant enhancement of $\tau _{xy}$ due to increased $\tau _{xy}^c$ . At $St=0.50$ , the overall magnitude of $\tau _{xy}$ decreases compared with $St=0.25$ and the relative contributions from Importantly, across the three characteristic scenarios, $\tau _{xy}^s + \tau _{xy}^c$ closely approximates the $\tau _{xy}$ directly computed from LES, and thus confirms the correctness of the present two-component RS modelling approach ( $\S$ 2.5). This validation extends consistently across all investigated AWC approaches, forcing frequencies and helix directions, though additional results are omitted for brevity.

3.3. Model validation

This section validates the proposed self-consistent model through comparison with LES results for both coherent fluctuation and mean wake characteristics. Comparison with the linear resolvent method is also conducted to demonstrate the improvement achieved by the proposed self-consistent model. Given the similarity among different control strategies, the flow field comparisons focus on the DIPC-forced wake. The section concludes with a quantitative analysis of AWC-induced wake recovery across various control strategies, forcing frequencies and helix directions.

For both linear and self-consistent models, a two-dimensional domain ( $xOr$ ) extending $x \in [-2D,10.5D]$ and $r \in [0,3D]$ is employed, similar to that of LES. The grid interval is uniform in the streamwise direction with $\Delta x = D/20$ . In the radial direction, the mesh is uniform in the wake with $\Delta r = D/40$ , and is gradually stretched out for $r\gt 1.2D$ , resulting in $N_x \times N_r$ = $251 \times 71$ grid points. The mesh convergence is studied in Appendix A.2.

Figure 16.

Coherent structures at hub height plane predicted by: (a,d,g) the linear resolvent model; (b,e,h) the self-consistent model; (c,f,i) LES; (a,b,c) the streamwise velocity, (d,e,f) the spanwise velocity and (g,h,i) the vertical velocity. The AWC approach: DIPC at $St=0.25$ with corotate helix.

Figure 17.

Radial profiles of fluctuating velocity magnitudes of the same coherent structures shown in figure 16 at downstream locations $x/D \in \{4,5,6,7,8\}$ : (a) streamwise component, (b) radial component, (c) azimuthal component.

Figure 16 presents a comparative analysis of coherent wake structures at the hub-height plane, comparing predictions from the self-consistent model, linear model and LES, for DIPC-forced wake at $St=0.25$ with corotating helix. Figure 17 compares coherent velocity oscillation magnitudes through radial profiles at multiple downstream positions within $4D \leqslant x \leqslant 8D$ . While both models capture the fundamental oscillatory patterns observed in LES, the self-consistent model achieves superior prediction accuracy in both amplitude and spatial distributions across all velocity components. The streamwise velocity components (figure 16 a–c) exhibit alternating positive and negative velocity regions, asymmetrically distributed about the wake centreline. Conversely, both spanwise (figure 16 d–f) and vertical (figure 16 g–i) velocity components display symmetric patterns about the centreline. The linear model exhibits known limitations (Li & Yang Reference Li and Yang2024). As shown in figure 16(a,d,g), the linear model underestimates the wavelength of coherent structures, particularly in the far wake, since it does not account for the nonlinear effects leading to baseflow modification. Additionally, the stronger shear in the unmodified baseflow causes the linear model to overestimate velocity fluctuation magnitude, as clearly demonstrated in figure 17. The self-consistent model effectively addresses these limitations through several key improvements. Figure 16(b,e,h) shows improved streamwise wavelength agreement with LES, reducing diamond-shaped deformation in streamwise velocity fluctuations compared with the linear model. Figure 17 demonstrates superior prediction of fluctuating velocity magnitudes across all components, better capturing both qualitative and quantitative wake dynamics. Some discrepancies are observed in the contour plots between the self-consistent model and LES, particularly in the near wake, which may be attributed to the modelling simplification and the effect of nacelle in LES.

Figure 18.

The time-averaged wake predicted by different methods: (a,d,g,j) the linear resolvent model, (b,e,h,k) the self-consistent model and (c,f,i,l) LES; (a,b,c) streamwise velocity, (d,e,f) spanwise velocity and (g,h,i) vertical velocity, (j,k,l) $\tau _{xy}$ . The AWC approach: DIPC at $St=0.25$ with corotate helix.

Figure 19.

Radial profiles of the same time-averaged wake characteristics shown in figure 18 at downstream locations $x/D \in \{4,5,6,7,8\}$ .

Figure 18 presents a comprehensive comparison of the time-averaged wake characteristics at the hub-height. The comparison encompasses the streamwise (figure 18 a–c), spanwise (figure 18 d–f) and vertical (figure 18 g–i) velocity components, along with the in-plane RS distribution (figure 18 j–l). Figure 19 plots the same data using radial profiles at downstream locations $4D \leqslant x \leqslant 8D$ . For the linear model, the velocity field is equal to the reference baseline wake, since the linear model does not include any mean wake modification, the in-plane Reynolds shear stresses ( $\tau _{xy}$ and $\tau _{xr}$ ) for the linear model are postprocessed using (2.41)–(2.43) but not employed to correct the mean wake. As seen, the self-consistent model successfully captures the modification made by AWC on the time-averaged field across all velocity components, whereas the linear model is inherently incapable. For the streamwise velocity, the self-consistent model accurately reproduces the apparent expansion and increased wind speed in the far wake, showing excellent agreement with LES results (figures 18 b,c and 19 a). For the spanwise component, an increased converging velocity from the free stream towards the wake centre is observed around $x=5D$ in both the prediction of the self-consistent model and LES (figures 18 e,f and 19 b), indicating that the enhanced flow entrainment to the wake region is well captured by the proposed model. Moreover, an improved agreement of the out-plane velocity ( $\overline {u}_z$ in Cartesian coordinates and $\overline {u}_\theta$ in the cylindrical coordinates) is found between the self-consistent model and LES at $x\gt 5D$ , indicating that the proposed model captures the increased rotational velocity due to the corotating helix (figures 18 h,i and 19 c). Lastly, $\tau _{xy}$ ( $\tau _{xr}$ ) predicted by the self-consistent model achieves improved agreement with LES in both magnitude and spatial distribution (figures 18 k,l and 19 d), explaining the model’s success in time-averaged wake prediction, since it is the key factor driving wake recovery in the mean momentum (2.9).

Figure 20.

Wake response predicted by the self-consistent model at different frequencies: $St = 0.10$ (a,d,g), 0.25 (b,e,h) and 0.50 (c,f,i). The AWC approach: DIPC with corotating helix.

Figure 21.

Comparison of radial profiles of the total fluctuating velocity magnitude $|\hat {\textbf{u}}| = \sqrt {\hat {u}^2_x + \hat {u}^2_r +\hat {u}^2_\theta }$ predicted by linear model (figure 6), self-consistent model (figure 20) and LES (figure 11) at downstream locations $x/D \in \{4,5,6,7,8\}$ : (a) $St = 0.10$ , (b) $St = 0.25$ , (c) $St = 0.50$ . The AWC approach: DIPC with corotating helix.

Figure 20 presents wake responses predicted by the self-consistent model at $St = 0.10$ (figure 20 a,d,g), $St = 0.25$ (figure 20 b,e,h) and $St = 0.50$ (figure 20 c,f,i). Figure 21 compares the magnitude of the coherent velocity fluctuation predicted by the linear model (figure 6), the self-consistent model (figure 20) and the LES (figure 11) by plotting radial profiles of the fluctuating velocity magnitude ( $|\hat {\textbf{u}}| = \sqrt {\hat {u}^2_x + \hat {u}^2_r +\hat {u}^2_\theta }$ ) at downstream locations in the far wake ( $4D \leqslant x \leqslant 8D$ ). At $St=0.10$ , the difference between linear model, self-consistent model and LES is minimal due to weak AWC-induced wake response and minimal nonlinearity. For the case with $St=0.25$ , where the large-scale meandering is the most prominent, good agreement between the self-consistent model and LES is found. For $St = 0.50$ , although discrepancies between the self-consistent model and LES exist, the overestimation of fluctuating velocities is alleviated compared with the linear model.

Figure 22.

The time- and rotor-area-averaged streamwise velocity predicted by the self-consistent model and LES for all considered AWC approaches: DIPC (a,d), translate (b,e), precession (c,f). The helix direction is corotating and counter-rotating for (a–c) and (d–f), respectively. The uncontrolled baseflow is provided for reference. The results are plotted at $x=6D$ .

Finally, the predictive capabilities on the enhanced wake recovery of the self-consistent model are quantitatively assessed. Here, the self-consistent model is employed to predict wake velocity over a broad range of forcing frequencies ( $0.1 \leqslant St \leqslant 0.6$ ), for all three forcing strategies, and with both helix directions. Figure 22 presents the time- and rotor-area-averaged streamwise velocity ( $\widetilde {U}_x$ ) at $6D$ downstream the rotor, comparing predictions from the linear resolvent model (without mean flow correction), the self-consistent model and LES. Compared with the baseline state, the self-consistent wake model successfully predicts the enhanced wake recovery following the same frequency dependency trend showed by LES results: the highest recovery occurs around $St \approx 0.35$ and diminishes when the frequency moves towards both extremities for all considered AWC approaches. Moreover, quantitative agreement is found between the self-consistent model and LES for all the considered cases, with the largest relative discrepancy less than 8 % (the maximum error occurs for the counter-rotating DIPC case at $St=0.35$ ). The discrepancies between the self-consistent model and LES may be associated with the simplifications in the modelling, and shows the room for further improvements. Even in this least favourable scenario, the velocity prediction accuracy is still remarkably improved against the linear model, which excludes the mean flow modification. For the linear model, the time-averaged wake is always equal to the uncontrolled baseline wake, which shows a 28.0 % underestimation of the mean wake velocity. Overall, the self-consistent model demonstrates a reliable and accurate mean velocity prediction across diverse control strategies, demonstrating its potential as a computationally efficient alternative to LES for AWC designs.