Modeling the effects of active wake mixing on wake behavior through large-scale coherent structures

The use of active wake mixing (AWM) to mitigate downstream turbine wakes has created new opportunities for reducing power losses in wind farms. However, many current analytical or semi-empirical wake models do not capture the flow instabilities that are excited through the blade pitch actuation. In this work, we develop a framework, which accounts for the impacts of the large-scale coherent structures and turbulence on the mean flow, for modeling AWM. The framework uses a triple-decomposition approach for the unsteady flow field and models the mean flow and fine-scale turbulence with a parabolized Reynolds-averaged Navier–Stokes (RANS) system. The wave components are modeled using a simplified spatial linear stability formulation that captures the growth and evolution of the coherent structures. Comparisons with high-fidelity large eddy simulations (LESs) of the turbine wakes showed that this framework was able to capture the additional wake mixing and faster wake recovery in the far-wake regions for both the pulse and helix AWM strategies with minimal computational expense. In the near-wake region, some differences are observed in both the RANS velocity profiles and initial growth of the large-scale structures, which may be due to some simplifying assumptions used in the model.

1 Introduction

Wake propagation downstream of turbines in a wind farm is often dominated by the evolution of large-scale coherent structures. These features may arise from unsteady atmospheric conditions, such as the naturally occurring streaks in an atmospheric boundary layer (Zhang and Stevens, 2020), or result from various control strategies intentionally applied to the turbines upwind. Given the wide range of spatial and temporal scales involved, reduced-order models are required to efficiently capture the evolution of turbine wakes and enable evaluation of the impact and performance of various wake control strategies at the farm scale. Early models for this purpose were based on steady-state flow assumptions and are thus unable to capture critical dynamic aspects of wake evolution. In this paper, we develop an improved reduced-order wake model that incorporates time-dependent propagation physics. We show that this new model can be used to more effectively capture wake perturbation and recovery dynamics. Although we focus the present work on capturing coherent structures generated by the turbine control strategies, our methods are likely to be relevant to other time-dependent sources as well.

Wind farm flow control methods are designed primarily to reduce power losses in wind farms due to the effects of wakes on downstream turbines. Common approaches include static or dynamic adjustments to the induction factor (turbine derating), yaw angle (wake steering), or blade pitch (wake mixing) (Meyers et al., 2022) of upstream turbines. The present work focuses specifically on active wake mixing (AWM), which aims to excite flow instabilities in the wake that enhance the entrainment of mean velocity, thereby accelerating wake recovery.

Following Cheung et al. (2024a), AWM can be implemented by specifying a dynamic blade pitch, θ(t), on top of the baseline pitch set point, θ₀(t), as

$$\begin{array}{}\text{(1)}& \mathit{\theta }\left(t\right)={\mathit{\theta }}_{\mathrm{0}}\left(t\right)+A\cos ({\mathit{\omega }}_{\mathrm{e}}t-n\mathit{\psi }(t)+{\mathit{\psi }}_{\text{clock}}),\end{array}$$

where A is the pitching amplitude, ω_e is the excitation frequency, ψ is the azimuth position of the blade, ψ_clock is the clocking angle, and n is an azimuthal wavenumber. The parameter n controls the structure of the flow instabilities imparted on the wake and is often used to distinguish between different AWM strategies. Examples include the pulse method (n=0), which generates an axisymmetric instability in the flow through collective blade pitching (Goit and Meyers, 2015; Munters and Meyers, 2018), and the helix method ($n=-\mathrm{1}$), which uses individual pitch control to impart a helical structure on the wake that rotates in the direction opposite to the turbine rotor (Frederik et al., 2020a). The instabilities are actuated in the wake according to an excitation frequency, which can be specified as a function of the Strouhal number, St, the inflow velocity, U_inf, and turbine diameter, D, as ${\mathit{\omega }}_{\mathrm{e}}=\mathrm{2}\mathit{\pi }\mathit{St}{U}_{\text{inf}}/D$. Strouhal numbers based on the natural unsteady properties of the wake are typically sought (St≈0.3), leading to flow structures that are generated over much longer periods than a rotor period (Frederik et al., 2020b).

One type of existing reduced-order wake model is the steady-state analytical one. This type of model finds its roots in the Jensen (Jensen, 1983) and Ainslie (Ainslie, 1988) models, for instance, and also includes more recent sophisticated versions such as the cumulative curl model (Bastankhah and Porté-Agel, 2014). These models are available for optimizing wind farm performance through the FLORIS code (Sinner and Fleming, 2024). However, the inherently steady-state nature of their implementations in FLORIS, along with their reliance on empirical tuning, limits their applicability in scenarios where unsteady flow features are critical. Recent work has shown the importance of unsteady flow features for AWM by connecting the performance of different AWM strategies to the underlying fluid mechanics associated with the induced flow instabilities, particularly the interactions between unsteady coherent flow structures and wake recovery dynamics (Korb et al., 2023; Cheung et al., 2024a). Notably, Cheung et al. (2024a) introduced a spatial linear stability analysis to quantify the growth characteristics of initial flow disturbances based on the temporal forcing frequency and forced azimuthal wavenumber and showed a correlation between turbulent entrainment statistics in the wake and modal energy gain. These findings suggest that an accurate model for AWM should be capable of representing the unsteady effects of coherent structures on the flow.

Although wind farm optimization typically relies on steady-state models, several approaches for dynamic wake modeling have also been developed. These methods often rely on data-driven representation of the coherent flow structures. Proper orthogonal decomposition (POD), for instance, is a commonly used data analysis technique for identifying the energetic structures in a flow (Lumley, 1967), and it has been applied to a wide range of applications including characterizing the coherent structures in a wind farm (Bastine et al., 2015; Ali et al., 2017; Hamilton et al., 2018; Zhang and Stevens, 2020). In the context of AWM, Yalla et al. (2025) demonstrated that spectral POD provides a useful representation of the coherent structures induced by dynamic blade pitch actuation, connecting the frequency and wavenumber inputs used by the turbine controller to structures in the wake. Spectral POD is closely related to other data-driven reduced-order modeling techniques such as resolvant analysis and dynamic mode decomposition (DMD) (Towne et al., 2018). Li and Yang (2024) developed a resolvant-based model to represent the wake of floating offshore wind turbines subjected to dynamic platform motions, which produce similar wake responses as the pulse and helix forcing strategies. Gutknecht et al. (2023) developed an AWM DMD model for the wake of a single actuated turbine, which easily scaled to different wind speeds and forcing Strouhal numbers. While these data-driven methods are valuable, they can require a substantial amount of training data to provide accurate representations of coherent flow structures. Generating these data can be prohibitively expensive, and avoiding this expense is a primary motivation for reduced-order modeling. Moreover, this reliance on training data can limit the applicability of these models outside the specific conditions for which they were developed and can make adopting these models challenging. Although other dynamic modeling approaches, such as the dynamic wake meandering model (Larsen et al., 2007; Madsen et al., 2010), exist, they too often rely on external turbulence simulations to provide the dynamic components of the flow. In this work, we propose an analytical representation of the coherent flow structures in the wake and limit the training to model constants, which should enhance the robustness and adaptability of the model.

The behavior of large-scale coherent structures in various canonical shear flows is a well-studied problem with a vast body of existing literature. A number of previous works describe the formation and behavior of these structures in turbulent boundary layers (Hussain, 1986; Robinson, 1991), free shear layers (Ho and Huerre, 1984), jets (Crow and Champagne, 1971), and wakes (Fuchs et al., 1979). Of particular relevance to the current work are the theoretical and modeling approaches used to analyze such flows. Hussain and Reynolds (1970) introduced the concept of a triple-decomposition analysis to separate the mean flow, fine-scale turbulent components and wave components of flow, which was widely used in modeling jet (Iqbal and Thomas, 2007) and boundary layer (Kwon et al., 2016) flows. A number of previous studies have shown that the growth of the coherent structures in shear flows can be modeled by spatial stability theory (Cheung and Lele, 2009; Cheung and Zaki, 2011) and noted that the coupling of the growth of the large-scale structures to the mean flow evolution was critical to capturing the behavior of the flow. However, these modeling approaches have yet to be applied to the problem of turbine wakes, leading to a large gap between the currently available steady-state wake models and computationally expensive, high-fidelity simulations.

