<figure><img src="image_1.png" /></figure>

HAL
open science

# Harmonics generation and the mechanics of saturation in flow over an open cavity: a second-order self-consistent description

Philippe Meliga

► To cite this version:

Philippe Meliga. Harmonics generation and the mechanics of saturation in flow over an open cavity: a second-order self-consistent description. Journal of Fluid Mechanics, 2017, 826, pp.503 - 521. 10.1017/jfm.2017.439. hal-01585331

HAL Id: hal-01585331

https://hal.science/hal-01585331v1

Submitted on 19 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Harmonics generation and the mechanics of saturation in flow over an open cavity : a second-order, self-consistent description

P. Meliga¹

¹ Aix-Marseille Université, CNRS, Ecole Centrale Marseille, Laboratoire M2P2, Marseille, France

(Received ?; revised ?; accepted ?. - To be entered by editorial office)

The flow over an open cavity is an example of supercritical Hopf bifurcation leading to periodic limit cycle oscillations. One of its distinctive features is the **existence of strong** higher harmonics, which results in the time-averaged mean flow being strongly linearly unstable. For this class of flows, a simplified formalism capable of unraveling how exactly the instability grows and saturates is lacking. **This study** builds on previous work by Mantič-Lugo et al. (Phys. Rev. Lett. 113:084501, 2014) to fill in the gap using a parameterized approximation of the instantaneous mean flow, coupled in a quasi-static manner to multiple linear harmonic disturbances interacting nonlinearly with one another and feeding back on the mean flow via their Reynolds stresses. This provides a self-consistent modeling of the mean flow/fluctuation interaction, in the sense that all perturbation structures are those whose Reynolds stresses force the mean flow in a way such that the mean flow generates exactly the aforementioned perturbations. The first harmonic is sought as the superposition of two components, a **linear** component generated by the instability and aligned along the leading eigenmode of the mean flow, and a **nonlinear** orthogonal component generated by the higher harmonics, that progressively distorts
the linear growth rate and eigenfrequency of the eigenmode. Saturation occurs when the growth rate of the first harmonic is zero, at which point the stabilizing effect of the second harmonic balances exactly the linear instability of the eigenmode. The model does not require any input from numerical or experimental data, and accurately predicts the transient development and the saturation of the instability, as established from comparison to time and ensemble-averages of direct numerical simulation data.

# 1. Introduction

For flows whose unsteadiness proceeds from an intrinsic instability mechanism, linear and weakly nonlinear stability analysis provide a rigorous mathematical foundation to explain the initial growth of disturbances of a certain scale and frequency, and the leading-order nonlinear effects associated to small-amplitude disturbances. The quality of the prediction however deteriorates rapidly with increasing distance from the instability threshold, as the perturbative nature of these approaches imposes to build the fluctuation as successive-order corrections to the leading eigenmode of the base flow (i.e., the steady state sustaining the instability), whose spatial structure can differ considerably from that of the nonlinearly saturated oscillation (Dušek, Le Gal & Fraunié 1994; Noack, Afanasiev, Morzynski, Tadmor & Thiele 2003). As an illustration, the vortex street of shed vortices that form the wake of a circular cylinder sets in at Reynolds number $Re = 47$, but the inaccuracy on the frequency computed by linear stability analysis is already by 30% at values of $Re$ as low as 80 (Barkley 2006), while that on the weakly nonlinear oscillation amplitude exceeds 100% (Mantič-Lugo, Arratia & Gallaire 2014).

It is the time-averaged mean flow, not the base flow, that is key in understanding the growth of disturbances while encompassing the nonlinearity of the system (Maurel,
Pagneux & Wesfreid 1995; Zielinska, Goujon-Durand, Dušek & Wesfreid 1997; Noack
et al. 2003). A simplified model for the saturation of the cylinder flow has been proposed
recently by Mantič-Lugo et al. (2014) using a so-called self-consistent approach coupling
the mean flow and its leading eigenmode via formation of Reynolds stresses. The authors
show that saturation occurs when the eigenmode has zero linear growth rate, wherein
the model predicts accurately the mean flow, Reynolds stresses, frequency, amplitude and
structure of the fluctuation obtained by direct numerical simulation (DNS). The model
is also relevant to describe the transient development of the instability (Mantič-Lugo,
Arratia & Gallaire 2015), which formalizes the intuitive picture invoked to describe the
growth and saturation of disturbances, dating back to Malkus (1956) and Stuart (1971):
perturbations feed on the unstable base flow and grow first according to the linear theory.
When they reach a size large enough for nonlinearities to set in, they feed back via their
Reynolds stresses, and distort the base flow into an increasingly stable mean flow, up
to the point where the mean flow is linearly marginally stable and perturbations stop
growing and saturate. It also rationalizes the fact that the mean flow eigenfrequency
agrees well with the exact oscillation frequency, as early noticed in cylinder flows and
related bluff-body wakes (Hammond & Redekopp 1997; Pier 2002; Barkley 2006).

The above scenario anticipates a linearly marginally stable mean flow because it re-
duces the consequences of the nonlinearity to mean flow distortions, and neglects another
important mechanism that may alter the predictions of the linear theory: the generation
of harmonics. The model of Mantič-Lugo et al. (2014) relies on the same simplifying as-
sumption. It is thus relevant to the cylinder case, whose mean flow is linearly marginally
stable (Barkley 2006) only because vortex-shedding is almost monochromatic (Dušek
et al. 1994). For flows whose nonlinearity unveils strong higher harmonics, e.g., open cav-
ity flows (Sipp & Lebedev 2007), turbulent wakes (Meliga, Sipp & Chomaz 2009b; Meliga,
Pujals & Serre 2012; Meliga, Cadot & Serre 2016b), or certain regimes of thermosolu-
tal convection flows (Turton, Tuckerman & Barkley 2015), the mean flow is generally
strongly linearly unstable and a simplified formulation capable of similarly accurately
unraveling the mechanics of saturation is lacking. This research intends to fill in the gap
using augmented self-consistent models coupling the mean flow to multiple harmonic lin-
ear disturbances interacting nonlinearly with one another and feeding back on the mean
flow via formation of Reynolds stresses. Such an approach is applied to the flow over a
square cavity, for which we show that only the second harmonic needs to be retained.

The paper is organized as follows : the configuration is described in § 2, and the
unsteady cavity flow is characterized in § 3 using DNS and classical linear stability anal-
ysis. The second-order self-consistent model is derived in § 4, and its ability to describe
the mechanics of saturation is assessed in § 5 from exhaustive comparison to time and
ensemble-averages of DNS data. Connection with weakly nonlinear analysis is addressed
in § 6, where we also evidence the ability of the method to encompass the effect of
higher-order harmonics, and discuss how to anticipate an appropriate truncation order.

## 2. Flow configuration

We consider the two-dimensional, incompressible flow over the open square cavity
described in details by Sipp & Lebedev (2007). Suffice it to say here that a uniform
velocity field is imposed at the inlet boundary and that a laminar boundary layer starts
to develop between the inlet and the upstream edge of the cavity, the exact position
being adjustable via an appropriate choice of boundary conditions; see figure 1(a). The
free-stream velocity and cavity length (or height) are used to make all quantities non-
dimensional. We denote by **u** = (u, v) the velocity field of components u and v in the
A second-order self-consistent model of saturation in the open cavity flow 5
stream-wise x and cross-stream y directions and by p the pressure. The flow is governed
by the Navier-Stokes equations (NSE) written in compact form as

$$
\partial_t \mathbf{u} + \mathbf{N}(\mathbf{u}) = \mathbf{0}, \qquad (2.1)
$$

where

$$
\mathbf{N}(\mathbf{u}) = \mathbf{u} \cdot \nabla \mathbf{u} + \nabla p - Re^{-1} \nabla^2 \mathbf{u}, \quad (2.2)
$$

is the Navier-Stokes operator, whose dependency on p is omitted to ease the notation.
Because of incompressibility, it is understood that all velocity fields are divergence free,
therefore we will not write this condition explicitly.

