Low-order modelling of wake meandering behind turbines

Based on recent numerical simulations and field experiments, the mechanism behind wake meandering is increasingly accepted to be through the amplification of upstream disturbances owing to the convectively unstable nature of the flow. In this paper, we deduce a low-order phenomenological model for the far-wake region, which is based on a modified form of the complex Ginzburg–Landau (CGL) equation for flows that are in the amplifier regime, i.e. are only convectively unstable. The model reproduces the main qualitative features of wake meandering: (i) its origin via amplification of upstream structures, (ii) dependence of oscillation frequency on the upstream disturbance amplitude (higher amplitudes lead to lower frequencies), (iii) shift towards lower frequencies as the wake flow evolves in the streamwise direction and, to an extent, (iv) the transfer of energy from very low frequencies towards relatively higher frequencies. Additionally, the model also predicts the increase in the meandering amplitude and an advancement in its onset with increasing thrust coefficient. To our knowledge, this is the first low-order dynamical system in the literature that models wake meandering. The model coefficients are obtained from the mean flow local stability results that we show correctly account for the changing operating conditions and thus pave way for the prediction of wake meandering features. Its low order makes it suitable to use inside an energy farm design model, where it can help to mitigate the adverse effects of wake meandering.


Introduction
Far-wake regions behind tidal or wind turbines are usually reported to have low-frequency oscillations (Medici & Alfredsson 2006;Larsen et al. 2008;Chamorro et al. 2013).These oscillations give rise to a meandering wake pattern, hence the term wake meandering, and cause an increase in turbulence in the far-wake flow.In an energy farm, this increased turbulence level adversely affects the performance and load characteristics of the downstream turbines (Ainslie 1988;Vermeer, Sørensen & Crespo 2003;Larsen et al. 2008).Qualitative understandings of the origin of wake 536 V. Gupta and M. Wan predict them under different operating conditions can enable an energy farm designer to account for the adverse effects of wake meandering.
The existing models for turbine wake flows are kinematic in nature, such as the stochastic wake meandering model by Larsen et al. (2007) or the input-output correlation model by Hamilton et al. (2018).Such models are useful in different ways, the model of Larsen et al. (2007) is a comprehensive farm model that includes the effect of wake meandering and the model of Hamilton et al. (2018) provides quantitative data for wake evolution under complex inflow conditions.Kinematic nature of these models, however, means they do not account for the underlying mechanism and hence cannot explain the origin of the above mentioned features.Towards that purpose, we aim to develop a low-order model for the far-wake region that qualitatively reproduces the phenomenon of wake meandering.
We obtain wake flows behind a tidal turbine under uniform and sinusoidally varying inflow conditions by performing large-eddy simulations (LES) in § 2. The obtained wake flows exhibit the main qualitative features of wake meandering that are reported in the literature.We then perform a local linear stability analysis in § 3 to characterise the mean flow and use those results in the low-order model deduced in § 4. We discuss the achievements and limitations of the model and physical insights gained through it in § 5.

Methodology
We obtain the flow fields behind a turbine placed in a straight channel with a rectangular cross-section using LES, where the scales larger than the grid size are directly resolved by solving the spatially filtered Navier-Stokes equations given below.
where u i represents the filtered velocity ( u, v, w) in the (x, y, z)-directions (of the Cartesian coordinates) for i = 1, 2 and 3, respectively.The modified pressure is given as p * = p/ρ + 0.5 u i 2 , where p is the filtered pressure and ρ is the density.The sub-grid scales (SGS) are modelled as SGS stress terms ( τ ij ) according to Meneveau, Lund & Cabot (1996).The viscous term is neglected.An external time-dependent forcing (F x ) maintains a spatially uniform inlet flow velocity in the x-direction that is either constant or varying sinusoidally with time.Lastly, F T i represents the force imparted by the turbine on the flow.
The turbine geometry is based on the National Renewable Energy Laboratory's hypothetical 550 kW two-bladed machine that is also simulated in Churchfield, Li & Moriarty (2013).It has a rotor diameter (D) of 20 m, the blade sections are of NACA 63(1)-424 airfoil shape and the blade chord varies from approximately 1.7 m at the root to 0.6 m at the tip.We do not resolve the turbine geometry, instead we assume the blades as actuator lines.These lines are divided into a number of segments (N T ) and each segment imparts an aerodynamic force ( f T j i ) on the flow.This force is based on the pre-tabulated lift and drag coefficients of the blade geometry (Sørensen & Shen 2002).They are imposed on a cell located at (x, y, z) as a three-dimensional Gaussian distribution around the centre of each segment (x j , y j , z j ).The effect of the turbine on the flow (F T i ) is then calculated as i (x j , y j , z j ) Low-order modelling of wake meandering 537 where r j is the distance between (x j , y j , z j ) and (x, y, z), and is the Gaussian distribution parameter.We account for the tip loss through a pre-multiplication factor f g ≡ (2/π) cos −1 [exp(−(0.5D− d j /d j sin β j ))], where d j is the radial distance of the jth blade section from the blade root and β j is the angle between the local relative velocity and the rotor plane (Shen et al. 2005).Additionally, F T i also includes the effect of nacelle.Following Wu & Porté-agel (2011), we model it as a porous disc of diameter 0.1D and drag coefficient 1.2.We follow the numerical methodology from Churchfield et al. (2013) closely, and borrow the actuator lines code from NREL's software SOWFA (Churchfield et al. 2013).It is based on an open-source finite volume solver OpenFOAM (Jasak 1996).
The channel simulated here is 12.0D long in the streamwise direction (x) and 2.5D wide in the other two directions.The turbine is placed 2.0D downstream of the inlet and at the centre of the y-z plane.The location of its centre is called (0, 0, 0).The grid is a uniform hexahedral mesh with 480 × 100 × 100 cells, which means the cell size is 1/40D in each direction.The side, upper and lower walls are modelled as free-slip boundary conditions.Each turbine blade is divided into N T = 40 segments.The distribution parameter is set as two times the grid size (i.e.= 1/20D).The time step is set as t = 0.025 s in all turbine simulations.This keeps the maximum Courant number below 0.2 and does not allow the actuator lines to travel across more than a cell in a time step.Troldborg (2009) recommended these limits for and t for avoiding numerical instabilities and are widely followed in the literature.The power spectral density (PSD) of the w-velocity fluctuations (φ w ) is obtained using data from t = 200 s to 2200 s, unless stated otherwise, sampled at every 0.1 s, and by using the 'pwelch' command in MATLAB.It is not normalised.