The objective of the current work is to develop a physics-based, computationally efficient model that can capture the effects of active wake mixing on turbine wakes. While turbine wakes contain significant differences from the canonical jet flows discussed above, we show that, by using a triple-decomposition approach, we can still capture the mean flow using a parabolized Reynolds-averaged Navier–Stokes (RANS) model, and the large-scale structures can be modeled with a spatial linear stability formulation. In the following sections, we describe the mathematical formulation used in this study and the high-fidelity numerical simulations of the turbine wakes used to calibrate and evaluate the reduced-order model. We then show comparisons between LES calculations and the RANS with a linear stability model for different AWM strategies, and we conclude with a summary of the work and a discussion of future work in this area.

2 Methodology

2.1 Atmospheric and turbine conditions of interest

Though the current model is meant to be broadly applicable to all turbine wake flows from both onshore and offshore horizontal axis wind turbines, this work focuses on modeling AWM as applied to larger offshore wind turbines under low-turbulence, relatively steady atmospheric conditions. In these situations, the application of AWM can potentially lead to substantial wake benefits and noticeable annual energy production (AEP) gains. For offshore locations, the prevalence of these conditions can also lead to many situations where the turbine wakes are especially long and provides an opportunity to improve wind farm power performance.

Representative offshore conditions were selected for this study based on measured data from a floating lidar measurement campaign conducted off the New York bight (Mason, 2022; DNV, 2022). The floating lidar data, collected in 10 min intervals over a period of 1.6 years, provided velocity and turbulence intensity (TI) profile information for heights between 20 and 200 m. From these data, a selection process was undertaken to generate three representative wind speed profiles with relatively low TI (see Table 1). To generate the precursor simulations, small velocity and temperature perturbations were introduced near the surface to accelerate turbulence development. The low-TI conditions were produced by imposing negative ground surface temperature rates and adjusting the surface roughnesses, followed by 10 000 s of flow time. As such, the generated conditions were stable atmospheric boundary layers. More details of the precursor generation process and comparison versus the measured atmospheric data are described in Brown et al. (2025).

Table 1 Hub-height wind speed conditions used in the turbine wake study. All values are taken from the simulated atmospheric boundary layer.

These conditions corresponded to the likely operating range of the IEA 15 MW reference turbine where AWM strategies might be deployed. The IEA 15 MW reference turbine was used in this study due to its similarity to current offshore wind turbines being developed by major turbine original equipment manufacturers (OEMs). The details of this turbine are summarized in Table 2, with additional information available from Gaertner et al. (2020).

Table 2 Details of the IEA 15 MW reference turbine. TSR: tip speed ratio.

As discussed in Sect. 2.5.1, the selected wind conditions and turbine model were used to set up high-fidelity LES calculations for turbines with and without AWM activated. The LES data were then used to calibrate RANS closure model coefficients in the wake model and to compare the accuracy of the final outputs from the model.

2.2 Mathematical formulation

To model both the steady-state wake profiles and the unsteady dynamics of coherent structures that may be excited through AWM, we use the triple-decomposition approach pioneered by Hussain and Reynolds (1970) in their studies of boundary layer flows. The triple-decomposition formulation separates the flow velocity u(x,t) into the three components,

$$\begin{array}{}\text{(2)}& \mathbf{u}(\mathbf{x},t)=\stackrel{\mathrm{‾}}{\mathbf{U}}\left(\mathbf{x}\right)+\stackrel{\mathrm{̃}}{\mathbf{u}}(\mathbf{x},t)+{\mathbf{u}}^{\prime }(\mathbf{x},t),\end{array}$$

where $\stackrel{\mathrm{‾}}{\mathbf{U}}\left(\mathbf{x}\right)$ is the time-averaged mean flow, $\stackrel{\mathrm{̃}}{\mathbf{u}}(\mathbf{x},t)$ is the wave component of the velocity, and ${\mathbf{u}}^{\prime }(\mathbf{x},t)$ represents the fine-scale turbulent fluctuations. The mathematical operations required to compute the mean flow over an averaging time, T, are given by

$$\begin{array}{}\text{(3)}& \stackrel{\mathrm{‾}}{f\left(\mathbf{x}\right)}={\displaystyle \frac{\mathrm{1}}{T}}\underset{\mathrm{0}}{\overset{T}{\int }}f(\mathbf{x},t)\mathrm{d}t,\end{array}$$

and the phase average is defined by

$$\begin{array}{}\text{(4)}& \langle f\rangle ={\displaystyle \frac{\mathrm{1}}{N}}\sum _{m=\mathrm{0}}^{N}f(\mathbf{x},t+m\mathit{\tau }),\end{array}$$

for a given time period, τ, of the coherent structure and for a specified number of periods, N. Once the mean and the phase-averaged velocities are known, the wave component, $\stackrel{\mathrm{̃}}{\mathbf{u}}(x,t)$, of the flow field is defined as

$$\begin{array}{}\text{(5)}& \stackrel{\mathrm{̃}}{\mathbf{u}}(x,t)=\langle \mathbf{u}(x,t)\rangle -\stackrel{\mathrm{‾}}{\mathbf{U}}\left(\mathbf{x}\right),\end{array}$$

and the fine-scale fluctuating components can be calculated as

$$\begin{array}{}\text{(6)}& {\mathbf{u}}^{\prime }(\mathbf{x},t)=\mathbf{u}(\mathbf{x},t)-\stackrel{\mathrm{‾}}{\mathbf{U}}\left(\mathbf{x}\right)-\stackrel{\mathrm{̃}}{\mathbf{u}}(\mathbf{x},t).\end{array}$$

An example of a turbine wake that has been triply decomposed is shown in Fig. 1. In this case, the turbine flow field was calculated using LES and averaged according to Eqs. (3) and (4), which leads to clear depictions of the mean flow field features as well as the large-scale coherent structures that develop within the wake. In the current study, an averaging time of T=600 s is typically applied to the unsteady LES data.

Figure 1 An example of a triply decomposed flow field for a wind turbine wake. This case is from the HelixA4 case in the low-WS ABL conditions with 4° amplitude forcing. In each contour, the normalized streamwise velocity, $U/U_{\mathrm{\infty }}$, is plotted.

One advantage of using the triple-decomposition approach is that it allows computationally efficient models to be developed that can solve for each of the three components. Interactions among the different flow components can also be included, which shows how the large-scale coherent structures can impact the mean flow and vice versa. In the following sections, we describe how a parabolic RANS model can be used to efficiently capture the mean flow and fine-scale turbulent flow components. This is coupled to a linear stability model for the wave components of the flow, and we show that, as the large-scale coherent structures develop within the wake, the mean velocity profiles are impacted as well, leading to the desired wake mixing behavior in this application.

The current model described in the work applies to the wake, immediately downstream of the rotor, of a single turbine. The unsteady inflow effects and the rotor loading dynamics are not explicit in this formulation, and the behavior of more complicated phenomena, such as the merging of multiple wakes, is not considered here. With additional development, we intend to extend the current model to wind farm configurations with multiple turbines, but we focus on the single-turbine wake dynamics initially.

2.3 RANS model

In this section, a RANS model is formulated to couple the effects of coherent structures and turbulence on the evolution of the mean velocity field. The model is formulated here in the context of the standard k–ε RANS closure model (Jones and Launder, 1972), although it may be easily adapted to many turbulence closure modeling approaches. A few assumptions are made to simplify the model, as the focus is on developing a computationally efficient representation of the effects of active wake mixing on the mean flow. First, the dynamics are assumed to be axisymmetric, reducing the complexity of the model to two dimensions. A schematic of the streamwise, radial, and azimuthal coordinates relative to the turbine is shown in Fig. 2. Second, the boundary layer approximation is applied so that (1) second-order derivatives in the streamwise direction, x, are small relative to those in the radial direction, r; (2) the radial pressure gradient is decoupled from the velocity field; and (3) turbulent production is dominated by the radial mean streamwise velocity gradient. This leads to a parabolic system that can be marched in the downstream direction, as in Cheung et al. (2024b). Lastly, direct interactions between the coherent structures and the turbulence are neglected so that AWM only forces the evolution of the mean velocity. The resulting equations are

