Abstract

The present work aims to study the effect of veer, namely, the height-dependent lateral deflection of wind velocity due to Coriolis acceleration, on the coherent structures in the wake of the national renewable energy laboratory 5-MW reference wind turbine using the sparsity-promoting dynamic mode decomposition (SP-DMD) method for the detection of dynamically relevant flow structures. Large eddy simulation (LES) of the flow impacting the wind turbine is carried out by solving the filtered Navier–Stokes equations for incompressible flows. The effects of both tower and nacelle are included in the simulation using the immersed boundary method. Simulations are performed at Re108, generating the inlet velocity profile by a precursor simulation of the atmospheric boundary layer (ABL) subjected to Coriolis acceleration. The analysis of the dynamic mode decomposition (DMD) results allows to study the effects of the Coriolis acceleration on the most relevant dynamic modes in the turbine wake and to understand the basic mechanisms by which the wind veer significantly affects the wake recovery rate. Moreover, as a result of the SP-DMD methodology, the most relevant modes are extracted from the wake and a limited subset of relevant flow features that optimally approximates the original data sequence is identified. This small number of modes represent the kernel that can be employed for developing an accurate reduced order model of the wake.

1 Introduction

Wind energy is one of the pillars of the strategy for limiting global surface temperature increase to 1.5 °C above pre-industrial levels by the end of this century. According to a recent IRENA report, to achieve this task, we need to expand more than threefold the global installed renewable power generation capacity from 3382 GW in 2022 to 11,174 GW in 2030, the final contribution of wind energy being about 32% [1].

Improving the efficiency of new turbines and wind farms is mandatory to meet the above goal limiting the impact on landscape by reducing the occupancy of land and sea.

In the case of on-shore applications, single-unit power production of 5–6 MW is foreseen in the future, with respect to an average installed power of 2.6 MW in 2018, with diameters of 160–170 m. In the case of off-shore installations, a single wind turbine will achieve 15–20 MW in the next 10–20 years, with diameters of about 250 m and hub heights of 150 m [2].

In wind farms, the wind turbines influence each other in a complex way through the action of their wakes [3,4]. For instance, the presence of a turbine-placed upstream leads to an efficiency reduction of 40% for a single turbine in off-shore applications. This interaction becomes stronger when one attempts to increase the power density of the wind farm by reducing the distance between the turbines. Optimizing the design of the farm, for example, in terms of turbine layout and control, is therefore crucial, and their control is essential to ensure high efficiency, flexibility, and security of power generation.

Moreover, considering that the hub heights of innovative wind turbines are greater than 100 m, the rotors operate within the atmospheric boundary layer (ABL) and a streamwise velocity gradient along the vertical direction characterizes the incoming flow. Such a gradient influences the blade loading, increasing the load on the blades pointing upwards in the ABL [5]. Therefore, any attempt to optimize the performance of a wind farm, taking into account the interactions among the turbines, should include the effect of the ABL gradient and its modifications.

In fact, the properties of the ABL impacting the turbine significantly influence both the operational state of a turbine and the dynamics of its wake, including the characteristics of the turbulence generated by the turbine and the recovery of the wake [3]. Numerical [6] and experimental data [7,8] suggest that individual turbine wakes exhibit faster recovery when subjected to heightened levels of turbulence intensity within the incoming flows.

The standard nonuniform velocity profile of the ABL, which follows a logarithmic pattern under neutral conditions, along with its turbulence distribution, induces a higher level of shear production of turbulent kinetic energy in the upper part of the wake [9]. De-Cillis et al. [10] indicate that turbulence significantly alters the key flow structures that are most relevant in terms of dynamics within the wake of a wind turbine. For the laminar-inflow case, relevant dynamic structures in the wake have high wave-numbers and are spatially localized. Introducing turbulence at the inflow replaces the high-frequency modes within the wake with lower-frequency modes that align with the natural oscillations of the wake, spanning the entire domain. These lower-frequency modes are primarily those emphasized in the precursor simulation. This suggests that in the presence of inflow turbulence, coherent structures associated with internal mechanisms like tip and root vortices become less dynamically significant compared to structures induced by external turbulence. This underscores the importance of incorporating atmospheric turbulence into low-dimensional turbine wake models.

Moreover, it is widely recognized that the rate of wake recovery is influenced by the Coriolis acceleration resulting from Earth's rotation [11], which induces an Ekman spiral flow within the surface layer of the ABL. This induces a height-dependent spanwise velocity component of the incoming wind, called wind veer, which can significantly influence the wind farm power output, as shown by Gadde et al. [12] using a large eddy simulation (LES) of a wind farm with an actuator disk model. Earlier studies [13] have explored the impact of veer effects on wind turbines and led to the improvement of a Gaussian wake model by introducing a correction term to account for veer. This updated model effectively predicts the altered wake structure characterized by skewness and shear resulting from the spanwise wind veer.

Howland et al. [14] considered the impact of turbine yaw and veer on blade aerodynamics. They introduced a model that accounts for the power generation of a wind turbine when misaligned with the wind direction. This model utilizes aerodynamic blade elements and integrates changes in wind speed and direction across the rotor area. Narasimhan et al. [15] studied the evolution of turbine wakes in the presence of both yawing and veer, and included these effects in an analytical model of the turbine wake. They employed LES, considering an isolated turbine rotor modeled by the actuator disk method. However, the effect of veer on the frequencies and wavelengths of the most dynamically relevant coherent flow structures is still to be determined.

The present work aims to investigate the effect of the Coriolis acceleration in a turbulent neutral ABL on the wake of the national renewable energy laboratory (NREL)-5 MW reference wind turbine using a sparsity-promoting dynamic mode decomposition (SP-DMD) [16] for the detection of dynamically relevant flow structures. For the first time, unlike the previous numerical works on this topic, we consider the complete turbine geometry including the tower and the nacelle. Concerning the dynamic mode decomposition (DMD) method, the basic theory of DMD was first developed by Rowley et al. [17] for computing the Koopman expansion of a sequence of flow-field snapshots, and then improved by Schmid et al. [18]. The DMD method consists of finding the eigenvalues and the eigenvectors of a linear operator approximating the nonlinear dynamics embedded in the data sequence and it has been recently exploited for the formulation of reduced order models of wind-turbine flows [19]. Using DMD, time-resolved data, such as a series of snapshots collected from either numerical simulations or experimental investigations, are analyzed to identify the principal coherent structures governing the flow dynamics, known as the DMD modes. These modes are not orthogonal and are characterized by a singular time-frequency, allowing them to be categorized based on their spectral contribution and spatial growth rate. This can be particularly interesting for designing a wind farm layout since one can analyze which modes conserve a relevant energy contribution downstream of a wind turbine and which frequencies are involved in the dynamic load acting over downstream turbines.