Simulation results under uniform inflow conditions
The inlet flow is maintained as (U, 0, 0), where U = 1.85 m s −1 .The turbine is operated at the rotor frequency of 0.15-0.20 Hz, the corresponding Ω = 9.0 to 12.0 expresses it in terms of revolutions per minute.The flow corresponding to Ω = 10.5 r.p.m. is treated as the representative case for which figure 1(a) shows the instantaneous axial velocity field in the x-z plane passing through the turbine's centre.It should be noted that owing to the uniform laminar inflow conditions, the wake flow evolution is slower here as compared to in cases of turbulent inflow conditions, such as in Churchfield et al. (2013).
In the near-wake region (up to x/D ≈ 4.5), the velocity deficit is maximum in the core region (z/D ≈ 0.325) and is negligible in the area between the blade root and the nacelle.The tip and root vortices are present behind the blades in the immediate downstream region.Other vortices, called near-wake vortices here, also start to appear in the downstream region as indicated in the figure.The PSD in the near-wake region (panel b) shows peaks at f t = 0.350 Hz (root and tip vortices, f t is twice the rotor frequency) and f n ≈ 0.219 Hz (near-wake vortices).Qualitatively, the near-wake region matches well with the LES results in Kumar & Mahesh (2017) where a propeller's geometry is fully resolved.We do not scrutinise the origin of near-wake structures here, which is a complex area of research as evident from the numerous studies (Widnall 1972;Okulov & Sørensen 2007;Felli, Camussi & Di Felice 2011;Hong et al. 2014).In this paper, we are concerned only with the effect of the obtained near-wake structures on the far-wake region.The flow evolves and becomes turbulent in the far-wake region where a meandering pattern is observed (indicated in the figure).The PSDs in the far-wake region (panel c) show an increase in energy at all frequencies, which is due to the turbulence.More importantly, there is a peak at a low frequency ( f m ≈ 0.048 Hz) that corresponds to the observed wake meandering.It should be noted that this peak is broad, i.e. it comprises a band of frequencies as also shown in panel (d) and figure 2, which hints that it is a result of convective instability (i.e. via the amplification of the upstream disturbances) (Huerre & Monkewitz 1990).In order to trace the origin of the observed wake meandering, we plot φ w at different streamwise locations in panel (d).In the later part of the near-wake region (at x/D = 2.5), there are peaks at f ≈ 0.315 Hz, 0.219 Hz and their linear combinations.Successive φ w at x/D = 4.0, 6.0 and 9.0 show the peak at 0.048 Hz gets amplified and this, we speculate, gives rise to the observed wake meandering in the far-wake region.When the frequency ( f ) is nondimensionalised as St = f D/U, St m ≈ 0.52 is found to be in the range observed in Mao & Sørensen (2018).
The origin of wake meandering via the amplification of near-wake structures is consistent with (i) the new mechanism discussed in § 1, (ii) convectively unstable nature of the flow shown in § 3 and (iii) the evolution of far-wake vortical structures in flows behind cylinders.Williamson & Prasad (1993) showed that in the far-wake region behind a cylinder (which is much further downstream than a turbine's far-wake region), the flow amplifies a subharmonic component of the near-wake vortex shedding structure (particularly in the absence of other upstream perturbations) that gives rise to sustained vortical oscillations in the far-wake region.The results presented here show how near-wake structures can affect the far-wake region in flows behind turbines.This provides a plausible explanation for the hypotheses in Okulov & Sørensen (2007), Iungo et al. (2013) and Viola et al. (2014), where they suggested that instabilities in the near-wake region can lead (or contribute) to meandering in the far-wake region.
Figure 2(a,c,e,g) shows the non-dimensional instantaneous axial flow fields at increasing Ω from the top to bottom.The flow fields show some level of meandering in all cases.It is relatively weak in the Ω = 9.0 r.p.m. case and gets stronger with increasing Ω.This is consistent with the observations in Heisel et al. (2018), where wake meandering is not observed for highly derated turbines, i.e. when the thrust coefficient becomes too small.The corresponding PSDs at z/D = 0.325 at four streamwise locations (x/D = {2.5, 4.0, 6.0, 9.0}) are presented in the adjacent plots (b,d,f,h).Except in the Ω = 9.0 r.p.m. case, there is a broad peak in φ w at St ≈ 0.5.The fact that St is nearly constant in all the cases suggests that the wake meandering frequency scales with the turbine diameter and the mean incoming velocity, as observed in Foti et al. (2018b) and others.The weak dependence of the peak frequency on Ω shows the effect of the near-wake structures.This effect is expected to be small under turbulent inflow conditions.540 V. Gupta and M. Wan x/D 8 10 2.3.Simulation results under sinusoidally varying inflow conditions Turbines usually operate under turbulent inflow conditions where energetic structures from the boundary layer, upstream turbines, surface waves and other intermittent sources are present.Therefore, it is imperative to study the wake evolution under perturbed inflow conditions.In this paper, we limit ourselves to perturbations in the axial velocity that are spatially uniform and vary sinusoidally in time as (U(1 + A f sin(2πSt f (U/D)t)), 0, 0), where U is the same as before, A f and St f are the non-dimensional forcing amplitude and frequency, respectively.Although such a forcing is a poor approximation of turbulent inflow conditions for utility-scale turbines, many studies have shown harmonic response of a flow can provide valuable insights into instability generated flow structures.These include McKeon & Sharma (2010) for understanding the scaling of coherent structures in turbulent pipe flows and Garnaud et al. (2013) for understanding the preferred mode selection in turbulent axisymmetric incompressible jets.
We obtain the wake flow for the turbine operating at Ω = 10.5 r.p.m. under various St f and A f inflow conditions.Figure 3 presents the non-dimensional instantaneous axial velocity fields in the x-z plane.The four columns correspond to four forcing frequencies, increasing from the left to right, while the forcing amplitudes (increasing from the bottom to top) are mentioned in the respective panels.In the absence of sinusoidal forcing, the far-wake oscillations result from the amplification of background noise (which includes the decomposing near-wake structures and the resulting turbulence).Consequently, the wake meandering pattern is irregular (figure 1 and the lowest row here) and its spectrum has a broad peak at St m (see figure 1d).As A f is increased, the wake meandering pattern becomes regular (the top row except in the St f = 0.11 case) and the flow starts to exhibit periodic oscillations at St f (see figure 5).This phenomenon is similar to the lock-in observed in oscillators (Pikovsky, Rosenblum & Kurths 2001), where an oscillator starts to oscillate at the forcing frequency (provided the forcing amplitude is sufficiently high) in place of its natural frequency.We call it pseudo lock-in here because our flow is not an oscillator system, it is an amplifier of external disturbances.At pseudo lock-in, all other frequencies are suppressed (i.e. the given background noise is no longer amplified) and the wake flow starts to oscillate at the forcing frequency similar to a globally unstable flow.
We also see from figure 3 that the wake meandering pattern becomes regular at lower A f values when St f = 0.43 and 0.76 as compared to when St f = 0.11 and 0.97.Thus showing that pseudo lock-in is achieved earlier when St f is closer to St m .The variation in A f that is required to achieve pseudo lock-in with St f is shown in figure 4 for the Ω = 10.5 r.p.m. case.The colour map represents the ratio , where φ w (St f ) is the value of φ w at St f (representing the wake flow response to the sinusoidal forcing) and φ w [St m ] is the value of φ w at the broadband peak near St m (representing the wake flow response to the background noise) calculated at (x/D, z/D) = (9, 0.325).The black line represents the minimum A f required to achieve pseudo lock-in (r f ≈ 0.95).The pseudo lock-in curve has a parabola-like shape with a minimum at St f ≈ 0.54-0.65,which is close to St m ≈ 0.52, and a little tilt towards St f > St m .The hollow circles indicate the locations where the calculations are performed, φ w at many points in figure 4 are calculated using only 400 s of data.The purpose here is to show the trend of pseudo lock-in, for which the accuracy of φ w is not required.
Quantitative response of the wake flow to sinusoidal forcing is calculated in terms of where (|u|, |v|, |w|) are frequency spectra of the (ũ, ṽ, w) fluctuations, respectively.These calculations are performed using 1000 s of data sampled at 0.1 s. Figure 5 shows E f at (x/D, z/D) = (9, 0.325) for the Ω = 10.5 r.p.m. case, where the three columns correspond to three forcing frequencies (increasing from the left to right) and the three rows correspond to three forcing amplitudes (increasing from the bottom to top).A shift towards lower frequencies with increasing forcing amplitudes can be seen here.At the smallest forcing amplitude (bottom row), the maximum response is at the highest forcing frequency (i.e.St f = 0.76).At the medium and the largest forcing amplitudes, the maximum response shifts to St f = 0.54 and 0.32, respectively.This is in agreement with the findings of Mao & Sørensen (2018).
V. Gupta and M. Wan

Local linear stability analysis
Our purpose here is to find local stability characteristics of the turbine wake flows.These results guide the development of the low-order model in § 4. We perform stability analysis on the time-averaged flow profiles, which are not the stationary solutions of the Navier-Stokes equations.Nonetheless, such analysis is found to be effective in finding shear layer generated flow oscillations in the literature (Garnaud et al. 2013).

Local mean velocity profiles
We assume the mean flow to be axisymmetric around the (y, z) = 0 axis and to be varying slowly in the x-direction for the Wentzel-Kramers-Brillouin-Jeffrey (WKBJ) approximation to be valid.The mean axial (U) and azimuthal (W) velocities are then obtained by time averaging ũ(x, 0, z, t) and ṽ(x, 0, z, t), respectively, and are nondimensionalised by the mean incoming velocity.Figure 6(a) presents U − 1 (blue lines) and W (red lines) at various streamwise locations for the Ω = 10.5 r.p.m. case.The velocity scales are shown at the top of the horizontal axis.The background colours here indicate the division between the near-and far-wake regions.In the near-wake region (x/D = 0.0-4.5), the mean axial velocity profiles show two regions of deficits similar to Foti et al. (2018a) -the outer one due to the turbine blades and a smaller middle one due to the nacelle.Similarly, the azimuthal velocity has two regions as well -behind the blades it is in the opposite direction to the turbine rotation, while behind the nacelle it is smaller and in the same direction as the turbine rotation.The near-wake region is followed by the transition (x/D = 4.5-5.5),far-wake (x/D = 5.5-9.5) and buffer regions (x/D = 9.5-10.0).The mean velocity profiles in these regions do not have the middle wake and the shear layers are not as sharp as in the near-wake region.

Methodology
Because the mean flow is axisymmetric, the perturbation equations are written in (x, r, θ ) as axial, radial and azimuthal coordinates.The perturbations are assumed to be of the form (u x , iu r , u θ , p) exp(i(mθ + kx − ωt)), where u x , u r , u θ and p are perturbations to the (x, r, θ) velocities and pressure, respectively, m is the azimuthal wavenumber, k = (k r + ik i ) is the streamwise wavenumber and ω = (ω r + iω i ) is the angular frequency (subscripts r and i stand for real and imaginary parts, respectively).The linearised perturbations equations are derived from the Navier-Stokes equations as ) The stability analysis code is based on the Chebyshev spectral collocation method.The discretisation in the radial direction is performed on a Gauss-Lobatto-Chebyshev grid.It is mapped on a radially unbounded physical space r/D , where ζ is the Chebyshev grid from 0 to 1 and R max = 25 is an arbitrarily large number to represent an unbounded space.The number of collocation points used is N = 120.The change in the results on doubling R max and N (not shown here) is less than a per cent.