$$\begin{array}{}\text{(7a)}& {\displaystyle \frac{\partial \overline{U}}{\partial x}}+{\displaystyle \frac{\partial \overline{V}}{\partial r}}+{\displaystyle \frac{\overline{V}}{r}}=\mathrm{0},\text{(7b)}& {\displaystyle }\overline{U}{\displaystyle \frac{\partial \overline{U}}{\partial x}}+\overline{V}{\displaystyle \frac{\partial \overline{U}}{\partial r}}={\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial }{\partial r}}\left[r(\mathit{\nu }+{\mathit{\nu }}_t){\displaystyle \frac{\partial \overline{U}}{\partial r}}\right]+F_{\text{CS}},\text{(7c)}& {\displaystyle }\begin{array}{rl}\overline{U}{\displaystyle \frac{\partial k}{\partial x}}+\overline{V}{\displaystyle \frac{\partial k}{\partial r}}=& {\mathit{\nu }}_t{\left({\displaystyle \frac{\partial \overline{U}}{\partial r}}\right)}^{\mathrm{2}}-\mathit{\epsilon }\\ & +{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial }{\partial r}}\left[r(\mathit{\nu }+{\mathit{\nu }}_t/{\mathit{\sigma }}_k){\displaystyle \frac{\partial k}{\partial r}}\right],\end{array}\text{(7d)}& {\displaystyle }\begin{array}{rl}\overline{U}{\displaystyle \frac{\partial \mathit{\epsilon }}{\partial x}}+\overline{V}{\displaystyle \frac{\partial \mathit{\epsilon }}{\partial r}}=& {\displaystyle \frac{C_{\mathrm{1}\mathit{\epsilon }}\mathit{\epsilon }}{k}}\left[{\mathit{\nu }}_t{\left({\displaystyle \frac{\partial \overline{U}}{\partial r}}\right)}^{\mathrm{2}}\right]-{\displaystyle \frac{C_{\mathrm{2}\mathit{\epsilon }}}{k}}{\mathit{\epsilon }}^{\mathrm{2}}\\ & +{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial }{\partial r}}\left[r(\mathit{\nu }+{\mathit{\nu }}_t/{\mathit{\sigma }}_{\mathit{\epsilon }}){\displaystyle \frac{\partial \mathit{\epsilon }}{\partial r}}\right],\end{array}\end{array}$$

where $\overline{U}$ and $\overline{V}$ are the mean streamwise and radial velocity components, respectively. The effects of turbulence, u^′, on the mean flow are represented by the eddy viscosity, ${\mathit{\nu }}_t=C_{\mathit{\mu }}{k}^{\mathrm{2}}/\mathit{\epsilon }$, where k and ε are the turbulent kinetic energy and dissipation, respectively. To close the k–ε model, the RANS constants, C_1ε, C_2ε, C_μ, σ_k, and σ_ε, are calibrated based on LES data as discussed in Sect. 2.3.1. The term F_CS represents forcing of the mean flow by the wave component, $\stackrel{\mathrm{̃}}{\mathbf{u}}$, as

$$\begin{array}{}\text{(8)}& F_{\text{CS}}=-\stackrel{\mathrm{‾}}{\stackrel{\mathrm{̃}}{u}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{u}}{\partial x}}+\stackrel{\mathrm{̃}}{v}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{u}}{\partial r}}},\end{array}$$

and the coupling between the mean component and wave component is discussed further in Sect. 2.4.

Figure 2 Schematic of the wind turbine coordinate system used in this work. The downstream, radial, and azimuthal coordinates are given by x, r, and ψ, respectively, and the turbine rotor radius is given by R. Note that a negative value of ψ is depicted above.

Equations (7a)–(7d) are discretized on a uniform grid in the radial direction using a second-order centered difference method. The radial domain extends to r_max=5R with a uniform spacing of Δr=0.025R. In the x direction, the equations are discretized around the cell centers and a Crank–Nicolson method is used to march 20D downstream with uniform step sizes of Δx=0.1R. The resulting tridiagonal system is solved using an iterative solver, which advances the solution from one x location to the next.

For each variable, Neumann boundary conditions are applied at r=0,

$$\begin{array}{c}{\displaystyle \frac{\partial \overline{U}}{\partial r}}(r=\mathrm{0})=\mathrm{0},{\displaystyle \frac{\partial k}{\partial r}}(r=\mathrm{0})=\mathrm{0},{\displaystyle \frac{\partial \mathit{\epsilon }}{\partial r}}(r=\mathrm{0}),\\ {\displaystyle \frac{\partial \overline{V}}{\partial r}}(r=\mathrm{0})=\mathrm{0},\end{array}$$

and Dirichlet boundary conditions are applied at r=r_max,

$$\begin{array}{c}{\displaystyle }\stackrel{\mathrm{‾}}{U}(r=r_{\text{max}})=U_{\mathrm{\infty }},k(r=r_{\text{max}})=k_{\mathrm{\infty }},\\ {\displaystyle }\mathit{\epsilon }(r=r_{\text{max}})=\mathrm{0},\end{array}$$

where $k_{\mathrm{\infty }}/U_{\mathrm{\infty }}^{\mathrm{2}}=\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{3}}$ is specified based on the LES calibration data. Note that the continuity relation (Eq. 7a) only requires one boundary condition to be imposed on $\overline{V}$.

A hyperbolic tangent profile is used to model the initial condition for $\stackrel{\mathrm{‾}}{U}$ at the initial streamwise location, x=x₀:

$$\begin{array}{}\text{(9)}& \stackrel{\mathrm{‾}}{U}\left(x_{\mathrm{0}}\right)=\mathrm{0.5}(U_{\mathrm{\infty }}-U_{\mathrm{0}})\left(\mathrm{1}+\tanh \left({\displaystyle \frac{r-r_{\mathrm{e}}}{\mathrm{\Delta }}}\right)\right)+U_{\mathrm{0}},\end{array}$$

where the nondimensionalized values $U_{\mathrm{0}}/U_{\mathrm{\infty }}=\mathrm{0.5}$, $r_{\mathrm{e}}/R=\mathrm{1.2}$, and $\mathrm{\Delta }/R=\mathrm{0.05}$ were determined to provide good agreement with the azimuthally averaged velocity deficit profiles from the LES data near $x/D=\mathrm{2}$ (see Figs. 5 and 6). However, it should be noted that the effects of the nacelle in the near wake are not accounted for in the RANS formulation. The initial profile for k is taken to be proportional to the square of the mean velocity gradient, $k\left(x_{\mathrm{0}}\right)\sim (\partial \stackrel{\mathrm{‾}}{U}/\partial r{)}^{\mathrm{2}}+k_{\mathrm{\infty }}$, such that $\sqrt{\mathrm{3}max\left(k\right(x_{\mathrm{0}}\left)\right)/\mathrm{2}}=\mathrm{0.125}$, and the initial ε is taken to balance turbulent kinetic energy production, i.e., $\mathit{\epsilon }\left(x_{\mathrm{0}}\right)=\sqrt{C_{\mathrm{1}\mathit{\epsilon }}{k}^{\mathrm{2}}(\partial \stackrel{\mathrm{‾}}{U}/\partial r{)}^{\mathrm{2}}}$.

2.3.1 Calibration of RANS