Recently, Kleine et al. [20] have provided a numerical study of the tip vortex instability in two in-line turbines using the actuator-line method and the DMD approach. The unsteady flow over a two-blade wind turbine model was computed by Sun et al. [21] through LES. They investigated the wake flow dynamics using the DMD method, analyzing it in both absolute and relative frames of reference with respect to the rotor. Typically, only a restricted set of DMD modes holds significance for turbine-turbine interactions, enabling the simplification of interaction models for optimization purposes [22]. However, using nonalgorithmic criteria to sort the modes and choosing a subset of the most significant ones may overlook important dynamics of the phenomena being studied and result in low accuracy and predictive capability of the generated low-dimensional models.

Thus, in the present work, we use the sparsity promoting [16] algorithm for ranking the most relevant DMD modes. The data sequence is composed of snapshots obtained by the LES, employing the actuator line method to simulate the rotor and the immersed boundary method to model the tower and the nacelle [10,23]. The wind-turbine flow is simulated imposing a realistic incoming ABL in neutral conditions, generated by a precursor LES including the effect of the Coriolis acceleration.

2 Methodologies

2.1 Large Eddy Simulations.

The wake behind the NREL-5 MW wind turbine is computed by LES, solving the following set of Navier–Stokes equations:
(1)
(2)

where u is the filtered nondimensional velocity; Re=UrefDν is the Reynolds number; p is the pressure; τij is the modeled subgrid-scale stress tensor; i and j1,2,3 indicate the streamwise (x), wall-normal (y) and spanwise (z) nondimensional directions, respectively; fc is the Coriolis coefficient, chosen equal to 1.25×104Hz, typical of North Sea latitude; ugj is the imposed geostrophic velocity component; FPI,i is the correction set by the proportional-integral (PI) wind controller similar to that used in Refs. [12,13,15,24], which is used only for the precursor simulation (see  Appendix); Fi a forcing term representing aerodynamic forces per unit volume exerted by the turbine blades on the fluid. Re is defined using the main characteristics of the wind turbine: Uref is the inlet velocity at hub height, D the rotor diameter, and ν the kinematic viscosity of the fluid. To capture the influence of small-scale structures on the filtered flow, the Smagorinsky model [25] with a constant Cs = 0.09 is employed. The simulations discussed in the present paper were conducted using a code developed at UT Dallas, called UTD-WF [26]. The code has been validated [27] using the numerical results provided in Ref. [28] concerning the NREL 5 MW reference turbine and also with respect to experimental wind tunnel data [27]. The governing equations are discretized in space through second-order centered finite differences on a staggered Cartesian grid. Concerning time integration, a third-order-accurate Runge–Kutta scheme is employed for the nonlinear terms and the Crank–Nicolson scheme for the linear terms. Therefore, the time-step limitation comes from the maximum Courant-Friedrichs-Lewy (CFL) condition for the explicit scheme, which, for linear problems, is CFL3. Moreover, another constraint to the CFL derives from the blade tip motion, which requires that CFL1/TSR.

The actuator line method [29] is employed to simulate the aerodynamic forces, whose values are calculated using the lift, CL, and drag, CD, coefficients of the rotor blades. The rotor blades are approximated as rotating rigid lines divided into discrete segments. For each element, the lift and drag forces are estimated as
(3)
(4)

where ρ is the fluid density, urel is the relative inflow velocity, α is the angle of attack, and F is the modified Prandtl correction factor, which considers the aerodynamic effect of the tip and root vortices. A Gaussian distribution kernel is then used to distribute these forces across areas perpendicular to each actuator line.

The tower and the nacelle are described using the immersed boundary method, avoiding the use of a body-fitted grid and reducing the computational cost of the simulations [10,30,31].

2.2 Dynamic Mode Decomposition.

Dynamic mode decomposition [17,18] is a mathematical technique used in the analysis of dynamic systems, particularly in the fields of fluid dynamics. It assumes that the evolution of the system can be approximated by a linear transformation between successive snapshots, allowing one to identify dynamically relevant flow structures, known as dynamic modes, by analyzing the inherent patterns in the simulated data. Regarding the data collection, this study acquires snapshots of the flow field from numerical simulations. Each snapshot, denoted as qiRN, is gathered at a consistent sampling frequency. The dimension of these snapshots (N) is determined by the product of the number of observables (O), which in this case are the three dimensionless components of the velocity, and the number of measurement points (S), which is the number of the grid points selected for the DMD analysis (see Sec. 4.2). Assuming a linear time-invariant mapping A and M snapshots, it is possible to express the relationship between two successive snapshots as follows:
(5)
Thus, defining the snapshot matrix
(6)
and combining Eqs. (5) and (6), one has
(7)
In order to render the computational procedure robust and affordable, a low-dimensional representation of A is needed. Toward this aim, we project the snapshot matrix Q0 on a low-rank orthogonal basis found by a singular value decomposition Q0=UΣVT, known as proper orthogonal decomposition (POD). For a given POD subspace dimension r, the obtained optimal rank-r approximation of Q0 is defined as
(8)
where U˜ and V˜ denote the first r leading columns of U and V and Σ˜ contains the leading r × r sub-block of Σ. This corresponds to the solution of the following optimization problem:
(9)
where S˜ is defined by the similarity transformation
(10)
(11)
where S˜ is defined by the similarity transformation
(12)
and is finally obtained as
(13)
In order to extract the dynamic modes, the eigen decomposition of S˜ is first computed
(14)
and is finally obtained as
(15)
where yk and zk* are the left and right eigenvectors of S˜, respectively, normalized as yk*yk=1 and zk*yj=δij. The flow field in the POD subspace for the ith snapshot can be summed up as follows:
(16)
The DMD modes are computed as
(17)
so that Eq. (16) can be rewritten as
(18)
where αk is the amplitude of the kth DMD mode. Therefore, Q0 can be written in matrix form
(19)
where Vand is the Vandermonde matrix, which contains the r eigenvalues of S˜ for each snapshot. Using Eq. (8), the amplitudes αk are determined by minimizing the error of the approximation Eq. (19), defined by the following Frobenius norm:
(20)