Stability analysis results
We first heuristically define the concepts of convective and absolute instabilities and their relation to the global flow instability (for detailed reviews please see Huerre & Monkewitz (1990) and Chomaz (2005)).A parallel flow is linearly stable if ω i <0 for all values of k, otherwise it is linearly unstable.In a linearly stable parallel flow, all external perturbations eventually die everywhere in the domain, whereas in a linearly unstable parallel flow, there are two possibilities.The first is the convective instability where perturbations of certain wavenumbers can grow but the growth rate is slower than the advection rate, i.e. perturbations get advected away from their origin.Mathematically, this happens when ω i > 0 for some values of k, but the modes with zero group velocity (i.e.∂ω/∂k = 0) have ω i < 0. Such flows do not exhibit self-sustained oscillations but can greatly amplify external perturbations and are called amplifier flows.The second is the absolute instability where perturbations of certain wavenumbers can grow and the growth rate is faster than the advection rate (Tobias, Proctor & Knobloch 1998) (i.e. at least one mode with ∂ω/∂k = 0 has ω i > 0).Such flows exhibit self-sustained oscillations and are called oscillator flows.
These concepts of local stability analysis are applicable in weakly non-parallel flows where the basic flow varies on a longer length scale as compared to the instability wavelength (Monkewitz, Huerre & Chomaz 1993;Chomaz 2005).A weakly non-parallel flow is (i) stable -when it is locally linearly stable everywhere, (ii) an amplifier of external perturbations -when it is locally convectively unstable in some part of the domain, and (iii) an oscillator -when it is locally absolutely unstable in a sufficiently large region of the domain.A transition from stable to oscillator flow with increasing instability is shown in § 4.2.
Turbine wake flows are only locally convectively unstable (Iungo et al. 2013;Mao & Sørensen 2018) and thus belong to the second class (i.e.amplifier flow).In figure 6(b-d), therefore, we present the stability results corresponding to the locally most amplified modes (unless explicitly stated, the results are for Ω = 10.5 r.p.m. case and m = −1 mode).Panel (b) shows the variation of the most amplified frequency (ω r /2π) with x/D.The near-wake region of the flow amplifies higher frequencies (compared to the far-wake region) because the mean flow in this region has smaller-scale features (such as the middle wake) and sharper shear layers.The sudden changes in the frequency signify mode switchings (determined by the mean flow feature that is dominant locally).Again, we do not scrutinise the origin of different near-wake modes.In the transition region, there is a clear shift from the higher to lower frequencies.Finally, in the far-wake region, there is a slow but consistent decrease in the most amplified frequency.Panels (c) and (d) show the corresponding spatial growth rate (k i ), the group velocity (∂ω/∂k| r ), and the diffusion term (−10(∂ 2 ω/∂k 2 )| i ).The spatial growth rate has the same behaviour as ω sr , it decreases slowly in the far-wake region as the flow evolves and the shear layers become less sharp.The group velocity signifies the advection rate of the perturbations and is roughly equal to 0.75 times the mean incoming velocity.The negative value of the diffusion term signifies that very high-frequency perturbations decay in the flow.The imaginary part of the group velocity (∂ω/∂k| i ) at (ω s , k s ) is zero by definition, the linear dispersion term (∂ 2 ω/∂k 2 | r ) is negligibly small and the real wavenumber is given as k sr ≈ ω sr (∂ω/∂k) −1 (please see appendix A).Low-order modelling of wake meandering 545 Additionally, in panel (a), we show the normalised eigenfunction components at locations x/D = 6, 7, 8 and 9 for the m = 0 and −1 modes in terms of u 2 x + u 2 r + u 2 θ from r/D = 0 to 1.2.Both the modes initially show two peaks (at r/D ≈ 0.2 and 0.5), but later at x/D = 9.0, only the outer peak at r/D ≈ 0.5 survives.The eigenfunctions show the flow region in the radial direction where the wake meandering is influential (see (4.5)).In panel (b), we show the variation of the most amplified frequency in the Ω = 12.0 r.p.m. and 9.0 r.p.m. cases (only after their respective transition regions).
The results for these flow cases are qualitatively similar, implications of this similarity are explored in § 5.1.In panel (c), we show the spatial growth rates for the m = 1, 0 and −2 modes, which indicate that the local stability results are very similar for all the azimuthal modes (see § 4.2).