The coefficients C_μ, C_1ε, and C_2ε of the k–ε RANS closure model were calibrated to match the rotor-averaged velocities from the baseline LES discussed in 2.5. Since the RANS formulation does not account for the hub and nacelle region from the LES (see Fig. 1), the calibration was formulated to match the rotor-averaged velocities from a distance of $x/D=\mathrm{2.0}$ to $x/D=\mathrm{8.0}$. The cost function for this calibration was an ℒ₂ norm error between the RANS output and the LES output. The L-BFGS-B (Byrd et al., 1995; Zhu et al., 1997) algorithm as implemented in scipy was used for the calibration. The optimal values from this calibration are C_μ=0.0035, C_1ε=0.163, and C_2ε=2.86. It is important to note that these values are particular to the initial conditions and RANS closure model used and not a general guideline for wake predictions. It is also to be noted that the calibration is only performed for the baseline cases and not the AWM cases. The constants σ_k and σ_ε were not included in the calibration process; instead, the standard values of σ_k=1.0 and σ_ε=1.3 proposed by Jones and Launder (1972) were used. Figure 8 shows a close match between the baseline RANS and LES results, showing that these calibrated parameters are representative model constants for capturing the wake behavior in the baseline cases and will be used for all the RANS results presented in this work.

2.4 Linear stability model

In this work, we are primarily interested in evaluating the feasibility of using a wave component model to determine the impact of large-scale coherent structures on the turbine wakes. Many approaches have been used previously in the literature to capture the dynamics of large-scale structures in shear flows, including linear and nonlinear stability analysis (Cheung and Lele, 2009), non-modal stability analysis (Hack and Zaki, 2015), and global stability analysis (Schmid, 2007). These methodologies have been very well developed and successful in analyzing other canonical flows such as pipe flows, boundary layers, and jets.

As an initial step towards demonstrating the feasibility of this modeling approach, a simple parallel-flow, inviscid, spatial linear stability analysis was chosen for this work. The focus of the analysis is to model the growth of the large-scale coherent structures and capture resulting changes to the mean flow of the turbine wake with minimal computational effort. Additional effects not captured in this analysis, such as the effects of shear, veer, swirl, or atmospheric stratification, will be included in future analyses.

2.4.1 Piecewise constant velocity profile

Analytic solutions to the spatial linear stability problem are possible if we assume the initial wake profile remains axisymmetric and roughly follows a piecewise constant velocity profile. In the current work, we adopt the two-step profile shown in Fig. 3, which is defined by

$$\begin{array}{}\text{(10)}& {\stackrel{\mathrm{‾}}{U}}_{\text{pw}}\left(r\right)=\left\{\begin{array}{ll}{U}_{\mathrm{0}},& r<r_{\mathrm{1}}\\ U_{\text{half}},& r_{\mathrm{1}}\le r\le r_{\mathrm{2}}\\ U_{\mathrm{\infty }},& r>r_{\mathrm{2}},\end{array}\right.\end{array}$$

where U₀ is the centerline velocity, U_∞ is the freestream velocity, and $U_{\text{half}}=\frac{\mathrm{1}}{\mathrm{2}}(U_{\mathrm{0}}+U_{\mathrm{\infty }})$ is the averaged velocity of the wake shear region from $r_{\mathrm{1}}\le r\le r_{\mathrm{2}}$. During the solution process, the r₁ and r₂ parameters can be chosen so that the displacement, δ, and momentum, δ_θ, areas of the ${\stackrel{\mathrm{‾}}{U}}_{\text{pw}}$ profile match the same displacement and momentum areas calculated from the RANS mean flow profiles, ${\stackrel{\mathrm{‾}}{U}}_{\text{RANS}}$. We use the following definitions for δ and δ_θ:

$$\begin{array}{}\text{(11)}& {\displaystyle }\mathit{\delta }\left(\stackrel{\mathrm{‾}}{U}\right)=\mathrm{2}\mathit{\pi }\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}\left(\mathrm{1}-{\displaystyle \frac{\overline{U}\left(r\right)}{U_{\mathrm{\infty }}}}\right)r\mathrm{d}r,\text{(12)}& {\displaystyle }{\mathit{\delta }}_{\mathit{\theta }}\left(\stackrel{\mathrm{‾}}{U}\right)=\mathrm{2}\mathit{\pi }\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\displaystyle \frac{\overline{U}\left(r\right)}{U_{\mathrm{\infty }}}}\left(\mathrm{1}-{\displaystyle \frac{\overline{U}\left(r\right)}{U_{\mathrm{\infty }}}}\right)r\mathrm{d}r,\end{array}$$

where r₁ and r₂ are then found by solving the following algebraic system.

$$\begin{array}{}\text{(13a)}& {\displaystyle }\mathit{\delta }\left({\stackrel{\mathrm{‾}}{U}}_{\text{pw}}\right)=\mathit{\delta }\left({\stackrel{\mathrm{‾}}{U}}_{\text{RANS}}\right)\text{(13b)}& {\displaystyle }{\mathit{\delta }}_{\mathit{\theta }}\left({\stackrel{\mathrm{‾}}{U}}_{\text{pw}}\right)={\mathit{\delta }}_{\mathit{\theta }}\left({\stackrel{\mathrm{‾}}{U}}_{\text{RANS}}\right)\end{array}$$

Figure 3 Comparison of the step profile with the hyperbolic tangent profile.

An additional simplification is possible if we assume that the wake shear region remains small relative to the size of the rotor diameter. In this case, we can decompose the ${\stackrel{\mathrm{‾}}{U}}_{\text{pw}}$ profile into

$$\begin{array}{}\text{(14)}& {\stackrel{\mathrm{‾}}{U}}_{\text{pw}}\left(r\right)={\stackrel{\mathrm{‾}}{U}}^{\left(\mathrm{0}\right)}\left(r\right)+{\stackrel{\mathrm{‾}}{U}}^{\left(\mathrm{1}\right)}\left(r\right),\end{array}$$

where ${\stackrel{\mathrm{‾}}{U}}^{\left(\mathrm{0}\right)}$ is the Heaviside step function,

$$\begin{array}{}\text{(15)}& {\stackrel{\mathrm{‾}}{U}}^{\left(\mathrm{0}\right)}\left(r\right)=\left\{\begin{array}{ll}{U}_{\mathrm{0}},& r\le r_{\mathrm{e}}\\ U_{\mathrm{\infty }},& r>r_{\mathrm{e}},\end{array}\right.\end{array}$$

and ${\stackrel{\mathrm{‾}}{U}}^{\left(\mathrm{1}\right)}$ is a small perturbation to the single step profile:

$$\begin{array}{}\text{(16)}& {\stackrel{\mathrm{‾}}{U}}^{\left(\mathrm{1}\right)}\left(r\right)=\left\{\begin{array}{ll}\mathrm{0},& r<r_{\mathrm{1}}\\ +\mathrm{\Delta }U,& r_{\mathrm{1}}\le r\le r_{\mathrm{e}}\\ -\mathrm{\Delta }U,& r_{\mathrm{e}}\le r\le r_{\mathrm{2}}\\ \mathrm{0},& r>r_{\mathrm{2}}.\end{array}\right.\end{array}$$

This assumption allows the analytical results of Batchelor and Gill (1962) to be directly applied with some minor modifications, as discussed in the following section.

2.4.2 Spatial linear stability formulation

For the flow variables $\stackrel{\mathrm{̃}}{\mathit{\varphi }}=[\stackrel{\mathrm{̃}}{u},\stackrel{\mathrm{̃}}{v},\stackrel{\mathrm{̃}}{w},\stackrel{\mathrm{̃}}{p}]$, where $\stackrel{\mathrm{̃}}{u}$, $\stackrel{\mathrm{̃}}{v}$, and $\stackrel{\mathrm{̃}}{w}$ are the streamwise, radial, and azimuthal velocities, respectively, and $\stackrel{\mathrm{̃}}{p}$ is pressure, we assume that they can be expressed in terms of the radial eigenfunctions, ${\widehat{\mathit{\varphi }}}_n$, and the complex exponential basis functions,

$$\begin{array}{}\text{(17)}& \stackrel{\mathrm{̃}}{\mathit{\varphi }}(x,r,\mathit{\psi },t)={\widehat{\mathit{\varphi }}}_n\left(r\right)e^{i\mathit{\alpha }x+in\mathit{\psi }-i\mathit{\omega }t},\end{array}$$

where $\mathit{\alpha }={\mathit{\alpha }}_r+i{\mathit{\alpha }}_i$ is the complex streamwise wavenumber, n is the azimuthal mode number, and ω is the temporal frequency. Here the values of n and ω are taken to match the azimuthal mode number and Strouhal frequencies used in the desired blade pitching strategy from Eq. (1) and described in Sect. 2.5.1 and Table 3. In this study, we consider the impact of a single instability wave at a single Strouhal number and a specific n on the turbine wake development. For the more general case, multiple instability wave components can be included in the analysis, and a summation over all wave components is required in Eq. (17).

Table 3 AWM parameters.

n/a: not applicable

Assuming an inviscid, parallel flow with a piecewise constant velocity profile, the governing equations for mass and momentum conservation of the wave components, $\stackrel{\mathrm{̃}}{\mathit{\varphi }}$, are

$$\begin{array}{}\text{(18a)}& {\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{u}}{\partial x}}+{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{v}}{\partial r}}+{\displaystyle \frac{\stackrel{\mathrm{̃}}{v}}{r}}+{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{w}}{\partial \mathit{\psi }}}=\mathrm{0},\text{(18b)}& {\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{u}}{\partial t}}+\stackrel{\mathrm{‾}}{U}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{u}}{\partial x}}=-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{p}}{\partial x}},\text{(18c)}& {\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{v}}{\partial t}}+\stackrel{\mathrm{‾}}{U}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{v}}{\partial x}}=-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{p}}{\partial r}},\text{(18d)}& {\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{w}}{\partial t}}+\stackrel{\mathrm{‾}}{U}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{w}}{\partial x}}=-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }r}}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{p}}{\partial \mathit{\psi }}}.\end{array}$$