2.3 Sparsity Promoting Dynamic Mode Decomposition.

This study employs a version of DMD known as SP-DMD, as proposed by Jovanović et al. [16]. Unlike conventional DMD, SP-DMD has a distinct emphasis on achieving sparsity in the extracted dynamic modes, namely, it seeks a concise representation of the system's behavior. In fact, by constrained error minimization, only a limited number of DMD modes are selected. The approach used in this paper consists of enforcing a balance between the number of modes extracted and the approximation error relative to the entire data sequence, using the following modified objective function:
(21)

where γ is a positive regularization parameter given by the user. Larger values of γ correspond to sparser solutions. The alternating direction method of multipliers is applied to solve problem Eq. (21) because it is a convex optimization problem.

3 Simulation Setup

3.1 Precursor Simulation.

A precursor simulation is carried out to mimic a turbulent ABL, including the Coriolis effect. Considering the characteristics of the turbine, Re108 is selected. The computational domain has dimensions 6×11×3 diameter units in the streamwise, vertical, and spanwise directions, respectively, discretized using a 512×512×512 grid. No-slip boundary conditions are applied at the bottom wall, and free-slip conditions are applied at the top wall. Periodicity is applied in both the streamwise and spanwise directions. Given the boundary condition in the spanwise direction, a box width of three diameters is to be understood as equivalent to a row of infinite side-by-side turbines, spaced by the same length. Bernardoni et al. [32] investigated wind turbines arranged in clusters. Their findings indicate that side-by-side turbines spaced by three diameters exhibit no impact on each other's performance.

Figure 1 shows the time-averaged wind velocity profiles obtained by the precursor simulation. The black solid line is the ABL velocity profile; the blue dash-dotted line and the orange dashed line indicate the streamwise and spanwise velocity components, respectively; finally, the dotted line represents the Ekman spiral. A transverse wind shear is obtained due to the presence of the Coriolis effect.

Fig. 1
Time- and space-averaged dimensionless ABL velocity obtained from the precursor simulation (black solid line). The blue (dark) dash-dotted line is the streamwise component, the orange (light) dashed line is the spanwise component and the black dotted line is the Ekman spiral. (Color version online.)
Fig. 1
Time- and space-averaged dimensionless ABL velocity obtained from the precursor simulation (black solid line). The blue (dark) dash-dotted line is the streamwise component, the orange (light) dashed line is the spanwise component and the black dotted line is the Ekman spiral. (Color version online.)
Close modal
A wind direction PI controller (see  Appendix) is used to avoid a yaw misalignment between the wind turbine rotor and the mean velocity profile at hub height. In this way, the geostrophic wind is gradually adjusted by the addition of a source term to Eq. (1) [24]. The time-varying turbulent velocity obtained by the precursor simulation is imposed as an inlet condition for the turbine flow simulation. A turbulence intensity of 12% is obtained at the hub height, whereas the turbulence integral length scale has been evaluated by the following procedure. In light of the primary objective of this study, namely, to examine the impact of Coriolis acceleration on the wake structure, the integral length scale Luz of the incoming flow has been calculated considering the correlation function in the spanwise z-direction of the streamwise wind speed component. The time interval over which the analysis was conducted is approximately 1100 s. In order to evaluate the integral length scale at a specific height y from the wall, first the correlation function of a time series is computed at abscissa x as follows:
(22)
where u=uu¯ is the streamwise velocity fluctuation, the overline indicates time average, zk are the spanwise positions, Nk is the number of grid points in the spanwise direction, and δ and Δz are the separation distance and the grid spacing along the z-direction, respectively. Then, this procedure is repeated for every streamwise position and ρuuz(y,δ) is obtained through a spatial averaging in the x-direction. Given that periodic boundary conditions were applied to the side walls, the calculation of the lateral distance has been restricted to values of δ between 0 and Lz/2. The integral length scale is obtained as follows:
(23)

where δ̂ is the value corresponding to a ρuuz=0.05 as suggested by Flay et al. [33]. In Fig. 2(a), the correlation function relative to the hub height is shown. The value of the turbulence integral length scale at hub height Lu,hubz is equal to 27.4 m (0.22 diameter units), which is in good agreement with the value obtained by Stanislawski et al. [34] for the neutral condition case. The vertical profile of Luz(y) is shown in Fig. 2(b), focusing on the wind turbine region. Notably, the turbulence integral length scale exhibits an upward trend with increasing height, indicating a progressive enlargement of turbulent structures as altitude rises. Finally, for completeness, the power spectral density of the turbulence obtained by the precursor simulation is provided in Fig. 2(c). The power spectral density has been evaluated at a point at hub height and is independent of the x and z positions indicating that the turbulent broadband spectral energy distribution is fully developed.

Fig. 2
(a) Correlation curve at hub height of u′ versus the separation distance δ (the dashed line is the threshold value ρuu=0.05; the solid line indicates ρuu=0), (b) profile of the integral length scale versus the height, and (c) power spectral density of the turbulence obtained by the precursor simulation
Fig. 2
(a) Correlation curve at hub height of u′ versus the separation distance δ (the dashed line is the threshold value ρuu=0.05; the solid line indicates ρuu=0), (b) profile of the integral length scale versus the height, and (c) power spectral density of the turbulence obtained by the precursor simulation
Close modal

3.2 Turbine Flow Simulation.