Low-order modelling of the far-wake region
For spatially developing open shear flows, such as jets and wakes, Chomaz and co-workers (Chomaz et al. 1988;Chomaz, Huerre & Redekopp 1990, 1991;Chomaz 1992;Le Dizès et al. 1996;Cossu & Chomaz 1997) pioneered the use of the complex Ginzburg-Landau (CGL) equation given as This equation governs the evolution of a hydrodynamic instability wave (in terms of the complex amplitude A = |A|e iφ , where φ is its phase) travelling downstream at the group velocity U g .The growth rate σ r is the driving term, σ i is the frequency shift, c dr (> 0) is the diffusive coupling coefficient, c nr (> 0) is the nonlinear saturation coefficient and c di and c ni are the linear and nonlinear dispersion coefficients, respectively.
4.1.Model deduction Strictly, the CGL equation is limited to describing instability waves in flows that are (i) weakly nonlinear and (ii) weakly non-parallel (Aranson & Kramer 2002;van Saarloos 2003).The first condition is satisfied when a flow is only marginally unstable (i.e. has just transited to the oscillator regime) and ensures the absence of the higher harmonics.The second condition ensures that instability waves satisfy the local dispersion relations as per the WKBJ approximation (Monkewitz et al. 1993).Thus, it provides a way to obtain the linearised CGL coefficients from the local stability results (Le Dizès et al. 1996) as where X ≡ x ( 1) is the slow scale at which the basic flow varies, ω 0 (X) and k 0 (X) are the complex frequency and wavenumber, respectively, at X corresponding to the zero group velocity modes (i.e.∂ω/∂k(ω 0 , k 0 ; X) = 0) and Ã is the linear counterpart of A. Following Crighton & Gaster (1976), Monkewitz et al. (1993) and Siconolfi et al. (2017), linear self-sustained oscillations in a weakly non-parallel flow can be approximated as where W g (X) is the slow amplitude variation, q g (r, θ ; X) and k g (X) are the local eigenfunction and wavenumber, respectively, corresponding to the complex global mode frequency ω g .In order to check the applicability of (4.2), we insert (4.3) in (4.2) and retain only the leading-order terms (at O( 0)).The resultant equation is the local linear dispersion relation at X up to the second-order term as Although, in general, flows have higher-order terms in their dispersion relations, they can be neglected because in a weakly non-parallel flow ω g is close to ω 0 (X) (Chomaz et al. 1991;Monkewitz et al. 1993;Le Dizès et al. 1996;Pier & Huerre 2001;Pier 2002).
In contrast to self-sustained oscillations modelled by (4.2), wake meandering arises from amplification of upstream disturbances (see § § 1 and 2).We argue that (i) the CGL equation can still be used to model wake meandering and (ii) its coefficients can still be obtained from local stability results.In support of the first argument, Cossu & Chomaz (1997), Tobias et al. (1998) and Bagheri et al. (2009) showed that flow structures arising from the amplification of upstream perturbations can be qualitatively modelled using the CGL equation.In support of the second argument, Gaster (1969) showed that spatially amplified modes can be modelled using a wave equation tuned by local mean flow parameters.More formally, Chomaz (2005) and Viola, Arratia & Gallaire (2016) showed that similar to (4.3), the linear response to sinusoidal forcing (say at frequency ω f ) in weakly non-parallel flows can be approximated as where W f is the slow amplitude variation, q f (r, θ; X) and k f (X) are the local eigenfunction and wavenumber, respectively, corresponding to ω f .The range of ω f that is of interest is where the spatial amplification is maximum, which is at ω s (X) (see § 3.3).The CGL equation to describe the spatially amplified waves, therefore, should be based on the dispersion relation around (ω s , k s ) (see appendix A) and can be written as It should be noted here that we replaced Ã by its nonlinear counterpart A and also included the nonlinear term (C n ≡ c nr + ic ni ).The conversion from the linear to nonlinear CGL equation in not trivial and is discussed in detail in § 5.2.Equation (4.6) is solved in the range X = 0-12, where X is equivalent to x/D in the wake flow.At the inlet (X = 0), the boundary condition is set as time-dependent forcing, which can be an impulse, sinusoidal or white noise forcing or their combination.At the outlet (X = 12), following Heaton, Nichols & Schmid (2009), the boundary condition is set as a convective outflow.We convert the partial differential equation to a set of ordinary differential equations by discretising the spatial coordinate using the fourth-order accurate central differencing scheme (dX = 1/32).We then use the implicit Crank-Nicolson scheme (dt = 0.04) for time marching and Picard's method for iterating the nonlinear term.
In the near-wake region (X < 4.5), the coefficients are extrapolated (as shown in the figure) such that the model response and turbine wake flow results roughly match in the region X = 2.0-3.0 (see § 5.1).In principle, the model coefficients should be fixed separately for different azimuthal modes and thus different CGL equations should be used for the forced response based on the azimuthal wavenumber of the inflow disturbances.However, the stability results are very similar at all m (as shown in figure 6c).We, therefore, assume that a single CGL equation should be able to qualitatively model the evolution of all m modes.Figure 8 shows the linear impulse response of the model (panel b) as well as for its variants (panels a, cand d).An impulse of magnitude 1 is applied at X = 0 at time t = 0, and the response amplitudes at X = {1, 2, 3, . . ., 10} with time are plotted.Panels (a-d) show transition from (a) the stable regime to (b) the amplifier regime to (c) being marginally close to the oscillator regime and finally to (d) the oscillator regime as the flow instability increases.In the stable regime in (a), a disturbance dies down monotonically.In the amplifier regime in (b), a disturbance shows a transient growth via a convective-type non-normality (Cossu & Chomaz 1997) before eventually decaying.The flow is still in the amplifier regime but is marginally close to the oscillator regime in (c).The transient growth in this case is much higher and the decay rate is much slower as compared to those in (b).Thus, even small external noise can sustain oscillations in such a flow (Huerre & Monkewitz 1990;Babcock, Ahlers & Cannell 1991).In the oscillator regime in (d), a disturbance grows exponentially until the nonlinearities saturate the oscillations to a finite amplitude.This figure confirms that, like the wake flow, the model (i) is in the amplifier regime and (ii) can spatially amplify the incoming disturbances.
The nonlinear saturation term (c nr ) represents the reduction in amplification of disturbances with their amplitude.It is observed in § 2.3, that this reduction is higher at the higher forcing frequencies (St f ).Thus, c nr should be proportional to a positive exponent of St f .Among the exponents tried (0.5, 1.0, 1.5 and 2.0), 1.5 gives the best results.Thus leading to c nr = 4(2πSt f ) 1.5 , where the pre-multiplication factor controls the level of the saturation.The nonlinear dispersion term (c ni ) represents the change in the wavenumber with the forcing amplitude and, similar to the linear dispersion term, is set as zero.When a number of frequencies are present in the boundary forcing, the nonlinear saturation term is set as c nr (X) = 4(2π(Σ j A fj St fj G Lj (X)/Σ j A fj G Lj (X))) 1.5 , where A fj and G Lj (X) are the forcing amplitude and the linear gain corresponding to the forcing frequency St fj (linear gains at some frequencies are shown in appendix B).The implications of the nonlinear coefficient and alternative ways for its mathematically rigorous calculation are discussed in § 5.2.