Here, $\stackrel{\mathrm{‾}}{U}\left(r\right)={\stackrel{\mathrm{‾}}{U}}_{\text{pw}}$, with the parameters r₁, r₂, U₀, and U_∞ chosen to match the characteristics of the RANS wake profiles. Inserting the representation (Eq. 17) into Eqs. (18) leads to the following spectral versions of the governing equations:

$$\begin{array}{}\text{(19a)}& {\displaystyle }i\mathit{\alpha }{\widehat{u}}_n+{\displaystyle \frac{{\stackrel{\mathrm{̃}}{v}}_n}{r}}+{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{v}}_n}{\partial r}}+{\displaystyle \frac{in}{r}}{\stackrel{\mathrm{̃}}{w}}_n=\mathrm{0},\text{(19b)}& {\displaystyle }\mathit{\xi }{\widehat{u}}_n+{\widehat{v}}_n{\displaystyle \frac{\mathrm{d}\stackrel{\mathrm{‾}}{U}}{\mathrm{d}r}}+{\widehat{f}}_n^{\left(\mathrm{1}\right)}=-i{\displaystyle \frac{\mathit{\alpha }}{\mathit{\rho }}}{\widehat{p}}_n,\text{(19c)}& {\displaystyle }\mathit{\xi }{\widehat{v}}_n=-{\displaystyle \frac{\mathrm{1}}{\mathit{\rho }}}{\displaystyle \frac{\mathrm{d}{\widehat{p}}_n}{\mathrm{d}r}},\text{(19d)}& {\displaystyle }\mathit{\xi }{\widehat{w}}_n=-i{\displaystyle \frac{n}{\mathit{\rho }r}}{\widehat{p}}_n,\end{array}$$

where $\mathit{\xi }\left(r\right)=\mathit{\alpha }{\stackrel{\mathrm{‾}}{U}}^{\left(\mathrm{0}\right)}\left(r\right)-\mathit{\omega }$. Note that the ${\stackrel{\mathrm{‾}}{U}}_{\text{pw}}$ profile has been decomposed according to Eq. (14) and the term ${\widehat{f}}_n^{\left(\mathrm{1}\right)}=i\mathit{\alpha }{\stackrel{\mathrm{‾}}{U}}^{\left(\mathrm{1}\right)}{\widehat{u}}_n$. The system of equations in Eq. (19) can be combined into the Rayleigh ordinary differential equation for the pressure, ${\widehat{p}}_n$,

$$\begin{array}{}\text{(20)}& {\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}r}}\left(r{\displaystyle \frac{\mathrm{d}{\widehat{p}}_n}{\mathrm{d}r}}\right)-\left[{\left({\displaystyle \frac{n}{r}}\right)}^{\mathrm{2}}+{\mathit{\alpha }}^{\mathrm{2}}\right]{\widehat{p}}_n=i\mathit{\rho }\mathit{\alpha }{\widehat{f}}_n^{\left(\mathrm{1}\right)}.\end{array}$$

Both the eigenfunctions, $\widehat{p}\left(r\right)$, and the eigenvalues, α, can be decomposed into zeroth-order and first-order components:

$$\begin{array}{}\text{(21a)}& {\displaystyle }{\widehat{p}}_n={\widehat{p}}_n^{\left(\mathrm{0}\right)}+{\widehat{p}}_n^{\left(\mathrm{1}\right)},\text{(21b)}& {\displaystyle }\mathit{\alpha }={\mathit{\alpha }}^{\left(\mathrm{0}\right)}+{\mathit{\alpha }}^{\left(\mathrm{1}\right)}.\end{array}$$

Here, both ${\widehat{p}}_n^{\left(\mathrm{1}\right)}$ and α⁽¹⁾ are assumed to be small relative to ${\widehat{p}}_n^{\left(\mathrm{0}\right)}$ and α⁽⁰⁾, respectively, and the solution can be found as part of an eigenvalue perturbation problem. Equation (20) can be similarly divided into the zeroth-order and first-order contributions, where only the appropriate order terms are included in the equations:

$$\begin{array}{}\text{(22a)}& {\displaystyle }\begin{array}{rl}\mathcal{L}\mathit{\left\{}{\widehat{p}}_n^{\left(\mathrm{0}\right)}\mathit{\right\}}& ={\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}r}}\left(r{\displaystyle \frac{\mathrm{d}{\widehat{p}}_n^{\left(\mathrm{0}\right)}}{\mathrm{d}r}}\right)-\left[{\left({\displaystyle \frac{n}{r}}\right)}^{\mathrm{2}}+({\mathit{\alpha }}^{\left(\mathrm{0}\right)}{)}^{\mathrm{2}}\right]{\widehat{p}}_n^{\left(\mathrm{0}\right)}\\ & =\mathrm{0},\end{array}\text{(22b)}& {\displaystyle }\mathcal{L}\mathit{\left\{}{\widehat{p}}_n^{\left(\mathrm{1}\right)}\mathit{\right\}}+\mathrm{2}{\mathit{\alpha }}^{\left(\mathrm{0}\right)}{\mathit{\alpha }}^{\left(\mathrm{1}\right)}{\widehat{p}}^{\left(\mathrm{0}\right)}=i\mathit{\rho }{\mathit{\alpha }}^{\left(\mathrm{0}\right)}{\widehat{f}}_n^{\left(\mathrm{1}\right)}.\end{array}$$

The zeroth-order solution to Eq. (22a) is given by the modified Bessel functions:

$$\begin{array}{}\text{(23)}& {\widehat{p}}_n^{\left(\mathrm{0}\right)}\left(r\right)=\left\{\begin{array}{ll}{C}_{\mathrm{1}}{I}_n\left({\mathit{\alpha }}^{\left(\mathrm{0}\right)}r\right),& r<r_{\mathrm{e}}\\ C_{\mathrm{2}}{K}_n\left({\mathit{\alpha }}^{\left(\mathrm{0}\right)}r\right),& r\ge r_{\mathrm{e}}.\end{array}\right.\end{array}$$

The constants, C₁ and C₂, are chosen so the pressure is continuous at r=r_e, and the kinematic condition for the displacement, η, of a material line at r=r_e is also satisfied:

$$\begin{array}{}\text{(24)}& {\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{\mathit{\eta }}}{\partial t}}+\stackrel{\mathrm{‾}}{U}{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{\mathit{\eta }}}{\partial x}}=\stackrel{\mathrm{̃}}{v}.\end{array}$$