All calculations are performed on the same computational domain made of 194,447
triangles, equipped with the L² inner product (·|·) and related norm ||·||. We perform
direct numerical simulations (DNS) of equations (2.1) with the finite-elements solver pre-
sented in Meliga, Boujo, Pujals & Gallaire (2014). Time integration is achieved using a
second-order accurate Crank-Nicholson scheme with time step Δt = 10⁻². Unless speci-
fied otherwise, all time-averaged quantities (as denoted by ⃗) are obtained by averaging
on-the-fly over 200 time units after the flow has settled down to a periodic regime, which
represents ~ 300 cycles depending on the Reynolds number. We also use the Newton
method to compute the base cavity flow u_b, i.e., the solution to the steady NSE

$$
\mathbf{N}(\mathbf{u}_b) = \mathbf{0}, \tag{2.3}
$$

whose stability is analyzed using the Arnoldi method. In the following, we denote by $\hat{\mathbf{u}}_{10b}$
the leading eigenmode of linear growth rate $\sigma_{0b}$ and eigenfrequency $\omega_{0b}$, hence solution
to the eigenvalue problem

$$
(\sigma_{0b} + i\omega_{0b})\hat{\mathbf{u}}_{10b} + \mathbf{L}(\mathbf{u}_b)\hat{\mathbf{u}}_{10b} = \mathbf{0}, \quad (2.4)
$$
<figure><img src="image_3.png" /><figcaption>FIGURE 1. (a) Iso-contours of vorticity showing the periodic ejection of vortices from the upstream edge of the cavity - Re = 6000. The stream-wise velocity spectrum in figure 3 is measured at the position (x, y) = (0.7, 0) of the grey symbol, and the boundary layer develops from the position (x, y) = (-0.4, 0) of the yellow symbol. (b) Oscillation frequency $\omega_{nl}$ against Re.</figcaption></figure>

where

$$L(\mathbf{u}_b)\hat{\mathbf{u}} = \mathbf{u}_b \cdot \nabla \hat{\mathbf{u}} + \hat{\mathbf{u}} \cdot \nabla \mathbf{u}_b + \nabla \hat{p} - Re^{-1}\nabla^2 \hat{\mathbf{u}}, \quad (2.5)$$

is the linearized Navier-Stokes operator (LNS). We also determine the linear stability of the time-averaged mean flow $\bar{\mathbf{u}}$, whose leading eigenmode $\hat{\mathbf{u}}_{10m}$ of linear growth rate $\sigma_{0m}$ and eigenfrequency $\omega_{0m}$ is solution to the eigenvalue problem

$$(\sigma_{0m} + i\omega_{0m})\hat{\mathbf{u}}_{10m} + L(\bar{\mathbf{u}})\hat{\mathbf{u}}_{10m} = \mathbf{0}. \quad (2.6)$$

### 3. Direct numerical simulation and eigenmode analysis

Figure 1(a) shows a snapshot of vorticity at Re = 6000. The shear layer rolls-up into large-scale vortices that ultimately impact the downstream edge of the cavity. A pressure wave forms, travels upstream via the recirculating flow in the cavity, and excites the shear layer at the upstream edge. This causes new perturbations to grow again into large-scale vortices via Kelvin-Helmholtz instability. This feedback loop leads to a linear global instability, and the flow settles on a periodic limit cycle (Åkervik, Ehrenstein, Gallaire & Henningson 2008; Barbagallo, Sipp & Schmid 2009). The curve of the oscillation frequency against Re shown in figure 1(b) presents a discontinuity at Re = 4400,
<figure><img src="image_4.png" /><figcaption>FIGURE 2. (a) Linear growth rate against $Re$. The value $\sigma_{0m}$ (resp. $\sigma_{0b}$) obtained by linear stability analysis of the time-averaged mean flow (resp. of the base flow) is shown by the grey line (resp. the dash-dotted grey line). (b) Same as (a) for the eigenfrequency $\omega_0$. The nonlinear frequency $\omega_{nl}$ extracted from the DNS is shown by the black line.</figcaption></figure>

which is because the base flow undergoes successive Hopf bifurcations at $Re_{c1} = 4114$ (in good agreement with the threshold value 4140 reported by Sipp & Lebedev 2007) and $Re_{c2} = 4348$. Both eigenmodes exchange dominance at $Re = 4567$, after which the mode selection follows the ‘largest growth rate’ criterion. This bifurcation sequence is robust in the sense that, while small changes in the set-up (especially the position at which the boundary layer starts to develop) yield small variations in the instability thresholds, the flow oscillations systematically end up being driven by the second eigenmode at sufficiently high Reynolds numbers (not shown here). The focus in the sequel is thus on the range $Re \ge 4400$, with default value set to $Re = 6000$.

The mean cavity flow fails to comply with the linear marginal stability criterion of Mantič-Lugo et al. (2014), as evidenced by the large, positive growth rate $\sigma_{0m} = 0.117$ of the leading eigenmode at $Re = 6000$. Figure 2(a) unveils a similar behavior regardless of the Reynolds number, which is because the nonlinearity generates strong higher harmonics, as evidenced in figure 3 by the spectrum of the stream-wise velocity measured in the shear layer at $(x, y) = (0.7, 0)$, almost midway between both edges of the cavity: beyond the peak at the fundamental frequency $\omega_{nl} = 11.19$, we notice a substantial
<figure><img src="image_5.png" /><figcaption>FIGURE 3. Spectrum of the stream-wise shear-layer velocity at $(x, y) = (0.7, 0)$ - $Re = 6000$.</figcaption></figure>

All results normalized for the fundamental to peak at unity.

second harmonic at $2\omega_{nl}$, and even third and fourth harmonics (although the peaks at $3\omega_{nl}$ and $4\omega_{nl}$ are of much smaller magnitudes). Interestingly, the eigenfrequency of the mean flow $\omega_{0m} = 11.22$ does predict well the nonlinear oscillation frequency (unlike that $\omega_{0b} = 10.69$ of the base flow; see also figure 2(b)), which is a little surprising because, according to Sipp & Lebedev (2007), the mean flow eigenfrequency is expected to miss on the exact frequency in the presence of strong higher harmonics. Such an interpretation however issues from an asymptotic analysis valid only close to the instability threshold, while linearly unstable mean flows whose eigenfrequency somehow manages to do well on the frequency prediction have been documented in the context of strongly nonlinear turbulent wakes (Meliga et al. 2009b, 2012, 2016b). It is part of the objectives of this research to show that the nonlinearity induced by the higher-order harmonics can help straighten out this apparent contradiction.

# 4. Self-consistent model

## 4.1. *Theoretical formalism*

The self-consistent theory is meant to describe the development and saturation of the instability based on a parameterized approximation of the mean flow/fluctuation interaction. The first step is to apply Reynolds decomposition and to split the instantaneous
A second-order self-consistent model of saturation in the open cavity flow 9

solution into $\mathbf{u} = \langle \mathbf{u} \rangle + \mathbf{u}'$, where $\langle \cdot \rangle$ denotes the ensemble average, $\langle \mathbf{u} \rangle$ is the instantaneous mean flow, and $\mathbf{u}'$ is the zero-mean fluctuation. Substituting in the NSE yields