Model behaviour
Turbines operate under conditions comprising of small-scale turbulence as well as coherent structures arising from a number of sources.Therefore, we study the model behaviour when inflow conditions (at X = 0) are set as combinations of white Gaussian noise (representing small-scale turbulence) and sinusoidal forcing (representing coherent structures).Although simplistic, such input-output analyses are shown to be useful in understanding the origin of complex structures arising in turbulent flows (see § 2.3).
Figure 9 shows the model response to white Gaussian noise forcing (of normalised amplitude A w = 0.0004) in terms of the PSDs at X = {1, 2, 3, . . ., 10}.The PSDs evolve towards lower frequencies and have a broad peak at St ≈ 0.5 in the downstream region, which matches well with the wake flow PSDs shown in figure 1(d).The peak is broader here because the background noise is white, unlike in turbine far wake where it consists of the near-wake structures.In addition to the white noise, which is present in all the subsequent simulations, we add an upstream sinusoidal forcing term (of normalised amplitude A f and frequency St f ) at X = 0. Figure 10 shows the model responses to forcing at varying St f and A f in terms of where φ A (St f ) and φ A (St max ) are the PSD values (at X = 9) at the forcing frequency (St f ) and at the frequency (St max ) where the response is maximum but St max = St f .The ratio r fm is nearly zero when the background noise is dominant and approaches one as the sinusoidal forcing suppresses the effect of the background noise.The black line shows the pseudo lock-in curve (r fm ≈ 0.95), i.e. the minimum A f at a given St f where the effect of the background noise is suppressed.Similar to that for the wake flow in figure 4, pseudo lock-in is achieved at lower A f when St f ≈ 0.54-0.65,which is close to St m ≈ 0.52, and the pseudo lock-in curve is parabola like in shape with a tilt towards Figure 11 presents the model responses (in terms of the frequency spectrum at X = 9) to sinusoidal forcing at varying St f (increasing from the left to right) and A f (increasing from the bottom to top).The maximum response shifts from St f = 0.54 to 0.32 as the forcing amplitude increases.This shows a shift towards lower frequencies  The wake meandering frequency shifts to lower values with increasing (i) x/D and (ii) A f .with increasing A f , as also seen for the wake flow in figure 5 and Mao & Sørensen (2018).