Assuming the functional form $\stackrel{\mathrm{̃}}{\mathit{\eta }}(x,\mathit{\psi },t)=\widehat{\mathit{\eta }}{e}^{i\mathit{\alpha }x+in\mathit{\psi }-i\mathit{\omega }t}$, enforcing above conditions leads to the following nonlinear relation which can be used to solve for α⁽⁰⁾ at every frequency ω:

$$\begin{array}{}\text{(25)}& {\displaystyle \frac{\mathit{\xi }(r_{\mathrm{0}}{)}^{\mathrm{2}}}{\mathit{\xi }(r_{\mathrm{\infty }}{)}^{\mathrm{2}}}}={\displaystyle \frac{K_n^{\prime }\left({\mathit{\alpha }}^{\left(\mathrm{0}\right)}{r}_{\mathrm{e}}\right)I_n\left({\mathit{\alpha }}^{\left(\mathrm{0}\right)}{r}_{\mathrm{e}}\right)}{K_n\left({\mathit{\alpha }}^{\left(\mathrm{0}\right)}{r}_{\mathrm{e}}\right)I_n^{\prime }\left({\mathit{\alpha }}^{\left(\mathrm{0}\right)}{r}_{\mathrm{e}}\right)}}.\end{array}$$

Up to this point, the analysis follows that of Batchelor and Gill (1962) for piecewise constant velocities and is shown to be valid for infinitely sharp, top-hat velocity profiles. However, in the current work, capturing the effects of the wake spreading are important to the growth and evolution of large-scale structures. This can be accomplished by including a small perturbation to the ${\stackrel{\mathrm{‾}}{U}}^{\left(\mathrm{0}\right)}$ profile and calculating the corresponding perturbation to the growth rates. Once ${\widehat{p}}_n^{\left(\mathrm{0}\right)}$ and α⁽⁰⁾ are known, the perturbation, α⁽¹⁾, to the wavenumber can be found by applying the inner product,

$$\begin{array}{}\text{(26)}& \langle f,g\rangle =\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}f\left(r\right)g\left(r\right)r\mathrm{d}r,\end{array}$$

to Eq. (22b), leading to

$$\begin{array}{}\text{(27)}& \begin{array}{rl}\langle \mathcal{L}\mathit{\left\{}{\widehat{p}}_n^{\left(\mathrm{1}\right)}\mathit{\right\}},{\widehat{p}}_n^{\left(\mathrm{0}\right)}\rangle +\langle \mathrm{2}{\mathit{\alpha }}^{\left(\mathrm{0}\right)}{\mathit{\alpha }}^{\left(\mathrm{1}\right)}{\widehat{p}}^{\left(\mathrm{0}\right)},& {\widehat{p}}_n^{\left(\mathrm{0}\right)}\rangle =\\ & \langle i\mathit{\rho }{\mathit{\alpha }}^{\left(\mathrm{0}\right)}{\widehat{f}}_n^{\left(\mathrm{1}\right)},{\widehat{p}}_n^{\left(\mathrm{0}\right)}\rangle .\end{array}\end{array}$$

Because ${\widehat{p}}_n^{\left(\mathrm{0}\right)}$ is self-adjoint and satisfies Eq. (22a), the terms in Eq. (27) can be rearranged into the following expression for α⁽¹⁾:

$$\begin{array}{}\text{(28)}& {\mathit{\alpha }}^{\left(\mathrm{1}\right)}=-{\displaystyle \frac{({\mathit{\alpha }}^{\left(\mathrm{0}\right)}{)}^{\mathrm{2}}{C}_{\mathit{\alpha }}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\overline{U}}^{\left(\mathrm{1}\right)}{\displaystyle \frac{{\widehat{p}}_n^{\left(\mathrm{0}\right)}{\widehat{p}}_n^{\left(\mathrm{0}\right)}}{{\mathit{\alpha }}^{\left(\mathrm{0}\right)}{\stackrel{\mathrm{‾}}{U}}^{\left(\mathrm{0}\right)}-\mathit{\omega }}}r\mathrm{d}r}{\mathrm{2}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\widehat{p}}_n^{\left(\mathrm{0}\right)}{\widehat{p}}_n^{\left(\mathrm{0}\right)}r\mathrm{d}r}},\end{array}$$

where C_α is a complex calibration constant. The full eigenvalue, α, can then be reconstructed through Eq. (21b). The real part of the wavenumber, α_r, determines the streamwise wavelength of the large-scale coherent structures, while the imaginary component, α_i, dictates the spatial growth of the structures.

To examine the accuracy of this asymptotic, analytic approach with a piecewise constant velocity profile, a comparison of the linear stability solution using a continuous hyperbolic tangent profile (Eq. 16) and the ${\stackrel{\mathrm{‾}}{U}}_{\text{pw}}$ profile from Eq. (10) is shown in Fig. 4. As expected, the dispersion relation, α=α(ω), shows excellent agreement between the approaches when the profile width, Δ, is relatively small. For larger values of Δ, the wavenumber, α_r, calculated using the piecewise constant profile approximation U_pw, still agreed with the hyperbolic tangent profile, although there were some discrepancies visible for the growth rate, α_i. However, for the lower frequencies of interest, the agreement between the approaches is still reasonably accurate.

Figure 4 Comparison of linear stability theory between the analytic tanh profiles (Eq. 16) and the piecewise step profiles (Eq. 10). The first helical mode, n=1, is shown in all cases.

2.4.3 Coupled solution process

The spatial linear stability formulation described in Sect. 2.4.2 can be easily integrated into the RANS solution process discussed in Sect. 2.3. At every streamwise location, x, the RANS velocity profile is first computed assuming F_CS=0. This velocity profile is then used in the linear stability formulation to compute the velocity eigenfunctions, $\widehat{u}$ and $\widehat{v}$, and the corresponding wavenumbers, α. The evolution of the wave component of the flow variables can be calculated using the formula

$$\begin{array}{}\text{(29)}& \begin{array}{rl}\stackrel{\mathrm{̃}}{u}(x,r,& \mathit{\psi },t)=\\ & a_n{\widehat{u}}_n\left(r\right)\exp \left(i\int \mathit{\alpha }\left(x\right)\mathrm{d}x-in\mathit{\psi }-i\mathit{\omega }t\right),\end{array}\end{array}$$

where the integral over α in the exponential accounts for the slow changes in growth rate as the mean flow evolves. The initial amplitude of the wave component is given by a_n. Once $\stackrel{\mathrm{̃}}{u}$ and $\stackrel{\mathrm{̃}}{v}$ are known, the mean flow correction term, F_CS, can be calculated, and a new RANS velocity profile can be computed for the same x location. This process is repeated until the RANS velocity profiles meet specified convergence criteria (for which we use the Frobenius norm of the difference between two successive solutions being less than 10⁻⁷), after which the streamwise marching process proceeds to the next location at x+Δx. This process results in a two-way coupled model of the mean flow and the coherent flow structures, which differs from other dynamic approaches, such as the dynamic wake meandering model, that use a pre-determined mean flow to drive wake dynamics.

The initial formulation of both the RANS model and the linear stability model was implemented in Python and run on workstations with a single CPU for all cases. For typical cases, which used 200 grid points in the radial direction and 200 streamwise points, the baseline RANS calculation took 1–2 s to compute, and in cases with the RANS model coupled to the linear stability model, the total solve time was approximately 11–12 s.

2.5 AMR-Wind LES calculations

To generate the data necessary to calibrate the RANS model coefficients and evaluate the accuracy of the coupled RANS and linear stability approach, a series of LESs of turbine wakes was performed. These were done with the AMR-Wind code (Sharma et al., 2024; Sprague et al., 2020; Kuhn et al., 2025), a massively parallel, block-structured adaptive-mesh, incompressible flow solver for wind turbine and wind farm simulations. AMR-Wind solves the incompressible and low-Mach formulations of the Navier–Stokes equations, as well as temperature, subgrid-scale kinetic energy, and other scalar equations necessary for LESs of wind farms. AMR-Wind solves the discretized equations using a second-order finite method and second-order temporal integration. AMR-Wind includes all the necessary physics modules to simulate atmospheric boundary layers (ABLs). Included in this effort are ABL forcing, Boussinesq buoyancy, Coriolis forcing, body forcing to maintain the precursor-derived inflow condition in the presence of the turbine's blockage, and body forcing from coupling to OpenFAST (Jonkman et al., 2018; NREL, 2023) for turbine representation using actuator line models (these are the same forcing terms used in Brown et al., 2025, and Hsieh et al., 2025, for instance). AMR-Wind leverages AMReX for data structures, parallelism abstractions, and performance portability on heterogeneous architectures (Zhang et al., 2019). This framework has shown computational performance across a wide range of systems and applications (Fedeli et al., 2022; Henry de Frahan et al., 2022, 2024).