Using the turbulent velocity obtained by the precursor simulation as inflow, we simulate the flow impacting the NREL-5 MW wind turbine, which has a rotor diameter of D=126m and a hub height of h=90m [35]. The wind turbine operates at its rated conditions, namely, with a tip-speed ratio TSR=7, resulting in a constant angular frequency of the rotor Ω=12.1 RPM and a reference velocity Uref=11.4m/s at the hub height. The computational domain has dimensions 12×11×3 diameter units in the x, y, and z directions, respectively. The computational grid is composed of 2048×512×512 grid points in x, y, and z, respectively. The grid is uniform along the streamwise and spanwise directions, whereas the grid points are stretched in the vertical direction. No-slip conditions at the bottom wall and free-slip conditions at the top wall are applied. Periodicity is applied in the spanwise direction. A radiative condition, using a uniform convection velocity c=0.9U is used as outlet condition [36]. At inlet points, the previously computed turbulent inflow, mimicking an ABL with the Coriolis effect, is imposed. Tower and nacelle are represented by the immersed boundary method [37], whereas the actuator line model is implemented to simulate turbine blades [31]. As shown in Fig. 3, the rotor is located at 4D from the inflow boundary and is centered in the spanwise direction. Finally, eight diameters downstream of the wind turbine are considered sufficient for allowing the wake to develop.

Fig. 3
Layout of the computational domain, in which the DMD subdomain is highlighted
Fig. 3
Layout of the computational domain, in which the DMD subdomain is highlighted
Close modal

4 Analysis of the Numerical Results

4.1 Mean Flow.

Through the analysis of the time-averaged velocity field, it is possible to evaluate the spatial development of the wind turbine wake structure. When the incoming wind is aligned in the streamwise direction, the wake is nearly symmetrical with respect to the vertical axis. Under geostrophic forcing, the rotor wake is highly skewed as it is advected downstream, due to the effect of a significant cross-flow component [38], as shown in Fig. 4. Moreover, the wake of the tower undergoes a deflection in the positive spanwise direction as one proceeds toward the outlet. This results in a further asymmetrical deformation of the rotor wake.

Fig. 4
Time-averaged streamwise velocity field obtained by LES, in the y–z planes at x = 5, 7, 9, and 11 diameter units downstream of the wind turbine
Fig. 4
Time-averaged streamwise velocity field obtained by LES, in the y–z planes at x = 5, 7, 9, and 11 diameter units downstream of the wind turbine
Close modal

Another approach to evaluate the development of the wake is to compute the velocity loss, defined as Δu=uinu, where uin is the time-averaged streamwise inflow velocity. This variable is provided in Fig. 5 showing that, as the wake proceeds along the streamwise direction, the velocity loss acquires a quasi-elliptical shape.

Fig. 5
Time-averaged velocity loss field obtained by LES, in the y–z planes at x = 5, 7, 9, and 11 diameter units downstream of the wind turbine
Fig. 5
Time-averaged velocity loss field obtained by LES, in the y–z planes at x = 5, 7, 9, and 11 diameter units downstream of the wind turbine
Close modal

Compared to the results of Abkar et al. [13], a stronger deformation is observed as the wake progresses along the domain, probably due to the presence of the tower.

4.2 Sparsity-Promoting Dynamic Mode Decomposition Results.

The SP-DMD analysis uses a data sequence consisting of 2370 three-dimensional snapshots of the velocity field, 1400 × 283 × 239 grid points in the DMD subdomain shown in Fig. 3, and a POD subspace of dimension r =400 for the projection step. The snapshots are sampled with a time interval Δt corresponding to 10deg rotation of the rotor. A three-dimensional subdomain is used to reduce the computational cost of this postprocessing. Sampling points are contained in a box with the following extent: [3.811.9]×[01.55]×[0.82.2] in the x, y, z directions. From now on, the origin of the reference system in the x direction is equal to the x coordinate of the turbine hub. The time-averaged velocity field is subtracted from each snapshot.

The eigenvalues μ of the linear operator S˜, computed as in Eq. (14), are shown in the left panel of Fig. 6. Results of the standard DMD algorithm are in orange, while those obtained by the sparsity-promoting algorithm are given by light blue circles for γ = 13,500, black circles for γ = 33,500, and red triangles for γ = 45,300. In all cases, as expected, the eigenvalues are distributed along the unitary circle. Moreover, the lower the γ, the higher the number of selected DMD modes and the consequent accuracy of the reconstruction of the flow field.

Fig. 6
Left panel: eigenvalues, μ, computed by the DMD (orange crosses) and by the SP-DMD for different sparsity levels, γ = 13,500 (light blue circles), γ = 33,500 (black circles) and γ = 45,300 (red triangles). Right panel: logarithmic mapping of the eigenvalues, ω, as defined in Eq. (24), for the same values of γ. The real part of ω, namely, R(ω), is the angular frequency, and the imaginary part, I(ω), is the growth rate of the DMD modes.
Fig. 6
Left panel: eigenvalues, μ, computed by the DMD (orange crosses) and by the SP-DMD for different sparsity levels, γ = 13,500 (light blue circles), γ = 33,500 (black circles) and γ = 45,300 (red triangles). Right panel: logarithmic mapping of the eigenvalues, ω, as defined in Eq. (24), for the same values of γ. The real part of ω, namely, R(ω), is the angular frequency, and the imaginary part, I(ω), is the growth rate of the DMD modes.
Close modal
In the right panel of Fig. 6, the logarithmic mapping of the eigenvalues μ is shown, in which ω is defined as follows:
(24)

where i is the imaginary unit and Δt is the time interval between two succeeding snapshots. The real part of ω is the angular frequency (temporal wavenumber) of the DMD mode, whereas the imaginary part is its growth rate. Every mode is also associated with a spatial structure, that will be shown later on, which oscillates in time at the given frequency ω. As the sparsity coefficient increases, the subset focuses more on low-frequency modes, selecting the coherent structures with higher amplitude α, as shown in Fig. 7. This suggests that high-frequency modes, such as tip vortices, are less significant for the flow dynamics in the presence of atmospheric turbulence.

Fig. 7
Absolute value of DMD amplitudes αi versus the angular frequency R(ω). The results are obtained by DMD (orange crosses) and by SP-DMD for different sparsity levels: γ = 13,500 (light blue circles), γ = 33,500 (black circles), and γ = 45,300 (red triangles).
Fig. 7
Absolute value of DMD amplitudes αi versus the angular frequency R(ω). The results are obtained by DMD (orange crosses) and by SP-DMD for different sparsity levels: γ = 13,500 (light blue circles), γ = 33,500 (black circles), and γ = 45,300 (red triangles).
Close modal