Achievements of the low-order model
The match between figures 10 and 11 with figures 4 and 5, respectively, shows how well the CGL equation models the convectively amplified flow structures that give rise to wake meandering.Thus, justifying the use of the CGL equation (4.6) for modelling wake meandering.Here, we present a more detailed comparison between the model and the wake flow results.It should be noted that we only qualitatively model the wake from x/D = 4.5 to 9.5.The quantitative match between the two arises because the model is designed to roughly match the wake flow in the region x/D = 2.0 to 3.0.The model output |A| is equivalent to the wake flow response amplitude calculated as the summation of E f at St f and its harmonics (see figure 5) and then averaged over z/D = 0.0-0.5.

Variations in the wake meandering features with the forcing parameters
Figure 12 shows the variations of the wake flow and the model response amplitudes with the streamwise distance from the turbine.The results are presented for three forcing frequencies (increasing from left to right) and two forcing amplitudes (a,c,e-A f = 0.016, b,d,f -A f = 0.135).There are two main points to note from this figure.The first is the change in the dominant frequency with the downstream distance.The top row shows that the response at St f = 0.32 increases in the streamwise direction up to x/D ≈ 10, while the responses at St f = 0.54 and St f = 0.76 initially increase up to x/D ≈ 9 and x/D ≈ 7, respectively, and then decrease.Consequently, up to x/D ≈ 6, the wake flow response at St f = 0.76 is dominant while after that the response at St f = 0.54 is dominant.The bottom row shows qualitatively similar behaviour, where the wake flow response at St f = 0.54 is dominant initially (up to x/D ≈ 5) while the Low-order modelling of wake meandering 551 response at St f = 0.32 is dominant later.This shows that as the wake evolves in the downstream direction, the frequency content shifts toward the lower values as also observed in Foti et al. (2018b).The physical origin of such behaviour is the mean flow evolution, which spreads in the radial direction with the downstream distance.Thus, the frequency of the locally most unstable mode decreases with x/D (as shown in figure 6b).Because the model coefficients are directly based on the local stability results, the model, apart from minor quantitative differences, mimics this wake flow behaviour correctly.
The second point to note is the change in the dominant frequency with the forcing amplitude.The figure shows that the wake flow responses at all three St f increase from the top to the bottom row as A f increases.This increment in the response is maximum at St f = 0.32 and is minimum at St f = 0.76.Consequently, in the top row, the model response is maximum at St f = 0.76 (for x/D < 6) and at St f = 0.54 (for x/D > 6).While, in the bottom row, the model response is maximum at St f = 0.54 (for x/D < 5) and at St f = 0.32 (for x/D > 5).This shows that as the disturbance amplitude increases, the wake flow response shifts towards lower frequencies as also observed in Mao & Sørensen (2018).Again the model captures this behaviour correctly, thus, justifying our choice of the nonlinear coefficient (see § 4.2).The figure also shows that at the higher A f (bottom row), the wake flow response fluctuates around its curve fit.These fluctuations are because of the presence of the higher harmonics and their interaction with the fundamental mode.The higher harmonic was also observed in figure 5(a), but is not present in the CGL equation based model results (such as in figure 11a).This is discussed in appendix B. 5.1.2.Variations in the wake meandering features with the turbine operating conditions Figure 13 presents the same results as in figure 12 but for the wake flows corresponding to Ω = 12.0 r.p.m. and 9.0 r.p.m. cases.A comparison of the wake flow behaviour for different Ω cases reveal three important features.The wake meandering (i) occurs earlier in the streamwise direction, (ii) has higher amplitude and (iii) has nearly constant peak Strouhal number range as Ω increases.The first two features are due to the fact that thrust coefficient (C t ) increases with Ω (C t = {0.59,0.74, 0.86} for Ω = {9.0r.p.m., 10.5 r.p.m., 12.0 r.p.m.}, respectively).While the third feature is because wake meandering frequency scales with the incoming velocity and turbine diameter, the operating conditions play only a marginal role.Foti et al. (2018a) performed LES on a turbine operating at optimal and suboptimal (lower thrust) conditions and simulated them with and without a nacelle.They observed that wake meandering is delayed when turbine is operating at suboptimal condition as well as when simulated without a nacelle (similar observations are also reported in Kang, Yang & Sotiropoulos (2014)).They attributed this observation to the turbine drag force that determines the wake deficit and affects the mean velocity evolution.The local stability results in figure 6(b) also support this reasoning, where transition to the far-wake region occurs earliest in the Ω = 12.0 r.p.m. case and latest in the Ω = 9.0 r.p.m. case.Consequently, this observation is also reflected in the model, where (like the wake flow) the peak amplitude in each A f and St f case is achieved earliest in the Ω = 12.0 r.p.m. case and latest in the Ω = 9.0 r.p.m. case.
Higher drag force (or wake deficit) is also responsible for higher wake meandering amplitude.This is first shown in Yang et al. (2015), where the added turbulent kinetic energy in the far-wake region is shown to increase with the thrust coefficient, and is later also confirmed in Heisel et al. (2018)   in local stability results here, the spatial growth rate (−k si ) also increases with Ω, particularly in the region up to x/D = 6, and hence is also observed in the model results.
The results in figures 2, 12 and 13 show that the peak frequency does not vary much with the turbine's operating conditions.The small increase in the peak frequency with increasing Ω in figure 2 was attributed to the effect of the near-wake structures.In figures 12 and 13, there is a small reduction in the peak frequency with increasing Ω.This reduction is attributed to the earlier onset of wake meandering in higher Ω cases as is also reflected in the stability results in figure 6(b).These results thus corroborate that wake meandering frequency scales with the turbine diameter and the mean velocity.