2.5.1 Turbine simulation parameters

Simulations using the IEA 15 MW reference turbine and the atmospheric conditions listed in Sect. 2.1 were performed in AMR-Wind using the one-equation k_sgs LES model (Moeng, 1984) and an actuator line model coupled to OpenFAST to represent the turbine blade forces. The simulation domains were either $\mathrm{4.5}\mathrm{km}\times \mathrm{2}\mathrm{km}\times \mathrm{1}\mathrm{km}$ (med-WS case) or $\mathrm{6.7}\mathrm{km}\times \mathrm{2}\mathrm{km}\times \mathrm{1}\mathrm{km}$ (low-WS and high-WS cases). In all cases, a background mesh resolution of 5 m was used, which was refined to a resolution of 2.5 m in a region 4.75D upstream to 12D downstream of the rotor, leading to mesh sizes of 179 M and 309 M, respectively. A time step of 0.02 s was used in the turbine simulations, and all simulations had a total runtime of at least 1000 s to allow the initial transients to dissipate and the wake structures to fully develop in the flow.

OpenFAST is a conglomeration of models that characterize the whole-system dynamics of a wind turbine including the aerodynamic and structural loading, drivetrain, and controls. Each of these individual models were tuned to match the specifics of the IEA 15 MW turbine. In addition, there are various submodels and parameters within the primary OpenFAST modules that can be adjusted to affect the fidelity, discretization, and properties of the simulations. For example, the aerodynamics module in OpenFAST consists of six submodels for rotor–wake induction, blade airfoil aerodynamics, tower influence on the blade, tower and nacelle drag, aeroacoustics, and buoyancy effects for floating substructure components. For OpenFAST, there were no wake, induction, hub-loss, or tip-loss models enabled as those phenomena were modeled within AMR-Wind. The blade airfoil structural dynamics were calculated through ElastoDyn and the unsteady aerodynamics were represented by the Beddoes–Leishman unsteady model with the Minnema–Pierce extension. Blade definition files are provided to the aerodynamics module to define the blade nodal discretization, geometry, chord, twist, airfoil identifiers, and buoyancy properties. A total of 50 and 20 nodes were used to represent the aerodynamics on the turbine blades and the tower, respectively. Similarly, the aerodynamic effects of the hub and nacelle were represented using actuator forcing with a drag coefficient of 0.5 and a representative area of 49.5 m². Precise subsampling by OpenFAST was performed with a subcycling ratio between OpenFAST and AMR-Wind of 4:1.

Additional simulation details, including the turbine model parameters and the AMR-Wind input files, can be found in the Supplement.

2.6 Blade pitch actuation for AWM

The earlier work of Cheung et al. (2024a) explored the connection between blade pitch actuation parameters used in AWM strategies and the instability modes considered in the current analysis. In that study, different blade pitch actuation strategies were applied to an OpenFAST turbine model simulation using different pitch amplitudes, azimuthal mode numbers, and the desired Strouhal frequency of St=0.30. An analysis of the resulting blade loads showed that there was a corresponding fluctuating streamwise blade force that appeared at the same azimuthal mode number and Strouhal frequency. Through a spectral POD analysis, we can see that these fluctuating streamwise blade forces then excite a similar response in the near wake (Yalla et al., 2025). Although there may be differences in the radial profiles between the streamwise forces induced by fluctuating blade motions and the eigenfunctions of the Rayleigh equation (Eq. 22), it is sufficient to pitch the blades at the specified azimuthal mode number n and Strouhal number to excite the desired instability mode.

In the baseline turbine simulations, no AWM strategy was employed and the wake was allowed to develop naturally without any blade pitch actuation. These baseline cases were compared to simulations in which the helix and pulse AWM strategies were used (Table 3). All AWM strategies used a single actuation frequency of St=0.30, which is consistent with the Strouhal forcing used in previous studies (Cheung et al., 2024a). The blade pitch amplitudes were set to either 2° or 4° in both the helix and the pulse AWM strategies to determine the relative effectiveness of each actuation strategy.

3 Results

To evaluate the accuracy of the RANS and linear stability wake model, we compare the modeled wake behavior with the corresponding wake behavior from the AMR-Wind simulations. Results are shown first for the baseline cases where no AWM strategy was used, which allows us to evaluate the underlying RANS model without any coupling to the linear stability model. This is followed by a discussion of the AWM cases with helix and pulse actuation and an evaluation of the full RANS plus linear stability model.

3.1 Baseline wake behavior

Comparisons of the hub-height velocity profiles between the RANS model and the AMR-Wind LES calculation for the med-WS and high-WS cases are shown in Figs. 5 and 6, respectively, for various downstream distances. The baseline wake behavior for the low-WS case was very similar to the med-WS case because the turbine was operating at the same thrust coefficient, so the low-WS comparisons are not shown below for the sake of brevity.

Figure 5 Comparisons of the normalized hub-height velocity profiles from the RANS model and AMR-Wind calculations for the baseline med-WS case.

Figure 6 Comparisons of the normalized hub-height velocity profiles from the RANS model and AMR-Wind calculations for the baseline high-WS case.

In the medium- to far-wake regions, for downstream distances of $x/D>\mathrm{3.0}$, good agreement is seen between the wake profiles from the RANS model and the AMR-Wind calculations. The general evolution of the wake deficit and the wake spreading behavior is well captured by the parabolized RANS model. The AMR-Wind wake profiles show evidence of veer effects, which causes asymmetry in the LES wake profiles. This effect is not captured by the RANS model due to the axisymmetric formulation, but the overall match between the methods remains high.

Very close to the turbine rotor some differences between the wake profiles are noticeable. For streamwise distances of $x/D<\mathrm{3.0}$, we see the influence of the hub and nacelle on LES wake profiles, which is not captured in the RANS model. The actuator line representation of the turbine in AMR-Wind more accurately models the aerodynamics near the hub and root sections of the blades, leading to a small recirculation zone immediately downstream of the nacelle. The simplified nature of the initial RANS profiles neglects these effects, as well as any asymmetries due to the interactions of swirl with shear and veer in addition to speedup of the ambient flow from wake blockage. However, despite these approximations, the RANS model still accurately captures the velocity shear near the wake edges.

A comparison of the centerline and rotor-averaged velocities, shown in Figs. 7 and 8, provides a similar picture of the RANS model's accuracy for the baseline wake cases. In the far-wake region, the RANS model accurately predicts the recovery of the centerline and rotor-averaged velocities. Very close to the rotor, the RANS model assumes the presence of a potential core region in the wake, which is not realistic, so it is unsurprising that the centerline velocities do not agree until $x/D\approx \mathrm{4.0}$.

Figure 7 Comparison of the normalized centerline velocities for the baseline turbine wakes.

Figure 8 Comparison of the normalized rotor-averaged velocities for the baseline turbine wakes.

Some differences are observed in the RANS behavior between the high-WS and the low/med-WS cases. One noticeable difference is that the potential core region is correctly modeled in the low- and medium-WS cases but overestimated in the high-WS case. This is reflected in the comparisons of Figs. 7 and 8, as well as the wake recovery in the velocity profiles in Fig. 6, which show lower centerline and rotor-averaged velocities for the RANS high-WS case in the far wake. We believe that these discrepancies can be reduced through improvements in the RANS model and additional calibration across a wider variety of wake cases in future work. It is acknowledged that a comparison to an axisymmetric LES of just the turbine rotor would have yielded a direct comparison with RANS. However, the objective of this work is to provide a usable, proof-of-concept framework that illustrates how a RANS model with a linear stability model can capture most of the phenomena of interest in the LES data. The long-term goal is to build it up from common principles towards being able to capture increasing physics complexity, such as shear effects, veer, and asymmetry. The differences between the RANS and LES discussed here and in other sections of the paper highlight that the current approach performs well in comparison with complex LES data while also pointing towards future improvements.