For γ = 33,500, the sparsity-promoting algorithm picks seven pairs of DMD modes. The presence of the Coriolis effect results in a significant wind veer. Therefore, it is interesting to analyze both streamwise and spanwise velocity isosurfaces of the seven selected modes. Some of the extracted modes are shown in Figs. 8 and 9, respectively. Table 1 provides for each DMD mode pair the amplitude, the frequency, and the main streamwise wavelength for the x and z component of the velocity, namely, λxu and λxw, respectively. The latter have been recovered by taking the larger-amplitude spatial wavenumber of the Fourier spectra in the streamwise direction, for the streamwise, u, and spanwise, w, components of the velocity, respectively. Such spectra are obtained by averaging the corresponding Fourier spectra computed at several wall-normal and spanwise planes.

Fig. 8
Streamwise velocity fluctuation isosurfaces for four modes selected by SP-DMD with γ = 45,300. The isosufaces colored in red (dark) indicate a positive fluctuation u = 0.0008, whereas, the gray (light) ones represent a negative fluctuation u=−0.0008: (a) pair 1, (b) pair 3, (c) pair 6, and (d) pair 7. (Color version online.)
Fig. 8
Streamwise velocity fluctuation isosurfaces for four modes selected by SP-DMD with γ = 45,300. The isosufaces colored in red (dark) indicate a positive fluctuation u = 0.0008, whereas, the gray (light) ones represent a negative fluctuation u=−0.0008: (a) pair 1, (b) pair 3, (c) pair 6, and (d) pair 7. (Color version online.)
Close modal
Fig. 9
Spanwise velocity fluctuation isosurfaces for four modes selected by SP-DMD with γ = 45,300. The isosufaces colored in red (dark) indicate a positive fluctuation w = 0.0008, whereas, the gray (light) ones represent a negative fluctuation w=−0.0008: (a) pair 1, (b) pair 3, (c) pair 6, and (d) pair 7. (Color version online.)
Fig. 9
Spanwise velocity fluctuation isosurfaces for four modes selected by SP-DMD with γ = 45,300. The isosufaces colored in red (dark) indicate a positive fluctuation w = 0.0008, whereas, the gray (light) ones represent a negative fluctuation w=−0.0008: (a) pair 1, (b) pair 3, (c) pair 6, and (d) pair 7. (Color version online.)
Close modal
Table 1

Dimensionless amplitudes, temporal frequencies, and main wavelengths in the streamwise and spanwise directions, respectively, of the selected complex conjugate dynamic mode pairs, ordered according to their amplitude |α|

|α|ωλxuλxw
Pair 171.760.857.457.45
Pair 252.560.997.456.21
Pair 349.221.843.393.39
Pair 426.193.721.551.50
Pair 525.222.872.072.07
Pair 617.953.351.621.60
Pair 78.8942.0000
|α|ωλxuλxw
Pair 171.760.857.457.45
Pair 252.560.997.456.21
Pair 349.221.843.393.39
Pair 426.193.721.551.50
Pair 525.222.872.072.07
Pair 617.953.351.621.60
Pair 78.8942.0000