System nonlinearity and transfer of energy from low to high frequencies
As mentioned at the end of § 4.1, the conversion from the linear to nonlinear CGL equation is not trivial.There are various options that exist in the literature.These include using (i) the nonlinear dispersion relations to obtain the CGL coefficients (Pier et al. 1998), (ii) a weakly nonlinear analysis to obtain the nonlinear term from the Navier-Stokes equations (Sipp & Lebedev 2007), (iii) an input-output analysis to find the nonlinear terms as well as their order (Lee et al. 2019) and (iv) a self-consistent framework to obtain the variation in the nonlinear gain with the forcing amplitude (Mantič-Lugo & Gallaire 2016).The first two methods are limited to self-sustained oscillations, while the latter two methods are computationally intensive as well as mathematically involved and hence are out of the scope of the present study.Here, we determine the nonlinear term to be proportional to ω nonlinear coefficient, as it correctly captures the change in wake meandering behaviour with the forcing amplitude.To further justify our choice of the nonlinear coefficient, we refer to the results in Keppens et al. (1999) where they showed the saturation amplitude to be higher at the lower wavenumbers (than the linearly most unstable wavenumber) in Kelvin-Helmholtz-type instabilities.
Physically, the origin of nonlinear saturation is the distortion of the base flow, which results in a mean flow profile that, via Reynolds stress, saturates the instability modes at finite amplitudes.This is first formulated by Stuart (1960) and calculated recently by Mantič-Lugo, Arratia & Gallaire (2014) and Meliga (2017) for oscillator flows and, more relevant to the present case, by Mantič-Lugo & Gallaire (2016) for the response to harmonic forcing in an amplifier flow.Although, we do not obtain the nonlinear term from rigorous calculations.Nevertheless, based on the present results, we speculate that in the far-wake region higher-frequency flow structures distort the wake flow faster and create more Reynolds stress as compared to lower-frequency flow structures.
Another role of nonlinearity is to transfer energy between the different frequencies that are otherwise linearly decoupled.Heisel et al. (2018) observed that very low-frequency disturbances transfer energy to relatively higher frequencies.We perturb the wake flow and the model by low-frequency disturbances (at St f = 0.11), and the spectra of their responses are shown in figure 14.Panel (a) shows that in the wake flow, energy is transferred to the higher frequencies and several peaks at the higher harmonics appear in the far-wake region.Panel (b) shows that in the model, energy from St f = 0.11 is also transferred to higher frequencies but there are no higher harmonic peaks in the far-wake region.This is because the CGL equation is inherently unable to capture the higher harmonics.An alternative to the CGL equation, in the form of a wave equation, that can capture the higher harmonics is formulated and discussed in appendix B.

Limitations of the low-order model
The main limitation of the present model that hinders quantitatively accurate predictions arises from the limitations in our understanding of the near-wake region.In contrast to the far-wake region, the flow instabilities and their interactions in the near-wake region are very sensitive to the turbine design and support structures (Kumar & Mahesh 2017).Reliable modelling and understanding the origin of flow instabilities in this region, therefore, may require turbine-geometry resolved high-fidelity simulations.This is beyond the scope of this study and is the reason that V. Gupta and M. Wan we avoided speculation on the origin of the near-wake structures and mode switchings observed in § § 2 and 3, respectively.The model coefficients in the near-wake region, therefore, are based only on extrapolation of the far-wake region (see § 4.2) and a rough match with the near-wake region (see § 5.1).Other limitations in the model arise from the strong nonlinearity (as discussed in § 5.2 and appendix B) and from lack of coupling of different azimuthal wavenumber modes.
Apart from that, there are limitations that arise from the simplifications in the simulations as compared to realistic conditions in utility-scale turbines.These include spatially uniform inflow conditions in place of a turbulent boundary layer mean flow profile or mean incoming flow at a yaw angle, which disrupt the axisymmetric nature of the flow.Other complexities, such as the spatial shape of turbulent structures, presence of support structures and further variations in the operating conditions increase the number of parameters to be considered in the model.To account for them, the present model may need data from input-output models developed in the literature (Hamilton et al. 2018).