$$ \partial_t \langle \mathbf{u} \rangle + \mathbf{N}(\langle \mathbf{u} \rangle) = - \langle \zeta(\mathbf{u}') \rangle, \quad (4.1a) $$

$$ \partial_t \mathbf{u}' + \mathbf{L}(\langle \mathbf{u} \rangle)\mathbf{u}' = -\zeta(\mathbf{u}') + \langle \zeta(\mathbf{u}') \rangle, \quad (4.1b) $$

where we denote by $\zeta(\mathbf{u}) = \mathbf{u} \cdot \nabla \mathbf{u}$ the nonlinear advection operator. The second step is to approximate the fluctuation by a Fourier series expansion. We elaborate further on this matter in § 6.2, but we settle here for a second-order expansion, as we anticipate from figure 3 that the effect of the higher-order harmonics will be of a much smaller magnitude. This yields

$$ \mathbf{u}'(t) = A(t)\mathbf{u}_1(t)e^{i\omega t} + A^2(t)\mathbf{u}_2(t)e^{2i\omega t} + \text{c.c.}, \quad (4.2) $$

where $\omega$ is the fundamental oscillation frequency, $\mathbf{u}_1$ is the structure of the first harmonic, parameterized by its (real) amplitude $A$, $\mathbf{u}_2$ is the structure of the second harmonic, coming (with no loss of generality) with amplitude $A^2$, and we allow both the amplitude and the structure of the harmonics to vary in time to reflect the distortions of the fluctuation throughout the development of the instability.

As explained in Mantič-Lugo et al. (2015), the premise of the self-consistent theory is that there exists a separation of scales between the fast time scale on which the flow oscillates, and the slow time scale on which the oscillations grow in amplitude (similar to that assumption made in weakly nonlinear analysis to compute the amplitude of the unstable eigenmode from an amplitude equation). For the cavity flow, the relevance of this assumption is supported by the results documented in § 5.2, therefore we invoke a quasi-static approximation to slave the mean flow $\langle \mathbf{u} \rangle$ and the harmonics $\mathbf{u}_1$ and $\mathbf{u}_2$ to their forcing by the self-interaction of the fluctuation, as described by the right-hand sides (RHSs) of system (4.1). In practice, this amounts to assuming at each time instant a
steady mean flow $\langle \mathbf{u} \rangle$ and purely oscillating harmonics $\mathbf{u}_1$ and $\mathbf{u}_2$, for which the ensemble average is equivalent to a temporal average over a period. Substituting in (4.1) and neglecting all higher interactions yields

$$ \mathbf{N}(\langle \mathbf{u} \rangle) = -A^2 \psi(\mathbf{u}_1) - A^4 \psi(\mathbf{u}_2), \quad (4.3a) $$

$$ (\sigma + i\omega)\mathbf{u}_1 + \mathbf{L}(\langle \mathbf{u} \rangle)\mathbf{u}_1 = -A^2\phi(\mathbf{u}_1^*, \mathbf{u}_2), \quad (4.3b) $$

$$ 2(\sigma + i\omega)\mathbf{u}_2 + \mathbf{L}(\langle \mathbf{u} \rangle)\mathbf{u}_2 = -\zeta(\mathbf{u}_1), \quad (4.3c) $$

where $\phi(\mathbf{u}, \mathbf{v}) = \mathbf{u} \cdot \nabla \mathbf{v} + \mathbf{v} \cdot \nabla \mathbf{u}$ is the linearized advection operator, * indicates complex conjugation, and we note $\psi(\mathbf{u}) = \phi(\mathbf{u}, \mathbf{u}^*)$. In equation (4.3b), $\sigma$ is the growth rate of the first harmonic (assumed to be relative neither to a particular point in space, nor to a particular component of the solution), and we take the second harmonic to grow at the same rate $2\sigma$ as the source term $\zeta(\mathbf{u}_1)$ in (4.3c), hence

$$ \partial_t (A \mathbf{u}_1) = \sigma A \mathbf{u}_1, \quad \text{and} \quad \partial_t (A^2 \mathbf{u}_2) = 2\sigma A^2 \mathbf{u}_2. \quad (4.4) $$

Equation (4.3a) now defines the mean flow as a solution to the steady NSE forced by the Reynolds stresses of the fluctuation, while equations (4.3b)-(4.3c) conversely define the harmonics as forced solutions to the NSE linearized about the mean flow, with the RHS of equation (4.3b) (resp. the RHS of equation (4.3c)) accounting for the self-interaction of the fluctuation at frequency $\omega$, that feeds back on the first harmonic (resp. the self-interaction at frequency $2\omega$, that generates the second harmonic).

System (4.3) does provide a self-consistent description of the mean flow/fluctuation interaction, in the sense that all perturbation structures are those whose Reynolds stresses force the mean flow in a manner such that the mean flow generates exactly the aforementioned perturbations. However, it cannot be solved as is because both $\sigma$ and $\omega$ are part of the unknowns. Therefore, we further decompose the first harmonic into

$$ \mathbf{u}_1 = \hat{\mathbf{u}}_{10} + \mathbf{u}_{1\perp}, \quad (4.5) $$
A second-order self-consistent model of saturation in the open cavity flow 11
where $\hat{\mathbf{u}}_{10}$ is the leading eigenmode of the mean flow, of linear growth rate $\sigma_0$ and
eigenfrequency $\omega_0$, hence solution to the eigenvalue problem

$$ (\sigma_0 + i\omega_0)\hat{\mathbf{u}}_{10} + \mathbf{L}(\langle\mathbf{u}\rangle)\hat{\mathbf{u}}_{10} = \mathbf{0}, \quad (4.6) $$

and $\mathbf{u}_{1\perp}$ is the projection on the hyperplane spanned by the remaining eigenmodes. It
is worth insisting that equation (4.6) reflects by no means the time-dependency of the
eigenmode. The model takes both $\hat{\mathbf{u}}_{10}$ and $\mathbf{u}_{1\perp}$ to have the same growth rate $\sigma$ and
frequency $\omega$ as the first harmonic, as we shall clarify that those precisely differ from the
values $\sigma_0$ and $\omega_0$ obtained by linear stability analysis on behalf of the resonance between
the first harmonic and its nonlinear interaction with the second harmonic. In the follow-
ing, we term $\mathbf{u}_{1\perp}$ the orthogonal component of the first harmonic, which is a small misuse
of language because, owing to the non-normality of the LNS operator (Chomaz 2005),
$\mathbf{u}_{1\perp}$ is orthogonal not to the leading eigenmode, but to the *adjoint* leading eigenmode of
the mean flow, i.e., the solution $\hat{\mathbf{u}}_{10}^\dagger$ to the eigenvalue problem

$$ (\sigma_0 - i\omega_0)\hat{\mathbf{u}}_{10}^\dagger + \mathbf{L}^\dagger(\langle\mathbf{u}\rangle)\hat{\mathbf{u}}_{10}^\dagger = \mathbf{0}, \quad (4.7) $$

where

$$ \mathbf{L}^\dagger(\langle\mathbf{u}\rangle)\hat{\mathbf{u}} = -\langle\mathbf{u}\rangle \cdot \nabla\hat{\mathbf{u}} + \hat{\mathbf{u}} \cdot \nabla\langle\mathbf{u}\rangle^T + \nabla\hat{p} - Re^{-1}\nabla^2\hat{\mathbf{u}}, \quad (4.8) $$

is the adjoint of the LNS operator. The governing equation for $\mathbf{u}_{1\perp}$ is obtained subtracting (4.6) from (4.3b), hence

$$ (\sigma_0 + i\omega_0)\mathbf{u}_{1\perp} + \mathbf{L}(\langle\mathbf{u}\rangle)\mathbf{u}_{1\perp} = -A^2\phi(\mathbf{u}_1^*, \mathbf{u}_2) - \gamma\mathbf{u}_1, \quad (4.9) $$

where the detuning parameter

$$ \gamma = \gamma_r + i\gamma_i = (\sigma - \sigma_0) + i(\omega - \omega_0) = -A^2 \frac{(\hat{\mathbf{u}}_{10}^\dagger | \phi(\mathbf{u}_1^*, \mathbf{u}_2))}{(\hat{\mathbf{u}}_{10}^\dagger | \hat{\mathbf{u}}_{10})}, \quad (4.10) $$

is made explicit taking the inner product of (4.9) with $\hat{\mathbf{u}}_{10}^\dagger$ and integrating by parts.
We can now recast (4.3) into the alternative, closed system

$$ \mathbf{N}(\langle \mathbf{u} \rangle) = -A^2 \psi(\hat{\mathbf{u}}_{10} + \mathbf{u}_{1\perp}) - A^4 \psi(\mathbf{u}_2), \quad (4.11a) $$

$$ (\sigma_0 + i\omega_0)\hat{\mathbf{u}}_{10} + \mathbf{L}(\langle\mathbf{u}\rangle)\hat{\mathbf{u}}_{10} = \mathbf{0}, \quad (4.11b) $$

$$ (\sigma_0 - i\omega_0)\hat{\mathbf{u}}_{10}^{\dagger} + \mathbf{L}^{\dagger}(\langle\mathbf{u}\rangle)\hat{\mathbf{u}}_{10}^{\dagger} = \mathbf{0}, \quad (4.11c) $$

$$ (\sigma_0 + i\omega_0)\mathbf{u}_{1\perp} + \mathbf{L}(\langle\mathbf{u}\rangle)\mathbf{u}_{1\perp} = -A^2\phi(\hat{\mathbf{u}}_{10}^* + \mathbf{u}_{1\perp}^*, \mathbf{u}_2) - \gamma(\hat{\mathbf{u}}_{10} + \mathbf{u}_{1\perp}), \quad (4.11d) $$

$$ 2(\sigma_0 + i\omega_0 + \gamma)\mathbf{u}_2 + \mathbf{L}(\langle\mathbf{u}\rangle)\mathbf{u}_2 = -\zeta(\hat{\mathbf{u}}_{10} + \mathbf{u}_{1\perp}), \quad (4.11e) $$

$$ \gamma = -A^2 \frac{(\hat{\mathbf{u}}_{10}^{\dagger} | \phi(\hat{\mathbf{u}}_{10}^{*} + \mathbf{u}_{1\perp}^{*}, \mathbf{u}_2))}{(\hat{\mathbf{u}}_{10}^{\dagger} | \hat{\mathbf{u}}_{10})}, \quad (4.11f) $$

for $\langle \mathbf{u} \rangle$, $\{\hat{\mathbf{u}}_{10}, \sigma_0, \omega_0\}$, $\mathbf{u}_{1\perp}$, $\mathbf{u}_2$, $\hat{\mathbf{u}}_{10}^{\dagger}$ and $\gamma$. Those are the equations of our second-order self-consistent model, whose sole free parameter for given Reynolds number is the amplitude $A$ of the first harmonic (of course, if $\mathbf{u}_2 = \mathbf{0}$, then $\gamma$ and $\mathbf{u}_{1\perp}$ are trivially zero and the model consistently reduces to the first-order model of Mantič-Lugo et al. 2014). As further developed in the following, solving the model equations for increasing values of the amplitude shapes an implicit relationship, e.g., $\sigma_0 = \sigma_0(A)$, $\omega_0 = \omega_0(A)$, $\gamma = \gamma(A)$, etc., that links the obtained family of solutions to the development in time of the instability.

## 4.2. Numerical resolution

System (4.11) can be solved using a simple nested-loop scheme: assuming relevant guesses are at hand for all quantities of interest, an outer loop solves iteratively $\langle \mathbf{u} \rangle$ to a precision of $10^{-12}$ (in $L^2$ norm) with the Newton method, then updates $\sigma_0, \omega_0, \hat{\mathbf{u}}_{10}$ and $\hat{\mathbf{u}}_{10}^{\dagger}$ with the Arnoldi method and normalizes the eigenmodes using $||\hat{\mathbf{u}}_{10}|| = 1$ and $(\hat{\mathbf{u}}_{10}^{\dagger} | \hat{\mathbf{u}}_{10}) = 1$. This is done using the first-order, self-consistent solver presented in Meliga, Boujo & Gallaire (2016a). An inner loop then solves iteratively $\mathbf{u}_{1\perp}$ and $\mathbf{u}_2$, also with the Newton method, only the inner convergence threshold starts from an $O(1)$ value and decreases down to $10^{-12}$ (also in $L^2$ norm) following a geometric-progression series. After the first five to ten outer passes, this was found to allow converging $\mathbf{u}_{1\perp}$ and $\mathbf{u}_2$ to the desired precision
<figure><img src="image_6.png" /><figcaption>FIGURE 4. Mean oscillation amplitude of the leading eigenmode $e_{10} = A||\hat{\mathbf{u}}_{10}||$, of the orthogonal component $e_{1\perp} = A||\mathbf{u}_{1\perp}||$, of the first harmonic $e_1 = A||\hat{\mathbf{u}}_{10} + \mathbf{u}_{1\perp}||$ and of the second harmonic $e_2 = A^2||\mathbf{u}_2||$. All amplitudes are measured by the root-mean-square of the corresponding fluctuating energy, and obtained solving the self-consistent system (4.11) for increasing values of the amplitude $A - Re = 6000$.</figcaption></figure>

within one single iteration of the inner Newton and thereby to speed-up considerably the resolution process. At each inner iteration, $\gamma$ is updated using the compatibility condition (4.11f), which is mandatory to ensure that equation (4.11d) admits a solution, as $(\sigma_0 + i\omega_0)\mathbf{I} + \mathbf{L}$ is otherwise non-invertible. The orthogonality between $\mathbf{u}_{1\perp}$ with $\hat{\mathbf{u}}_{10}^\dagger$ is ultimately enforced retaining the contribution $(\hat{\mathbf{u}}_{10}^\dagger | \mathbf{u}_{1\perp}) \hat{\mathbf{u}}_{10} - \mathbf{u}_{1\perp}$ of the obtained numerical solution ($\mathbf{u}_{1\perp}$ can be obtained in one go by moving the $\gamma \mathbf{u}_{1\perp}$ term in the LHS of equation (4.11d), as $(\sigma + i\omega)\mathbf{I} + \mathbf{L}$ is generally invertible, but the computational cost then soars because a different operator must be inverted at each inner iteration). Provided the corrections made at each outer iteration are under-relaxed, the present scheme achieves robust convergence within 40 outer iterations at $Re = 4500$ (using relaxation factors of $\sim 0.9$), up to 400 iterations at $Re = 7000$ (using relaxation factors of $\sim 0.15$). This dramatic increase in the numerical cost is due to the fact that strong nonlinearities considerably slow down the convergence rate of iterative methods. Mantič-Lugo et al. (2015) also report difficulties in solving their first-order model for values of $Re$ beyond 120, while Fornberg (1980) discusses similar issues computing steady solutions to the NSE at increasingly large Reynolds numbers.
<figure><img src="image_7.png" /><figcaption>FIGURE 5. (a) Linear growth rate $\sigma_0$ of the leading eigenmode and growth rate $\sigma$ of the first harmonic, as obtained solving the self-consistent system (4.11) for increasing values of the amplitude $A$. (b) Same as (a) for the eigenfrequency $\omega_0$ of the leading eigenmode and the oscillation frequency $\omega$ of the first harmonic - Re = 6000.</figcaption></figure>

Typical results obtained solving system (4.11) at Re = 6000 for increasing values of the amplitude are reported in figures 4 and 5. The amplitude of a given component is conveniently measured by the root-mean-square of its fluctuation energy, e.g., $e_1 = \sqrt{2A||u_1||}$ for the first harmonic, $e_2 = \sqrt{2A^2}||u_2||$ for the second harmonic, etc. For A = 0, the mean flow reduces to the unstable base flow $u_b$ and the structure of the first harmonic is given by the eigenmode $\hat{u}_{10b}$, whose growth rate $\sigma_{0b} = 0.55$ and eigenfrequency $\omega_{0b} = 10.69$ are those predicted by the linear theory. When A increases, the leading eigenmode and the orthogonal component grow steadily (with $e_{10} > e_{1\perp}$; see figure 4), and the harmonics follow suit (with $e_1 > e_2$, i.e., the first harmonic dominates even though $e_1 < e_{10}$, which stems from the non-orthogonality of $\hat{u}_{10}$ and $u_{1\perp}$). The growth rates $\sigma_0$ and $\sigma$ concurrently decrease while departing increasingly from each other (with $\sigma < \sigma_0$; see figure 5), meaning that the base flow is progressively distorted into an increasingly stable mean flow. Meanwhile, $\omega_0$ and $\omega$ increase but remain close to each other (with $\omega \lesssim \omega_0$). For $A = 2.87 \times 10^{-2}$, $\sigma$ is exactly zero and the fluctuation (4.2) is purely harmonic, just like the zero-mean, saturated fluctuation of the DNS. The underlying mean flow is strongly linearly unstable with $\sigma_0 = 0.124$, a value close to that $\sigma_{0m} = 0.117$ obtained
A second-order self-consistent model of saturation in the open cavity flow 15
by linear stability analysis of the time-averaged mean flow. The frequency $\omega = 11.21$ of
the first harmonic is very close to that $\omega_{nl} = 11.19$ extracted from the DNS, but so is
the leading eigenfrequency of the mean flow $\omega_0 = 11.23$, as further discussed in § 6.1.

The following picture comes out: in the linear regime, the first harmonic is parallel
to the leading eigenmode $\hat{u}_{10}$. When nonlinearities set in, the self-interaction of $\hat{u}_{10}$
generates a second harmonic $u_2$ whose nonlinear interaction with the eigenmode feeds
back and forces the growth of an orthogonal component $u_{1\perp}$. The latter distorts $\hat{u}_{10}$
into an increasingly modified first harmonic $u_1$, while turning its linear growth rate $\sigma_0$
and eigenfrequency $\omega_0$ into $\sigma$ and $\omega$. Both $u_1$ and $u_2$ feed back on the mean flow via
their Reynolds stresses, and the growth rate of the first harmonic decreases, up to the
point where it becomes exactly zero. At this stage, the base flow has been distorted into
a strongly linearly unstable mean flow ($\sigma_0 > 0$), but perturbations stop growing and
saturate because their linear growth is nonlinearly balanced by the stabilizing effect of
the second harmonic. It is thus the distorted growth rate $\sigma = \sigma_0 + \gamma_r$, not just the
linear growth rate $\sigma_0$, that is relevant to the mechanics of saturation. The oscillation
frequency is similarly given by the distorted frequency $\omega = \omega_0 + \gamma_i$, meaning that the
eigenfrequency $\omega_0$ alone is not expected to be predictive, consistently with the argument
of Sipp & Lebedev (2007).

## 5. Comparison with DNS results

### 5.1. Saturated regime

In the saturated regime, the ensemble average is equivalent to the time-average over
a period of the oscillation, therefore the model predictions can be compared to time-
averaged data from a single DNS. The total, self-consistent oscillation amplitude

$$e = \sqrt{e_1^2 + e_2^2}, \qquad (5.1)$$
<figure><img src="image_8.png" /><figcaption>FIGURE 6. Mean oscillation amplitude e in the saturated regime against Re: self-consistent (blue symbols) vs. DNS (black line).</figcaption></figure>

is shown in figure 6 to be almost identical to its DNS counterpart $||\overline{\mathbf{u}'}||^2$¹/². Similarly, the linear growth rate $\sigma_0$ of the leading eigenmode and the oscillation frequency $\omega$ fall onto the values $\sigma_{0m}$ and $\omega_{nl}$ extracted from the time-averaged mean flow; see figure 7.

The self-consistent mean flow $\langle \mathbf{u} \rangle$ reproduces remarkably well the time-averaged mean flow $\overline{\mathbf{u}}$; see figures 8(a)-(b), and the same striking agreement is observed in figures 8(e)-(f) for the Reynolds stress divergence, calculated either from the model as

$$A^2 \psi(\overline{\mathbf{u}}_{10} + \mathbf{u}_{1\perp}) + A^4 \psi(\mathbf{u}_2), \qquad (5.2)$$

or from the DNS as $\zeta(\overline{\mathbf{u}'})$. We also compare in figures 8(e)-(h) the self-consistent approximations of $\mathbf{u}_1$ and $\mathbf{u}_2$ to those harmonics obtained performing a fast Fourier transform (FFT) of the instantaneous velocity at each point of the computational grid and retaining the terms at frequencies $\omega_{nl}$ and $2\omega_{nl}$. For the first harmonic, the agreement is excellent, both in terms of magnitude and spatial distribution. At first sight, the structure resembles that of the leading eigenmode, whose self-consistent approximation $A\overline{\mathbf{u}}_{10}$ is documented in figure 8(i). There exist subtle differences, however, as the first harmonic precisely differs from the eigenmode by the orthogonal component $A\mathbf{u}_{1\perp}$ shown in figure 8(j), otherwise the cavity flow would comply with the linear marginal stability criterion. Locally, the magnitude of the orthogonal component can reach up to 50% of that of the leading eigenmode, but its structure is more concentrated in the vicinity of the downstream edge,
<figure><img src="image_9.png" /><figcaption>FIGURE 7. Same as figure 6 for the (a) linear growth rate of the leading eigenmode and (b) the oscillation frequency. The results obtained by linear stability analysis of the time-averaged mean flow (resp. of the base flow) are shown by the grey line (resp. the dash-dotted grey line).</figcaption></figure>

hence a ratio $e_{1\perp}/e_{10} \sim 25\%$ of the mean oscillation amplitudes, consistently with the results of figure 4. The spatial distribution of the second harmonic is also well predicted, but we notice a small discrepancy in magnitude, found to be slightly higher using the model. A possible explanation is that the FFT is performed over a time-span of 100 time units to avoid tremendous computational costs. This can alter the estimation of the higher Fourier coefficients, which is supported by the fact that the oscillation amplitude reconstructed from the first six FFT modes is lower by 6% than the value documented in figure 6, obtained averaging the DNS solution over 200 time units.

## 5.2. Transient regime

We compare now the self-consistent predictions to nonlinear data extracted from a stack of DNS, whose initial condition is made up of the base flow $\mathbf{u}_b$ and its unit-norm, leading eigenmode $\hat{\mathbf{u}}_{10b}$ with controlled amplitude $A_0 = 5 \times 10^{-6}$ but arbitrary phase $\phi_0$, i.e.,

$$ \mathbf{u}(t=0, \phi_0) = \mathbf{u}_b + A_0(\hat{\mathbf{u}}_{10b}e^{i\phi_0} + \text{c.c.}). \quad (5.3) $$

At each time instant, ensemble-averaged quantities are obtained averaging over 8 values of $\phi_0$ uniformly distributed in $[0; 2\pi]$. This suffices to converge meaningful averages since it has been checked that identical results are obtained using up to 16 values. It is easy
<figure><img src="image_10.png" /><figcaption>FIGURE 8. (a)-(b) Mean vorticity of the saturated cavity flow, as obtained (a) solving the self-consistent system (4.11) and (b) averaging in time the instantaneous DNS solution. (c)-(d) Same as (a)-(b) for the x-component of the Reynolds stress divergence. (e)-(f) Same as (a)-(b) for the first harmonic $Au_1$. (g)-(h) Same as (e)-(f) for the second harmonic $A^2u_2$. The fields in (f)-(h) have been obtained by fast Fourier transform of the DNS solution and retaining the terms at frequencies $\omega_{nl}$ and $2\omega_{nl}$. (i)-(j) Vorticity of the self-consistent (i) leading eigenmode $A\tilde{u}_{10}$, and (j) orthogonal component $Au_{1\perp}$ - Re = 6000. Only the real part of a complex field is shown.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>FIGURE 9. Transient evolution of the mean oscillation amplitude, as obtained solving the self-consistent system (4.11) for various amplitudes and integrating equations (5.4) in time (blue line) and ensemble averaging the instantaneous DNS solution over the phase of the initial perturbation (black line). The light grey shade marks the transient domain defined as the time-interval for which the amplitude is less than 99% of its saturated value. The grey dash-dotted line corresponds to the linear prediction of exponential growth at rate $\sigma_{0b}$. (a) $Re = 6000$. The yellow symbol marks the time $t = 15.8$ for which figure 10 compares the model and exact mean flows. (b) $Re = 5000$, together with a close-up on the overshoot region.</figcaption></figure>

to verify that the mean oscillation amplitudes $e_1$ and $e_2$ grow in time at rates $\sigma$ and $2\sigma$, respectively. The quasi-static evolution of the self-consistent amplitude is thus determined marching in time equations

$$d_t e_1 = \sigma e_1, \quad \text{and} \quad d_t e_2 = 2\sigma e_2, \qquad (5.4)$$

from the knowledge of $\sigma(A)$, $e_1(A)$ and $e_2(A)$, following which (5.1) is used to recover the total amplitude. For the results to be comparable, we start from initial condition $A(t = 0) = \sqrt{2}A_0$, and proceed with a fourth-order Runge-Kutta scheme, using polynomial interpolation at each time step to calculate $\sigma$ from the three closest neighbours.

Figure 9(a), shows that the obtained self-consistent amplitude matches remarkably well the exact DNS value $\langle\|\mathbf{u}'\|^2\rangle^{1/2}$, up to $t = 20$ where the amplitude has reached 99% of its saturation value. This defines the length of the transient, as indicated by the light grey shade in figure 9(a), and completes proof that the model properly captures the
(a)

(b)

<figure><img src="image_12.png" /><figcaption>FIGURE 10. (a)-(b) Mean vorticity of the saturated cavity flow, as obtained (a) solving the self-consistent system (4.11) for $A = 2 \times 10^{-2}$, and (b) ensemble averaging the instantaneous DNS solution at the corresponding time $t = 15.8$.</figcaption></figure>

main mechanisms responsible for the transient growth and saturation of the instability. During the transient, the flow achieves almost 30 flow oscillations (not shown here), which reasonably advocates the decoupling of scales assumption underlying the derivation of the model. The nonlinearity becomes apparent midterm, at $t \sim 10$, whereupon the self-consistent and DNS data deviate from the linear prediction of exponential growth at rate $\sigma_{0b}$ (grey dash-dotted line). Interestingly, there exist subtle differences in the dynamical paths taken by both solutions. For instance, we show in figure 10 that the self-consistent mean flow obtained for the amplitude $A = 2 \times 10^{-2}$ (that differs from the saturated solution by a lower magnitude of vorticity in the recirculating flow, and a thinner vorticity sheet close to the downstream edge) resembles a lot the DNS mean flow ensemble averaged at the corresponding time instant $t = 15.8$. The model however slightly overestimates the inner vorticity, which appears mandatory to better approximate the structure of the oscillation from the first two harmonics only.

Note, for $20 < t < 30$, the DNS data exhibits an overshoot that the self-consistent model fails to reproduce (probably because the quasi-static assumption is somehow faulted), whereupon the amplitude decreases while slowly oscillating towards its saturated value, reached at $t \sim 200$. Such oscillations may reflect the destabilization of
A second-order self-consistent model of saturation in the open cavity flow 21 additional eigenmodes, as the base cavity flow undergoes a third Hopf bifurcation at $Re_{c3} = 5698$, but similar oscillations show up at $Re = 5000$, albeit with an amplitude so small that they are barely visible in figure 9(b). Another possibility is that the system amplifies the energy of the initial condition by non-normal mechanisms (Chomaz 2005; Schmid 2007), which may cause transient, yet significant mean flow distortions (Pralits, Bottaro & Cherubini 2015). Either ways, the present model is meant to describe the growth of the sole bifurcating eigenmode, and is by construction unable to capture the mode interactions involved in these scenarios.

# 6. Discussion

## 6.1. Connection with weakly nonlinear analysis

There are several similarities (and also key differences) between the present approach and weakly nonlinear analyses meant to figure out the leading-order nonlinear effects associated to small-amplitude disturbances from the canonical Stuart-Landau amplitude

$$d_t A = \lambda \Delta A - (\mu_0 + \mu_2) A |A|^2, \quad (6.1)$$

where A is now the (complex) amplitude of the bifurcating eigenmode, $\Delta = Re_c^{-1} - Re^{-1}$ measures the departure from the instability threshold, $\lambda$ is the (complex) linear growth rate such that the leading eigenvalue is linearly approximated by $i\omega_{0b}(Re_c) + \lambda\Delta$, and $\mu_0$ (resp. $\mu_2$) is the (complex) Landau coefficient determined by the nonlinear interaction of the eigenmode with the mean flow distortion induced by its Reynolds stresses (resp. with the second harmonic). All coefficients can be determined numerically applying multiple time-scale analysis (Sipp & Lebedev 2007; Meliga, Chomaz & Sipp 2009a) at the second threshold of instability threshold $Re_{c2} = 4348$ (since we recall that
the first two eigenmodes to bifurcate exchange dominance at $Re = 4567$). This yields $\lambda = 7689 + 4645i$, $\mu_0 = 583.0 - 129.8i$ and $\mu_2 = 228.8 + 98.98i$.

Weakly nonlinear analysis also invokes a scale separation argument (fast time scale for the flow oscillations vs. slow time scale for the growth of the oscillation amplitude) to derive (6.1) from a hierarchy of equations at various orders in a small parameter, physically representing the order of magnitude of the flow disturbances. However, the mean oscillation amplitude computed from the (real) limit cycle amplitude

$$|A| = \sqrt{\frac{\lambda_r \Delta}{\mu_{0r} + \mu_{2r}}}, \qquad (6.2)$$

departs rapidly from the self-consistent and DNS data for $Re \gtrsim 4600$. This is because the structure of the weakly nonlinear first harmonic is fundamentally flawed far from threshold, as it reduces to the leading eigenmode $\hat{u}_{10b}(Re_c)$ of the marginally stable base flow. This pans out only if the growth rate of unstable disturbances is small enough (i.e., if the Reynolds number is close enough from threshold), otherwise the above results stress the necessity of encompassing the distortion of the leading eigenmode by the orthogonal component, as it is the flow response to this distortion that ultimately selects the oscillation amplitude and frequency. Another glitch concerns the fact that the Reynolds stresses of the weakly nonlinear second harmonic are discarded by virtue of their amplitude being two orders smaller than those of the leading eigenmode in the asymptotic expansion. There is no such hierarchy in self-consistent modeling, where the fluctuation is built as successive-order corrections to the leading eigenmode of the mean flow, and all effects are encompass rigorously regardless of their orders of magnitude, hence the ability to support an arbitrary level of supercriticality.

Valuable insights can still be gained from weakly nonlinear analysis. For instance, the above values of $\mu_{0r}$ and $\mu_{2r}$ indicate that both the mean flow distortion induced by the Reynolds stresses of the leading eigenmode and the second harmonic are stabilizing
<figure><img src="image_13.png" /><figcaption>FIGURE 11. (a) Mean oscillation amplitude in the saturated regime against Re: second-order self-consistent (blue symbols) vs. first-order self-consistent (green symbols) vs. DNS (black line) vs. weakly nonlinear analysis (grey dashed line). (b) Same as (a) for the oscillation frequency. The second-order self-consistent results feature both the nonlinear oscillation frequency ω (blue symbols), and the mean flow eigenfrequency ω₀ (smaller orange symbols). The value obtained by linear stability analysis of the base flow is shown as the grey dash-dotted line.</figcaption></figure>

(μ₀ₖ > 0 and μ₂ₖ > 0), and it can be inferred that both mechanisms contribute to saturate the instability since μ₂ₖ/μ₀ₖ = O(1). This is true also of the present self-consistent results: in the absence of a second harmonic, the stabilizing effect of the Reynolds stresses is evidenced by the mean flow becoming increasingly linearly stable when solving the self-consistent system of Mantič-Lugo et al. (2014) for increasing values of the amplitude Aₛₐ₁ (not shown here). As for the stabilizing effect of the second harmonic, it shows in figure 11(b) through the fact that the amplitude Aₛₐ₂ obtained with the present model is systematically smaller than Aₛₐ₁, just like the amplitude (6.2) is systematically smaller than that |A₀| = √λₜΔ/μ₀ₖ obtained retaining only the mean flow distortion.

The weakly nonlinear limit-cycle frequency is given by

$$ \omega_{wnl} - \omega_{0b} = -\frac{\mu_{0i} + \mu_{2i}}{\mu_{0r} + \mu_{2r}} \lambda_r \Delta = \underbrace{-\mu_{0i}|A_0|^2}_{(i)} - \underbrace{-\mu_{0i}(|A|^2 - |A_0|^2)}_{(ii)} - \underbrace{-\mu_{2i}|A|^2}_{(iii)}, \quad (6.3) $$

where the nonlinear frequency correction in the RHS brings out three contributions: a first one (i) due to the mean flow distortions induced by the Reynolds stresses of the leading eigenmode, a second one (ii) due to the second harmonic changing the amplitude
of the Reynolds stresses from $|A_0|^2$ to $|A|^2$, and a third one (iii) due to the nonlinear interaction between the second harmonic and the leading eigenmode. The above values of $\mu_{0i}$ and $\mu_{2i}$ suggest that the mean flow distortion increases the frequency ($\mu_{0i} < 0$), while the stabilizing effect of the second harmonic tends to decrease it ($|A| < |A_0|$), and the nonlinear interaction between harmonics yields additional frequency reduction ($\mu_{2i} > 0$). Quite fortuitously, all three effects balance almost exactly each other because $|\mu_{0i} + \mu_{2i}|/|\mu_{0r} + \mu_{2r}| \ll 1$, hence the weakly nonlinear frequency is almost identical to the leading eigenfrequency of the base flow $\omega_{0b}$, as can be seen comparing the grey dashed and dash-dotted lines in figure 11(b) (we insist that the nonlinearity itself is significant, otherwise the weakly nonlinear frequency would fall onto the self-consistent and DNS data regardless of the value of Re). These interpretations essentially carry over to the strongly nonlinear regime provided the self-consistent frequency $\omega_{sc2}$ is recast into

$$ \omega_{sc2} - \omega_{0b} = (\omega_{0sc2} - \omega_{0sc1}) + (\omega_{0sc1} - \omega_{0b}) + \gamma_i, \quad (6.4) $$

where $\omega_{0sc1} = \omega_{sc1}$ is the mean flow eigenfrequency predicted by the first-order model of Mantič-Lugo et al. (2014). The mean flow distorsion induced by the Reynolds stresses of the leading eigenmode is seen in figure 11(b) to increase the oscillation frequency ($\omega_{0sc1} > \omega_{0b}$). The modification of the Reynolds stresses induced by the second harmonic conversely decreases it ($\omega_{0sc1} > \omega_{0sc2}$), which is similar to the interpretation of the Landau coefficients, although we reiterate that the physics captured by weakly nonlinear analysis on this is very limited. Finally, the nonlinear interaction of the harmonics with the leading eigenmode yields additional frequency reduction for Re < 6500 ($\omega_{sc2} \lesssim \omega_{0sc2}$) but conversely an additional frequency increase for Re > 6500 ($\omega_{sc2} \gtrsim \omega_{0sc2}$). This is barely visible in figure 11(b) because of the smallness of the frequency detuning $\gamma_i$, that explains why the mean flow eigenfrequency ends up doing well on the frequency prediction. Anyhow, this is likely to be fortuitous, as Turton et al. (2015) report specific
A second-order self-consistent model of saturation in the open cavity flow 25
regimes of thermosolutal convection in which the eigenfrequency is not predictive, which
we infer is because the related value of $\gamma_i$ is $O(1)$.

The nonlinearity at play in given flow can be identified from the linear stability prop-
erties of the time-averaged mean flow: a **strongly linearly unstable** leading eigenmode
or one whose eigenfrequency approximates poorly the nonlinear frequency reflects both
**mean flow distortions** and generation of higher-order harmonics. The main features of
the related oscillation are best predicted using a multi-oder model to encompass the $O(1)$
value of the underlying detuning parameter. A leading eigenmode close to marginal linear
stability and whose eigenfrequency approximates well the nonlinear frequency conversely
reflects **pure mean flow distortions**, in which case relevant predictions can be obtained
with the model of Mantič-Lugo et al. (2014) because the detuning parameter is close to
zero. Both conditions are required, as the present cavity flow serves as a **reminder** that
the mean flow eigenfrequency can quite fortuitously predicts the nonlinear frequency
(if $\gamma_i \ll 1$) even though the nonlinearity does not restrict to mean flow distortion (if
$\gamma_r = O(1)$). This **generalizes** the distinction introduced by Sipp & Lebedev (2007) from
their analysis of the amplitude equation, the weakly nonlinear equivalent of the present
model being an amplitude equation with $|\mu_2|/|\mu_0| = O(1)$, while that of the model
of Mantič-Lugo et al. (2014) is an amplitude equation with $|\mu_2|/|\mu_0| \ll 1$.

6.2. *Effect of the truncation order*

Caution should be used in generalizing the present results to other cases involving signif-
icant generation of higher harmonics. If the flow is not within the scope of the first-order
model, then we expect a valuable (albeit incomplete) description of the main nonlinear
mechanisms at play to be obtained from a second-order model. Still, it may be necessary
to resort to a higher-order of Fourier expansion of the fluctuation to obtain more accu-
rate, quantitative predictions, which requires to augment the model with linear equations
<figure><img src="image_14.png" /><figcaption>FIGURE 12. (a) Amplitude $e_j$ of the first six harmonics, as obtained by fast Fourier transform of the DNS solution. (b) Contribution $\gamma_j$ of each mode to the detuning parameter, as defined by equation (6.5) - Re = 6000. All results normalized for the second harmonic to peak at unity.</figcaption></figure>

for the higher harmonics forced by relevant nonlinear terms describing the interactions
of the lower harmonics. In this eventuality, an interesting point to discuss is whether or
not it is possible to anticipate the proper number of harmonics to be retained. Spectra
can be useful to unravel the existence (or the absence) of harmonics, but they are not
very insightful for the intended purpose because they give only a local, partial vision of
the nonlinearity. For instance velocity spectra measured in the recirculating flow unveil
essentially a single peak at the fundamental frequency (not shown here), which merely
reflects the low magnitude of the second harmonic in this specific flow region; see figure 8.
While a more global perspective is achieved resorting to FFT to extract the structure
of the various harmonics (as has been done for the first and second harmonics in § 5.1)
it may not be relevant to rest on the sole oscillation amplitude $e_j = \|\mathbf{u}_j\|$, as we show
in figure 12(a) that the amplitude of the third (and even the fourth) mode of the cavity
is not necessarily obviously negligible. This is because the nonlinear correction induced
by the higher-order harmonic via the detuning parameter $\gamma$ depends not only on the
structure of the harmonics, but also on their recombination via the advection operator
(that yields the forcing terms acting at each frequency) and on the orientation of these
forces with respect to the adjoint eigenmode. Assuming that $n$ harmonics are available
$$
\gamma = \sum_{j=2}^{n} \gamma_j, \quad \text{with} \quad \gamma_j = - \frac{(\hat{\mathbf{u}}_{10}^{\dagger} | \phi(\mathbf{u}_{j-1}^{*}, \mathbf{u}_{j}))}{(\hat{\mathbf{u}}_{10}^{\dagger} | \mathbf{u}_{1})}, \qquad (6.5)
$$

where each individual detuning $\gamma_j$ accounts for additional distortions to the linear growth
rate and eigenfrequency of the leading eigenmode induced by the $j^{th}$ harmonic. We have
computed the individual detuning coefficients of the first six FFT modes, using the
leading adjoint eigenmode of the time-averaged mean flow to evaluate numerically all
inner products. The results shown in figure 12(b) clearly stress that only the contribution
$\gamma_2$ of the second harmonic needs to be retained, with marginal correction of the third
harmonic, and no effect of the higher-order harmonics. This is as close as one can get to
identify a relevant order of truncation without any input from the self-consistent model.

As a proof of feasibility (and also as an attempt to assess the relevance of the above
criterion), we show in figure 13 that all values obtained solving the third-order self-
consistent model are almost identical to those presented hereinabove, hence confirming
that the second-order results come with little to no truncation effect. Of course, the
computational cost increases with the order of the model because an inner loop needs
to be added in the iterative algorithm for the computation of each additional harmonics
(hence three nested loops for the third-order model, and so on). For our case, converging
the third harmonic to the same precision yields an increase in the number of iterations
by roughly 30%. Finding ways to improve the numerical resolution of the self-consistent
equations in the foreseeable future is thus a major stake, as increasing the number of
harmonics may be a much needed step towards extending the scope of the method to
more complex turbulent regimes.
<figure><img src="image_15.png" /><figcaption>FIGURE 13. (a) Mean oscillation amplitude in the saturated regime against $Re$: second-order self-consistent (blue symbols) vs. third-order self-consistent (red symbols) vs. DNS (black line). (b)-(c) Same as (a) for (b) the linear growth rate of the leading eigenmode and (c) the oscillation frequency. The grey line in (b) is the value obtained by linear stability analysis of the time-averaged mean flow.</figcaption></figure>

### 6.3. Concluding remarks

The present model provides a consistent picture of the saturation process in unsteady flows whose nonlinearity involves not only mean flow distortions but also production of higher harmonics. Saturation occurs when the growth rate of the first harmonic is zero, at which point the second harmonic balances exactly the non-zero linear growth rate of the leading eigenmode. This rationalizes the existence of strongly linearly unstable mean flows, and paves the way for new analysis techniques of such flows, including model reduction and flow control, just like the model of Mantič-Lugo et al. (2014) did for systems whose mean flow is linearly marginally stable (Meliga et al. 2016a). It is worth noticing that the relevance of this type of models, in which perturbation equations are linearized around a mean flow determined simultaneously from the Reynolds stress of the perturbations, goes well beyond limit cycles. In particular, we may hope that the method could be generalized to address extrinsic unsteady dynamics, for which a first step has been taken by Mantič-Lugo & Gallaire (2016) using their first-order model. Also, several papers in the recent literature have used similar coupled models retaining the Reynolds
A second-order self-consistent model of saturation in the open cavity flow 29
stresses of transiently unstable linear and weakly nonlinear disturbances to explain the
emergence of coherent structure in otherwise linearly stable, parallel shear flows (Farrell
& Ioannou 2012; Pralits et al. 2015; Beaume, Chini, Julien & Knobloch 2015; Farrell,
Ioannou, Jiménez, Constantinou, Lozano-Durán & Nikolaidis 2016).

This work is supported by the «Investissements d'Avenir» French Government pro-
gram, managed by the French National Research Agency (ANR) through the A*MIDEX
grant (ANR-11-IDEX-0001-02) and the LABEX MEC project (ANR-11-LABX-0092).

REFERENCES

ÅKERVIK, E., EHRENSTEIN, U., GALLAIRE, F. & HENNINGSON, D.S. 2008 Global two-dimensional stability measures of the flat plate boundary-layer flow. *Eur. J. Mech. B-Fluid* **27** (5), 501–513.

BARBAGALLO, A., SIPP, D. & SCHMID, P.J. 2009 Closed-loop control of an open cavity flow using reduced-order models. *J. Fluid Mech.* **641**, 1–50.

BARKLEY, D. 2006 Linear analysis of the cylinder wake mean flow. *Europhys. Lett.* **75**, 750–756.

BEAUME, C., CHINI, G.P., JULIEN, K. & KNOBLOCH, E. 2015 Reduced description of exact coherent states in parallel shear flows. *Phys. Rev. E* **91**, 014101.

CHOMAZ, J.-M. 2005 Global instabilities in spatially developing flows: Non-normality and non-linearity. *Annu. Rev. Fluid Mech.* **37**, 357–392.

DUŠEK, J., LE GAL, P. & FRAUNIÉ, P. 1994 A numerical and theoretical study of the first Hopf bifurcation in a cylinder wake. *J. Fluid Mech.* **264**, 59–80.

FARRELL, B.F. & IOANNOU, P.J. 2012 Dynamics of streamwise rolls and streaks in turbulent wall-bounded shear flow. *J. Fluid Mech.* **708**, 149–196.

FARRELL, B.F., IOANNOU, P.J., JIMÉNEZ, J., CONSTANTINOU, N.C., LOZANO-DURÁN, A. & NIKOLAIDIS, M.-A. 2016 A statistical state dynamics-based study of the structure and mechanism of large-scale motions in plane poiseuille ?ow. *J. Fluid Mech.* **809**, 1–26.

FORNBERG, B. 1980 A numerical study of steady viscous flow past a circular cylinder. *J. Fluid Mech.* **98**, 819–855.
HAMMOND, D. A. & REDEKOPP, L. G. 1997 Global dynamics of symmetric and asymmetric wakes. *J. Fluid Mech.* **331**, 231–260.

MALKUS, W. V. R. 1956 Outline of a theory of turbulent shear flow. *J. Fluid Mech.* **1**, 521–539.

MANTIČ-LUGO, V., ARRATIA, C. & GALLAIRE, F. 2014 Self-consistent mean flow description of the nonlinear saturation of the vortex shedding in the cylinder wake. *Phys. Rev. Lett.* **113**, 084501.

MANTIČ-LUGO, V., ARRATIA, C. & GALLAIRE, F. 2015 A self-consistent model for the saturation dynamics of the vortex shedding around the mean flow in the unstable cylinder wake. *Phys. Fluids* **27**, 074103.

MANTIČ-LUGO, V. & GALLAIRE, F. 2016 Self-consistent model for the saturation mechanism of the response to harmonic forcing in the backward-facing step flow. *J. Fluid Mech.* **793**, 777–797.

MAUREL, A., PAGNEUX, V. & WESFREID, J.E. 1995 Mean-flow correction as non-linear saturation mechanism. *Europhys. Lett.* **32**, 217–222.

MELIGA, P., BOUJO, E. & GALLAIRE, F. 2016a A self-consistent formulation for the sensitivity analysis of finite amplitude vortex shedding in the cylinder wake. *J. Fluid Mech.* **800**, 327–357.

MELIGA, P., BOUJO, E., PUJALS, G. & GALLAIRE, F. 2014 Sensitivity of aerodynamic forces in laminar and turbulent flow past a square cylinder. *Phys. Fluids* **26**, 104101.

MELIGA, P., CADOT, O. & SERRE, E. 2016b Experimental and theoretical sensitivity analysis of turbulent flow past a square cylinder. *Flow Turbul. Combust.* **97**, 987–1015.

MELIGA, P., CHOMAZ, J.-M. & SIPP, D. 2009a Global mode interaction and pattern selection in the wake of a disk: a weakly nonlinear expansion. *J. Fluid Mech.* **633**, 159–189.

MELIGA, P., PUJALS, G. & SERRE, E. 2012 Sensitivity of 2-D turbulent flow past a D-shaped cylinder using global stability. *Phys. Fluids* **24**, 061701.

MELIGA, P., SIPP, D. & CHOMAZ, J.-M. 2009b Elephant modes and low frequency unsteadiness in a high Reynolds number, transonic afterbody wake. *Phys. Fluids* **21**, 054105.

NOACK, B.R., AFANASIEV, K., MORZYNSKI, M., TADMOR, G. & THIELE, F. 2003 A hierarchy
*A second-order self-consistent model of saturation in the open cavity flow* 31
of low-dimensional models for the transient and post-transient cylinder wake. *J. Fluid Mech.* **497**, 335–363.

PIER, B. 2002 On the frequency selection of finite-amplitude vortex shedding in the cylinder wake. *J. Fluid Mech.* **458**, 407–417.

PRALITS, J.O., BOTTARO, A. & CHERUBINI, S. 2015 Weakly nonlinear optimal perturbations. *J. Fluid Mech.* **785**, 135–151.

SCHMID, P.J. 2007 Nonmodal stability theory. *Annu. Rev. Fluid Mech.* **39**, 129–162.

SIPP, D. & LEBEDEV, A. 2007 Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows. *J. Fluid Mech.* **593**, 333–358.

STUART, J.T. 1971 Non-linear stability theory. *Annu. Rev. Fluid Mech.* **3**, 347–370.

TURTON, S.E., TUCKERMAN, L.S. & BARKLEY, D. 2015 Prediction of frequencies in ther-
mosolutal convection from mean flows. *Phys. Rev. E* **91**, 043009.

ZIELINSKA, B.J.A., GOUJON-DURAND, S., DUŠEK, J. & WESFREID, J.E. 1997 Strongly non-linear effect in unstable wakes. *Phys. Rev. Lett.* **79**, 3893–3896.