3.2 Comparisons for AWM cases

With the application of an AWM strategy, we expect the turbine wake to mix faster due to the presence of the large-scale coherent structures. As shown in Figs. 9 and 10, the hub-height velocity profiles for the LES calculations and the RANS with a linear stability model indicate a faster wake recovery and increased mixing in the downstream wake. For the med-WS case in Fig. 9, there was qualitative agreement between the LES calculations and RANS with a linear stability model in predicting the changes to the wake width and centerline velocity for both the helix and pulse AWM cases and for both 2° and 4° pitch actuation. For the helix AWM case at the high-WS condition (Fig. 10), the LES calculations show more impact to the centerline velocity recovery, although the RANS with a linear stability model still shows the changes to the wake width due to AWM.

Figure 9 Comparisons of the hub-height velocity profiles from AMR-Wind and RANS with a linear stability model for the med-WS case with helix AWM (b, d, f) and pulse AWM (a, c, e) actuated.

Figure 10 Comparisons of the hub-height velocity profiles from AMR-Wind and RANS with a linear stability model for the high-WS case with helix AWM.

The rotor-averaged velocity comparisons in Figs. 11 and 12 show a similar level of agreement between the AMR-Wind and RANS with a linear stability model in the far wake. For downstream distances of $x/D>\mathrm{5.0}$, the RANS with a linear stability model qualitatively captured the wake recovery benefits for both the helix and pulse approaches. It was also seen that the helix AWM was not as effective as in the med-WS case as it was in the high-WS condition. This is attributed to the fact that the turbine operates at a lower thrust coefficient at the higher wind speeds, resulting in less initial wake deficit and lower velocity shear near the wake edges. The lower shear in the turbine wake translates to slower growth of the large-scale coherent structures, meaning that there is less opportunity for them to mix the turbine wake and impact the flow.

Figure 11 Comparisons of the normalized rotor-averaged velocities from the low-WS, med-WS, and high-WS cases using the helix AWM strategy with 2° and 4° actuation amplitudes.

Figure 12 Comparisons of the normalized rotor-averaged velocities from the low-WS and med-WS cases using the pulse AWM strategy with 2° and 4° actuation amplitudes.

Some differences between the LES calculations and the RANS with a linear stability model are observed in the near-wake region. For the low-WS and the med-WS cases, the growth of the coherent structures in the LES calculations is faster than the RANS with a linear stability model, so the wake benefits to the rotor-averaged velocity also appear earlier in the flow. However, in both the LES and the RANS with a linear stability model, the growth of the large-scale structures saturates at similar levels downstream, so the final wake benefits in the far wake remain comparable.

Note that the observed AWM model behavior for the high-WS case is consistent with earlier observations regarding the RANS model predictions for that condition. In Fig. 10, the larger potential core region in the RANS profiles limits the modifications from the coherent structures to the wake shear regions until father downstream in the wake. This leads to relatively minor changes to the centerline velocity for the high-WS case compared to the med- or low-WS cases (Fig. 9) and suggests that accurately capturing the mean flow is critical to modeling the impact of large-scale structures on wake behavior.

A quantitative measure for the accuracy of the RANS and linear stability approach to modeling AWM effects is provided in Table 4. Comparisons of the minimum hub-height streamwise velocities between the AMR-Wind calculations and the modeled wake predictions are shown at $x/D=\mathrm{8}$ and $x/D=\mathrm{9}$. For the med- and low-WS cases with helix and pulse forcing, the majority of the velocity errors are below 5 %, and, as expected, the largest differences compared to AMR-Wind occurred for the high-WS cases.

Table 4 Comparisons of the minimum hub-height streamwise velocity, $\text{min}\left(U\right)/U_{\mathrm{\infty }}$, between the AMR-Wind calculations and the RANS with a linear stability model at $x/D=\mathrm{8}$ and $x/D=\mathrm{10}$. The error is defined as $\mathit{ϵ}=U_{\text{RANS}}/U_{\text{AMRWIND}}-\mathrm{1}$.

3.2.1 Phase-averaged velocity

Some insight into the behavior of the large-scale structures can be seen in Fig. 13, where contours of the mean and phase-averaged velocity fields are shown for the low-WS case. As expected, the near-wake region of the RANS and linear stability model fails to capture some flow features immediately downstream of the rotor disk. However, in the far-wake region, the differences in the coherent structures between the pulse AWM forcing and the helix AWM forcing become apparent. The wave components in the pulse case form axisymmetric structures with a wavelength of approximately 2D, while a spiral pattern appears in the helix AWM cases with similar wavelength. In the AMR-Wind simulations the structures emerge earlier in the turbine wake but eventually saturate and impact the wake in a qualitatively similar way compared to the RANS and linear stability model predictions.

Figure 13 Comparisons of the AMR-Wind and RANS with the linear stability mean velocity component (a, b) with the phase-averaged velocity for the HelixA4 (c, d) and PulseA4 (e, f) cases. In all cases the normalized streamwise velocity, $U/U_{\mathrm{\infty }}$, is plotted for the low-WS case.

4 Conclusions

In this work, a framework for modeling AWM was developed that accounts for the effects of both the large-scale coherent structures and the turbulence on the mean flow. Using a triple-decomposition approach, the turbine wake flow was separated into a time-averaged mean flow, fine-scale turbulent, and phase-averaged components, and a computationally efficient method for solving these components was formulated. An axisymmetric, parabolized k–ε RANS model was used to solve for the mean flow and fine-scale turbulence components. To model the wave components of the flow, a simplified, inviscid, parallel-flow, linear spatial stability analysis was used. The linear stability modes were coupled with the RANS model to capture the interactions between the coherent structures and the mean flow.

Comparisons between the RANS with a linear stability model and high-fidelity LES calculations of the turbine wakes showed that this framework was able to capture the wake modifications due to AWM actuation, particularly in the far-wake regions. Additional wake mixing and more rapid wake recovery were observed for both the pulse and helix AWM strategies. Some differences are also observed in the near-wake region of the flow. The high-fidelity LES calculations include non-axisymmetric features immediately downstream of the rotor, which the baseline RANS model fails to capture, and the impact of the large-scale coherent structures is also more evident in the near-wake region of the LES cases.

There are several limitations associated with the current formulation of the RANS model that could be improved in future studies. One significant constraint is the axisymmetric assumption, which restricts the model's applicability in realistic environments, such as ABLs characterized by large veer. To address this limitation, a parabolic three-dimensional version of the model could be developed using the same principles established here by marching two-dimensional wake profiles downstream, similar to the improvements suggested by Cheung et al. (2024b). Furthermore, incorporating direct interactions between the coherent flow structures and the turbulence may result in a more accurate representation of the flow dynamics, although the effect on the mean flow from these higher-order interactions remains unknown. Lastly, this work has demonstrated the importance of calibrating the RANS model constants; however, further work is needed to establish robust values for these constants, including the use of additional data in the calibration process.

Additional improvements to the linear stability model are also possible. The current model ignores the effects of swirl, shear, and veer, which impacts the growth of the coherent structures and can help improve the comparisons with the high-fidelity simulations. It may also be possible to use the full RANS velocity profile in the Rayleigh equation instead of a piecewise constant approximation, which would help improve near-wake predictions. Finally, future work may also investigate the impact of mode-to-mode interactions in a nonlinear stability framework. For example, the interactions between the axisymmetric modes and helical modes may be crucial in determining the optimal forcing strategy and could be worth exploring in later studies. Finally, additional insight may be gained by comparing the coherent structures' behavior from linear stability theory with modes computed by spectral POD analysis (Yalla et al., 2025). This may provide some indications of which effects are important for the stability model to capture.