Conclusion
Several recent studies have shown that wake meandering observed behind turbines have broad spectral peaks centred around St = 0.3 (Heisel et al. 2018;Foti et al. 2018b).Based on these observations and stability analysis of the wake flow profile (Mao & Sørensen 2018), the mechanism behind wake meandering is concluded to be the amplification of upstream structures via shear flow instabilities in the far-wake region.In this paper, we obtain unsteady wake flow fields behind a tidal turbine in a uniform rectangular channel using LES (in § 2).The inflow in our simulations is maintained as spatially uniform and is either steady or sinusoidally varying in time.The main features of the wake flow, particularly of wake meandering in the far-wake region, are shown to agree well with the literature.We also discuss how the near-wake structures can lead or contribute to meandering in the far-wake region via convective instabilities (figure 1 The main contribution of this paper is the development of a low-order dynamical system to phenomenologically model wake meandering (in § 4).The model is based on a modified form of the CGL equation for flows that are in the amplifier regime.The model reproduces the main qualitative features of wake meandering, such as (i) its origin via amplification of upstream structures, (ii) dependence of oscillation frequency on the upstream disturbance amplitude (figure 12), (iii) shift towards lower frequencies as the wake flow evolves in the streamwise direction (figure 12), and, to an extent, (iv) transfer of energy from very low frequencies towards relatively higher frequencies (figure 14).Additionally, the model also predicts the increase in meandering amplitude and the advancement of the onset of meandering closer to the turbine with the increasing thrust coefficient (figure 13).The main reason that the model is able to predict these features is because its coefficients are based on the results from local linear stability analysis performed over the time-averaged mean flow (in § 3).
To our knowledge, this is the first low-order dynamical system that models wake meandering and account for the mechanism of its origin.The agreement between the model and the simulations is very encouraging, however, there are still some limitations in the model that need to be resolved for quantitatively accurate predictions.The main limitations include determination of (i) the model coefficients in the near-wake region, and (ii) the nonlinear coefficient that can account for higher harmonics and coupling of different azimuthal modes.Other limitations arise from the simplistic and limited operating conditions explored in this paper.Given that these limitations can be satisfactorily resolved, the small number of degrees of freedom makes this model an ideal candidate for application in energy farm models.It can thus help in designing farm layouts to minimise the adverse effects of wake meandering.
Here, we also propose a wave equation (given below) that is based on coupled oscillators.We show its equivalence with the CGL equation in linear and weakly nonlinear regimes, as well as its ability to capture the higher harmonics in strongly nonlinear regime Here, o, ε, b, µ, C di , u g and u gi are space-dependent real variables.We follow Nayfeh (1982) to derive an amplitude equation (i.e. the CGL equation) from the wave equation (see Lee et al. (2019) for more details).Towards that purpose, we first apply a variation of parameter and transform u(x, t) into variables |A(x, t)| and θ (t) as u = |A| cos(ω f t + θ), u = −|A|ω f sin(ω f t + θ ), (B 2a,b) where ω f is the global frequency (for oscillator flows) or the forcing frequency (for amplifier flows).Conditions in (B 2) also give the following two equations as 0 = |A| cos(ω f t + θ) − θ|A| sin(ω f t + θ ), (B 3) We insert (B 2) to (B 4) in (B 1) and apply the trigonometric identities to obtain two first-order differential equations in time as It should be noted, that until this point we have made no assumption about the system behaviour.We now make an assumption that the system is weakly nonlinear.Hence, the higher harmonics are weak and |A| and θ act like system's amplitude and phase, and vary sinusoidally with t at the fundamental frequency.The application of the method of averaging then gives the CGL equation governing the evolution of A = |A| exp(iω f t + iθ) as In order to make (B 1) equivalent to the model (4.6), we adjust its coefficients based on the form in (B 7) and the local stability results as o = 2ω f ω sr − ω 2 f , ε = −2(∂ω/∂k)(ω s , k s )k si , µ = −(∂ 2 ω/∂k 2 )| i , C di = 0, u g = 2(∂ω/∂k)| r and u gi = 0.The nonlinear term is set as εb = 32ω 1.5  f .The results of the wave equation (B 1) are compared with the CGL equation (4.6) in figure 16 where (a) linear gains (G L ) and (b-d) nonlinear responses (in terms of the spectra at {X = 4, 5, 6, 7, 8, 9} for A f = 0.135) at St f = 0.11, 0.54 and 0.97 are plotted.Panel (a) confirms the equivalence of the CGL and the wave equation in the linear limit.Even under the strong forcing, panels (b-d) show that the nonlinear responses match between the two models, except in the later regions in the St f = 0.11 case.In this case, the wave equation results show the emergence of a higher harmonic that is missing in the CGL equation results.
FIGURE 1. (Colour online) (a) The instantaneous axial velocity field in the x-z plane passing through y = 0 in Ω = 10.5 case.The PSDs of w-velocity fluctuations in the (b) near-wake region (the circles in panel a) show the tip and root vortices ( f t = 0.350 Hz) and other near-wake vortices ( f n ≈ 0.219 Hz), (c) far-wake region (the squares in panel a) show the far-wake oscillations have a broad peak at f m ≈ 0.048 Hz and (d) along the x-direction (the crosses at x/D = 0.0, 2.5, 4.0, 6.0 and 9.0 in panel a) shows that the far-wake oscillations generate from amplification of the near-wake structures at 0.048 Hz (St ≈ 0.52).
FIGURE 2. (Colour online) (a,c,e,g) The instantaneous axial velocity fields in the x-z plane at y = 0 and corresponding (b,d,f,h) φ w at x/D = [2.5, 4.0, 6.0, 9.0] and z/D = 0.325.Wake meandering is present in all the cases.Except for the Ω = 9.0 r.p.m. case, there is a peak in φ w at St ≈ 0.5.

FIGURE 3 .
FIGURE 3. The instantaneous axial velocity fields under sinusoidally varying inflow conditions St f = (a-c) 0.11, (d-f ) 0.43, (g-i) 0.76 and ( j-l) 0.97.The forcing amplitudes (increasing from the bottom to top) are mentioned in the respective panels.The wake meandering pattern becomes regular (except for St f = 0.11) as A f increases.
FIGURE 4. (Colour online) The wake flow response to sinusoidal forcing in terms of r f = φ w (St f )/(φ w (St f ) + φ w [St m ]), which is nearly zero when the effect of the background noise is dominant (i.e.φ w [St m ] φ w (St f )) and approaches one as the effect of the sinusoidal forcing becomes dominant (i.e.φ w (St f ) φ w [St m ]).The black line represents the pseudo lock-in curve (r f ≈ 0.95), above which the sinusoidal forcing has suppressed the effect of the background noise.Pseudo lock-in is achieved earlier (i.e. at lower A f ) when St f is closer to St m (≈ 0.52).

FIGURE 5 .
FIGURE 5.The wake flow response to sinusoidal forcing at frequencies (a-c) 0.32, (d-f ) 0.54 and (g-i) 0.76 and at varying A f (increasing from the bottom to top) in terms of the velocity fluctuations frequency spectra at (x/D, z/D) = (9, 0.325).As A f increases, the maximum wake flow response shifts from at St f = 0.76 to at St f = 0.32.Thus, showing a shift to lower frequencies with the increasing forcing amplitudes.
FIGURE 7. (Colour online) The coefficients are determined using least-square polynomial curve fits (red lines) to the local stability results at m = −1 (black lines) in the region X = 4.5-10.0.
FIGURE 9. (Colour online) The model response to white Gaussian noise forcing (A w = 0.0004) in terms of the PSDs (φ A ) at X = {1, 2, 3, . . ., 10}.The oscillations evolve towards lower frequencies with increasing X (indicated by the arrow) and finally have a broad peak at St ≈ 0.5.

FIGURE 11 .
FIGURE 10. (Colour online) The model response to sinusoidal forcing in terms of r fm , which is nearly zero when the background noise is dominant and approaches one as the sinusoidal forcing becomes dominant.The black line represents the pseudo lock-in curve (r fm ≈ 0.95), above which the sinusoidal forcing has suppressed the effect of the background noise.Pseudo lock-in is achieved at lower A f when St f is close to St m (≈ 0.52). https://doi.org/10.1017/jfm.2019.
FIGURE 12. (Colour online) Comparison of the wake flow and the model responses to forcing at frequencies (a,b) 0.32, (c,d) 0.54 and (e,f ) 0.76 and at two forcing amplitudes.The wake meandering frequency shifts to lower values with increasing (i) x/D and (ii) A f . https://doi.org/10.1017/jfm.2019.
FIGURE 13. (Colour online) Same results as in figure 12 but for the Ω = 12.0 r.p.m. and 9.0 r.p.m. cases.A comparison between different Ω cases shows that wake meandering (i) occurs earlier in the streamwise direction and (ii) has higher amplitude as Ω increases (because the thrust coefficient increases with Ω). (iii) The Strouhal number range, however, remains nearly unaffected by changing Ω.All these features are well captured by the low-order model.
FIGURE 14. (Colour online) The (a) wake flow and (b) model response to forcing at a low frequency (St f = 0.11).They both show a transfer of energy from the low frequency (at x/D = 4) to relatively higher frequencies at St = 0.2-0.8(at x/D = 9).The model, however, is unable to account for the transfer of energy to the higher harmonics.

FIGURE 15
FIGURE 15. (Colour online) A local dispersion relation and its second-order approximation around (ω s , k s ) in terms of the variations of the (a) real and (b) imaginary parts of the wavenumber with the real frequency.Equation (4.6) is based on the second-order approximation.
FIGURE 16. (Colour online) The CGL versus wave model (a) linear gain and (b-d) nonlinear responses (as the spectra at {X = 4, 5, 6, 7, 8, 9} for A f = 0.135) at St f = 0.11, 0.54 and 0.97.The CGL and wave models show a good agreement except for the nonlinear response in St f = 0.11 case, where the higher harmonics appear in the wave model response.