Regarding the isosurfaces of the streamwise velocity component in Fig. 8, white structures are referred to as u=0.0008. In contrast, red ones are relative to u =0.0008 (notice that the amplitude of each DMD mode's velocities has no physical meaning per se since they should be multiplied by the amplitudes αi for reconstructing the data sequence). Figure 8(d) shows the SP-DMD mode with R(ω)=42.00, corresponding to an angular frequency equal to three times the turbine rated one.

This mode, characterized by a high-frequency content, is associated with the tip vortices, as also observed by De-Cillis et al. [10] in the absence of the Coriolis effect. For lower frequency modes in the range R(ω)[0.851.84], it can be seen in Figs. 8(a) and 8(b) that the coherent structures cover the whole domain since they are associated with spatial structures with large length scale. Whereas, structures of moderate length scale are seen in the range R(ω)[1.843.73] (see Fig. 8(c)), according to the value of λxu in Table 1. A smaller γ results in a greater number of low-frequency SP-DMD modes, as discussed before. Such modes are mostly characterized by intermediate-size coherent structures, although some high-frequency modes are selected as well (see, for instance, the blue circle on the right of Fig. 6). In general, a distinct geometric pattern due to the presence of the turbine cannot be seen in these modes. An exception is the high-frequency modes linked to tip-vortices, such as that depicted in Fig. 8(d), where the coherent structures originate at the tip and root of the blades. Concerning the Strouhal numbers St of the low-frequency modes, one can wonder whether they match the wake meandering values obtained in recent studies [39,40]. For modes with R(ω)1, as also found by De-Cillis et al. [10] in the absence of the Coriolis effect, a Strouhal number St0.2 is found, which is consistent with the typical frequency ranges of wake meandering defined by the above-mentioned articles. This statement suggests that wake meandering is more likely caused by an external mechanism, linked to the turbulent nature of the ABL, rather than an internal one linked to the turbine itself. Concerning the spanwise velocity component, the SP-DMD mode's isosurfaces are shown in Fig. 9. In this case, an absolute value of w =0.0008 was chosen, using red for the positive value and white for the negative value. Like the streamwise component, the mode with R(ω)=42.00 in Fig. 9(d) shows the typical coherent tip-vortex structures. For modes with frequencies R(ω)1.84 (see Figs. 9(a) and 9(b)), the coherent structure spans the entire domain. This is because these structures are linked to spatial patterns with large length scales. In the range R(ω)[1.843.73], structures of moderate length scale can be observed in Fig. 9(c), according to the value of λxw in Table 1. Figure 9(b) shows the spanwise velocity component of the second SP-DMD mode. It is peculiar how moving away from the rotor region, the structures extend toward the left region of the computational domain. Moreover, the same structures in the far-wake region show that the wake of the tower is not centered in the z-direction, as one would expect when looking at Fig. 4. One can see that the structures have an oblique shape, reminiscent of the deformation dictated by the Coriolis effect. An interesting aspect is that unlike what is observed in the streamwise velocity component (see Figs. 8(a) and 8(b)), the low-frequency spanwise structures do not show the influence of the tower, nacelle, and of the rotor, indicating that the veer does not act strongly on the dynamics of the turbine but on the wake itself. For frequencies higher than R(ω)=3.35, there is a non-negligible influence of the presence of the tower and nacelle, as shown by the isocontours of the spanwise velocity component (see Figs. 9(c) and 9(d)).

In general, one can say that, due to the effect of the Coriolis acceleration, the spanwise component of the velocity strongly contributes to the formation of coherent structures in the wake of the wind turbine. The DMD mode analysis shows that such structures are comparable in amplitude and extension to those generated by the streamwise component of the velocity. Therefore, it appears that for a velocity profile with a significant cross-flow velocity component, the momentum exchange in the spanwise direction cannot be neglected. Concerning the wavelengths of the DMD modes, it is observed that the two components have values that are very close to each other (see λxu and λxw in Table 1). This is again an indication of how coherent structures in both directions have an important effect on the evolution of the wake within the considered subdomain.

5 Conclusions

This work provides an analysis of the effect of wind veer due to Coriolis acceleration in the atmospheric boundary layer on the coherent structures in the wake of the NREL 5-MW reference wind turbine using the SP-DMD method. Large eddy simulation of the flow through the wind turbines is carried on by solving the filtered Navier–Stokes equations for incompressible flows. The most relevant modes are extracted from the wake by the SP-DMD method and a limited subset of relevant flow features that optimally approximate the original data sequence is identified. The most dynamically relevant DMD modes show that the spanwise component of the velocity contributes to the formation of coherent structures that are comparable in amplitude and extension to those generated by the streamwise component of the velocity. In the far wake, such structures have an oblique shape, reminiscent of the deformation dictated by the Coriolis effect. An interesting aspect is that the far wake is dominated by low-frequency structures mainly influenced by the characteristics of the inflow turbulent boundary layer, whereas the coherent structures generated by the turbine (for instance, tip and hub vortices) are dynamically less relevant. The DMD modes obtained by the present analysis can be employed for developing accurate reduced models of the wake for wind turbines taking into account the characteristics of the incoming turbulence in the boundary layer (intensity and length scale) and the effect of Coriolis acceleration.

Acknowledgment

Diffuser augmented Wind Turbines for URBan environments (DWTURB) and NEST—Network 4 Energy Sustainable Transition pertain to the National Recovery and Resilience Plan (PNRR), Mission 4 Component 2 Investment 1.3, funded from the European Union—NextGenerationEU. This paper reflects only the authors' views and opinions, neither the European Union nor the European Commission can be considered responsible for them.

Funding Data

  • PRIN 20229YJP33_01, “Diffuser augmented Wind Turbines for URBban environments (DWTURB)”.

  • Ministero dell'Istruzione, dell'Università e della Ricerca (award ID): NEST - Network 4 Energy Sustainable Transition (D.D. 1243 02/08/2022, PE00000021; Funder ID: 10.13039/501100003407).

Nomenclature

CS =

constant of Smagorinsky model

DMD =

dynamic mode decomposition

fc =

Coriolis coefficient

Fi =

forcing term of aerodynamic forces of wind turbine blades

FPI =

correction of PI wind controller

LES =

large eddy simulations

p =

pressure

POD =

proper orthogonal decomposition

r =

rank of POD subset

SP – DMD =

sparsity promoting dynamic mode decomposition

TSR =

tip speed ratio

u =

streamwise velocity component

ug =

imposed geostrophic velocity

v =

vertical velocity component

w =

spanwise velocity component

x =

streamwise direction

y =

vertical direction

z =

spanwise direction

Greek Symbols
α =

amplitude of DMD

γ =

sparsity promoting regularization parameter

λ =

wavelength

τ =

modeled subgrid-scale stress tensor

ω =

temporal wavenumber

Ω =

rotor angular speed

Dimensionless Groups
Re =

Reynolds number, UD/ν

Superscripts and Subscripts
in =

inflow

ref =

reference value

=

free stream value

Appendix: Wind Direction Proportional-Integral Controller

In ABL dynamics, the wind direction has a vertical dependence in accordance with the Ekman spiral, for this reason, achieving the desired flow direction at the hub height becomes challenging. A possible solution is to correct the wind direction by acting directly on the geostrophic wind angle. In order to apply this action on the ABL, a correction source term FPI is added in the momentum equation (see Eq. (1))
(A1)
where ωPI is the control parameter. This term causes a rotation of the flow field around the vertical of
(A2)
where u3 and u1 are the horizontally averaged velocity components along z- and x-directions, respectively. Since the objective is to have a yaw angle ϕ=0,u3 and u1 are calculated at hub height. The overall control function is
(A3)

where KP and KI are the proportional and the integral terms, respectively; e(t) is the error, namely, the yaw angle ϕ. After tuning the PI controller, it is decided to use KP = 0.1 and KI = 1. The controller updates the streamwise and spanwise components of the geostrophic velocity, ug1 and ug3 while keeping constant the geostrophic velocity magnitude, G=(ug12+ug32)1/2. The controller remains active until the flow reaches a steady-state condition when the error is smaller than a given threshold (ϕthreshold=0.1 deg). Such a controller was previously devised by Sescu et al. [24] and then widely used in this research field [12,13,15].

References

1.
IRENA
,
2023
, “
COP28, IRENA and GRA (2023), Tripling Renewable Power and Doubling Energy Efficiency by 2030: Crucial Steps Towards 1.5 °C
,” International Renewable Energy Agency, Abu Dhabi, UAE, Report No.
978-92-9260-555-1
.https://www.irena.org/-/media/Files/IRENA/Agency/Publication/2023/Oct/COP28_IRENA_GRA_Tripling_renewables_doubling_efficiency_2023.pdf
2.
IRENA
,
2019
, “
IRENA (2019), Future of Wind: Deployment, Investment, Technology, Grid Integration and Socio-Economic Aspects (A Global Energy Transformation Paper)
,” International Renewable Energy Agency, Abu Dhabi, UAE, Report No.
978-92-9260-155-3
.https://www.irena.org/-/media/Files/IRENA/Agency/Publication/2019/Oct/IRENA_Future_of_wind_2019_summ_EN.pdf?la=en&hash=D07089441987EBABC7F4BED63B62C83820C18724
3.
Stevens
,
R. J.
, and
Meneveau
,
C.
,
2017
, “
Flow Structure and Turbulence in Wind Farms
,”
Annu. Rev. Fluid Mech.
,
49
(
1
), pp.
311
339
.10.1146/annurev-fluid-010816-060206
4.
Porté-Agel
,
F.
,
Bastankhah
,
M.
, and
Shamsoddin
,
S.
,
2020
, “
Wind-Turbine and Wind-Farm Flows: A Review
,”
Boundary-Layer Meteorol.
,
174
(
1
), pp.
1
59
.10.1007/s10546-019-00473-0
5.
Kleusberg
,
E.
,
Benard
,
S.
, and
Henningson
,
D. S.
,
2019
, “
Tip-Vortex Breakdown of Wind Turbines Subject to Shear
,”
Wind Energy
,
22
(
12
), pp.
1789
1799
.10.1002/we.2403
6.
Troldborg
,
N.
,
Sørensen
,
J. N.
, and
Mikkelsen
,
R.
,
2007
, “
Actuator Line Simulation of Wake of Wind Turbine Operating in Turbulent Inflow
,”
J. Phys.: Conf. Ser.
,
75
, p.
012063
.10.1088/1742-6596/75/1/012063
7.
Chamorro
,
L. P.
, and
Porté-Agel
,
F.
,
2010
, “
Effects of Thermal Stability and Incoming Boundary-Layer Flow Characteristics on Wind-Turbine Wakes: A Wind-Tunnel Study
,”
Boundary-Layer Meteorol.
,
136
(
3
), pp.
515
533
.10.1007/s10546-010-9512-1
8.
Medici
,
D.
, and
Alfredsson
,
P.
,
2006
, “
Measurements on a Wind Turbine Wake: 3D Effects and Bluff Body Vortex Shedding
,”
Wind Energy
,
9
(
3
), pp.
219
236
.10.1002/we.156
9.
Wu
,
Y.-T.
, and
Porté-Agel
,
F.
,
2012
, “
Atmospheric Turbulence Effects on Wind-Turbine Wakes: An LES Study
,”
Energies
,
5
(
12
), pp.
5340
5362
.10.3390/en5125340
10.
De Cillis
,
G.
,
Cherubini
,
S.
,
Semeraro
,
O.
,
Leonardi
,
S.
, and
De Palma
,
P.
,
2022
, “
The Influence of Incoming Turbulence on the Dynamic Modes of an NREL-5 MW Wind Turbine Wake
,”
Renewable Energy
,
183
, pp.
601
616
.10.1016/j.renene.2021.11.037
11.
Abkar
,
M.
, and
Porté-Agel
,
F.
,
2016
, “
Influence of the Coriolis Force on the Structure and Evolution of Wind Turbine Wakes
,”
Phys. Rev. Fluids
,
1
(
6
), p.
063701
.10.1103/PhysRevFluids.1.063701
12.
Gadde
,
S. N.
, and
Stevens
,
R. J.
,
2019
, “
Effect of Coriolis Force on a Wind Farm Wake
,”
J. Phys.: Conf. Ser.
,
1256
(
1
), p.
012026
.10.1088/1742-6596/1256/1/012026
13.
Abkar
,
M.
,
Sørensen
,
J. N.
, and
Porté-Agel
,
F.
,
2018
, “
An Analytical Model for the Effect of Vertical Wind Veer on Wind Turbine Wakes
,”
Energies
,
11
(
7
), p.
1838
.10.3390/en11071838
14.
Howland
,
M. F.
,
González
,
C. M.
,
Martínez
,
J. J. P.
,
Quesada
,
J. B.
,
Larranaga
,
F. P.
,
Yadav
,
N. K.
,
Chawla
,
J. S.
, and
Dabiri
,
J. O.
,
2020
, “
Influence of Atmospheric Conditions on the Power Production of Utility-Scale Wind Turbines in Yaw Misalignment
,”
J. Renewable Sustainable Energy
,
12
(
6
), p.
063307
.10.1063/5.0023746
15.
Narasimhan
,
G.
,
Gayme
,
D. F.
, and
Meneveau
,
C.
,
2022
, “
Effects of Wind Veer on a Yawed Wind Turbine Wake in Atmospheric Boundary Layer Flow
,”
Phys. Rev. Fluids
,
7
(
11
), p.
114609
.10.1103/PhysRevFluids.7.114609
16.
Jovanović
,
M. R.
,
Schmid
,
P. J.
, and
Nichols
,
J. W.
,
2014
, “
Sparsity-Promoting Dynamic Mode Decomposition
,”
Phys. Fluids
,
26
(
2
), p.
024103
.10.1063/1.4863670
17.
Rowley
,
C. W.
,
Mezić
,
I.
,
Bagheri
,
S.
,
Schlatter
,
P.
, and
Henningson
,
D. S.
,
2009
, “
Spectral Analysis of Nonlinear Flows
,”
J. Fluid Mech.
,
641
, pp.
115
127
.10.1017/S0022112009992059
18.
Schmid
,
P. J.
,
2010
, “
Dynamic Mode Decomposition of Numerical and Experimental Data
,”
J. Fluid Mech.
,
656
, pp.
5
28
.10.1017/S0022112010001217
19.
Debnath
,
M.
,
Santoni
,
C.
,
Leonardi
,
S.
, and
Iungo
,
G. V.
,
2017
, “
Towards Reduced Order Modelling for Predicting the Dynamics of Coherent Vorticity Structures Within Wind Turbine Wakes
,”
Philos. Trans. R. Soc. A
,
375
(
2091
), p.
20160108
.10.1098/rsta.2016.0108
20.
Kleine
,
V. G.
,
Kleusberg
,
E.
,
Hanifi
,
A.
, and
Henningson
,
D. S.
,
2019
, “
Tip-Vortex Instabilities of Two in-Line Wind Turbines
,”
J. Phys.: Conf. Ser.
,
1256
(
1
), p.
012015
.10.1088/1742-6596/1256/1/012015
21.
Sun
,
C.
,
Tian
,
T.
,
Zhu
,
X.
,
Hua
,
O.
, and
Du
,
Z.
,
2021
, “
Investigation of the Near Wake of a Horizontal-Axis Wind Turbine Model by Dynamic Mode Decomposition
,”
Energy
,
227
, p.
120418
.10.1016/j.energy.2021.120418
22.
Iungo
,
G. V.
,
Santoni-Ortiz
,
C.
,
Abkar
,
M.
,
Porté-Agel
,
F.
,
Rotea
,
M. A.
, and
Leonardi
,
S.
,
2015
, “
Data-Driven Reduced Order Model for Prediction of Wind Turbine Wakes
,”
J. Phys.: Conf. Ser.
,
625
, p.
012009
.10.1088/1742-6596/625/1/012009
23.
De Cillis
,
G.
,
Semeraro
,
O.
,
Leonardi
,
S.
,
De Palma
,
P.
, and
Cherubini
,
S.
,
2022
, “
Dynamic-Mode-Decomposition of the Wake of the NREL-5 MW Wind Turbine Impinged by a Laminar Inflow
,”
Renewable Energy
,
199
, pp.
1
10
.10.1016/j.renene.2022.08.113
24.
Sescu
,
A.
, and
Meneveau
,
C.
,
2014
, “
A Control Algorithm for Statistically Stationary Large-Eddy Simulations of Thermally Stratified Boundary Layers
,”
Q. J. R. Meteorol. Soc.
,
140
(
683
), pp.
2017
2022
.10.1002/qj.2266
25.
Smagorinsky
,
J.
,
1963
, “
General Circulation Experiments With the Primitive Equations: I. The Basic Experiment
,”
Mon. Weather Rev.
,
91
(
3
), pp.
99
164
.10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2
26.
Santoni
,
C.
,
Ciri
,
U.
,
Rotea
,
M.
, and
Leonardi
,
S.
,
2015
, “
Development of a High Fidelity CFD Code for Wind Farm Control
,” 2015 American Control Conference (
ACC
), Chicago, IL, July 1–3, pp.
1715
1720
.10.1109/ACC.2015.7170980
27.
Santoni
,
C.
,
Carrasquillo
,
K.
,
Arenas-Navarro
,
I.
, and
Leonardi
,
S.
,
2017
, “
Effect of Tower and Nacelle on the Flow Past a Wind Turbine
,”
Wind Energy
,
20
(
12
), pp.
1927
1939
.10.1002/we.2130
28.
Martinez
,
L.
,
Leonardi
,
S.
,
Churchfield
,
M.
, and
Moriarty
,
P.
,
2012
, “
A Comparison of Actuator Disk and Actuator Line Wind Turbine Models and Best Practices for Their Use
,”
AIAA
Paper No. 2012-0900.10.2514/6.2012-0900
29.
Sørensen
,
J. N.
, and
Shen
,
W. Z.
,
1999
, “
Computation of Wind Turbine Wakes Using Combined Navier-Stokes/Actuator-Line Methodology
,”
1999 European Wind Energy Conference
, Nice, France, Mar. 1–5, pp.
156
159
.
30.
Orlandi
,
P.
, and
Leonardi
,
S.
,
2006
, “
DNS of Turbulent Channel Flow With Two-and Three-Dimensional Roughness
,”
J. Turbul.
, 7(
1
), p.
N73
.10.1080/14685240600827526
31.
De Cillis
,
G.
,
Cherubini
,
S.
,
Semeraro
,
O.
,
Leonardi
,
S.
, and
De Palma
,
P.
,
2021
, “
POD–Based Analysis of a Wind Turbine Wake Under the Influence of Tower and Nacelle
,”
Wind Energy
,
24
(
6
), pp.
609
633
.10.1002/we.2592
32.
Bernardoni
,
F.
,
Ciri
,
U.
,
Rotea
,
M.
, and
Leonardi
,
S.
,
2020
, “
Real-Time Identification of Clusters of Turbines
,”
J. Phys.: Conf. Ser.
,
1618
(
2
), p.
022032
.10.1088/1742-6596/1618/2/022032
33.
Flay
,
R.
, and
Stevenson
,
D.
,
1988
, “
Integral Length Scales in Strong Winds Below 20 m
,”
J. Wind Eng. Ind. Aerodyn.
,
28
(
1–3
), pp.
21
30
.10.1016/0167-6105(88)90098-0
34.
Stanislawski
,
B. J.
,
Thedin
,
R.
,
Sharma
,
A.
,
Branlard
,
E.
,
Vijayakumar
,
G.
, and
Sprague
,
M. A.
,
2023
, “
Effect of the Integral Length Scales of Turbulent Inflows on Wind Turbine Loads
,”
Renewable Energy
,
217
, p.
119218
.10.1016/j.renene.2023.119218
35.
Jonkman
,
J.
,
Butterfield
,
S.
,
Musial
,
W.
, and
Scott
,
G.
,
2009
, “
Definition of a 5-MW Reference Wind Turbine for Offshore System Development
,” National Renewable Energy Laboratory, Golden, CO, Report No.
NREL/TP-500-38060
.https://www.nrel.gov/docs/fy09osti/38060.pdf
36.
Orlanski
,
I.
,
1976
, “
A Simple Boundary Condition for Unbounded Hyperbolic Flows
,”
J. Comput. Phys.
,
21
(
3
), pp.
251
269
.10.1016/0021-9991(76)90023-1
37.
de Tullio
,
M. D.
,
De Palma
,
P.
,
Iaccarino
,
G.
,
Pascazio
,
G.
, and
Napolitano
,
M.
,
2007
, “
An Immersed Boundary Method for Compressible Flows Using Local Grid Refinement
,”
J. Comput. Phys.
,
225
(
2
), pp.
2098
2117
.10.1016/j.jcp.2007.03.008
38.
Abkar
,
M.
, and
Porté-Agel
,
F.
,
2015
, “
Influence of Atmospheric Stability on Wind-Turbine Wakes: A Large Eddy Simulation Study
,”
Phys. Fluids
,
27
(
3
), p.
035104
.10.1063/1.4913695
39.
Mao
,
X.
, and
Sørensen
,
J. N.
,
2018
, “
Far-Wake Meandering Induced by Atmospheric Eddies in Flow Past a Wind Turbine
,”
J. Fluid Mech.
,
846
, pp.
190
209
.10.1017/jfm.2018.275
40.
Gupta
,
V.
, and
Wan
,
M.
,
2019
, “
Low-Order Modelling of Wake Meandering Behind Turbines
,”
J. Fluid Mech.
,
877
, pp.
534
560
.10.1017/jfm.2019.619