Extended mean flow analysis
of the circular cylinder flow

A new eigenvalue analysis is developed and applied to the circular cylinder laminar flow configuration to investigate the various mechanisms at play in the nonlinear saturation of perturbations yielding to limit cycles for supercritical values of the Reynolds number. Unlike the mean-flow analysis, which only accounts for the interaction of the first-harmonic of the time-periodic flow with its mean-flow, the so-called extended mean-flow analysis also accounts for its interaction with the second-harmonic. The results reveal the existence of two eigenvalues both having a growth rate exactly equal to zero. The high-frequency eigenmode gives better estimations of the frequency and spatial structure of the first-harmonic than the marginal modes obtained with the mean-flow analysis, especially when the Reynolds number is increased.

Self-sustained oscillations occur in open flows when an infinitesimal global perturbation gets unstable and then saturates to a limit-cycle state. Although the linear stability analysis of the underlying steady base flow provides a rigorous mathematical description of the initial instability growth mechanisms, its predictive capabilities rapidly reduce as flow oscillations grow in amplitude away from the bifurcation. Corresponding leading-order nonlinear effects can be taken into account by introducing a weakly nonlinear formalism, which however remains still limited by the perturbative nature of the involved approximation. A classical example is represented by the onset of the Bérnard-von Kármán vortex street in the wake of a circular cylinder[1, 2, 3] at Reynolds number $\textit{Re}\approx 47$. For this flow, the frequency drift from the linear stability estimate[4, 5] rapidly increases up to $30\%$ at $\textit{Re}=80$, while the error on the oscillation amplitude computed by a weakly nonlinear analysis[6] exceeds $100\%$, with the convergence radius of the corresponding normal-form series being vanishingly small[7].

The possibility to better account for the nonlinear effects and then for the saturation mechanism of the flow oscillations, obtaining an accurate prediction of their frequency and spatial structure, has moved the attention of the researchers to the stability properties of the resulting time-averaged mean flow. Indeed for the laminar cylinder wake[4] and other different flow configurations, either in the turbulent or laminar regime[8, 9, 10], the linear stability analysis heuristically performed on the top of the mean flow yields a mildly unstable eigenvalue at approximately the same frequency of the nonlinear unsteady flow, a result which is reminiscent of the seminal works by Malkus[11] and Stuart[12]. These authors first emphasized the role of the mean-flow distortion in the nonlinear saturation of linear perturbations. As disturbances, initially fed on the unstable base flow, linearly grow, reaching a finite-size amplitude, the resulting Reynolds stresses induce a distortion of the mean flow, which in turn reduces the temporal growth rate of the perturbations until the growth-rate is nearly zero. This conjecture leads to invoke the marginal stability property of the mean flow, for which the name ‘real-zero imaginary frequency’ (RZIF) has been introduced by Turton, Tuckerman and Barkley[10], showing that it exactly holds for pure monochromatic flow oscillations, and it is approximately satisfied when higher-order harmonics becomes negligible compared to the leading one, such as in the case of travelling waves in thermosolutal convection.

The RZIF property has been exploited by Mantič-Lugo, Arratia and Gallaire[6] to derive a self-consistent model of the saturation dynamics of the cylinder wake, by formalizing the above described phenomenological picture. In their model the nonlinear fluctuation is represented by the mean flow leading global mode, with the mean flow being determined as the solution of the steady Navier-Stokes equations forced by the Reynolds stresses resulting from the interaction of the leading mode with itself. The arbitrary mode amplitude is computed so that the mean flow distortion yields a zero growth-rate, thus assuming the RZIF property to close the model. The authors show that their model accurately predicts all the relevant features of the flow obtained from direct numerical simulations (DNS) up to $\textit{Re}=110$, including flow oscillation frequency, amplitude and spatial pattern, as well as the mean flow and the Reynolds stresses. Later this model has been employed by Meliga, Boujo and Gallaire[13] to perform an adjoint sensitivity analysis of the limit cycle frequency and amplitude, in order to identify the wavemaker regions of the nonlinear wake oscillations and the most sensitive regions for passive control design. However, despite these remarkable results, almost exclusively obtained for the cylinder wake, the RZIF property does not hold for all the flows. Well known counterexamples are represented by the oscillations in an open cavity[5], in turbulent wakes[14, 15] and and by standing waves in thermosolutal convection[10]. Indeed the feedback between the mean flow distortion and the leading mode amplification is not the unique mechanism explaining the nonlinear saturation but another fundamental mechanism exists, which is driven by the production of high-order harmonics. For this enriched dynamics, a strong deviation from the marginal stability of the mean flow is often observed, and even the frequency prediction is not always accurate. As an example in the case of an open cavity, Sipp & Lebedev[5] pointed out the role played by the interaction between the leading global mode and its second harmonic, resulting in an unstable mean flow. The aim of this paper is to introduce a new analysis that extends the mean-flow analysis as it takes into account both the effects of the mean flow distortion and the nonlinear interaction with the second harmonic.

The paper is organized as follow. The theoretical formulation is introduced in 2 including the eigenvalue analysis of the base flow, the mean-flow analysis and the extended mean-flow analysis. This new analysis is then applied to the circular cylinder flow configuration and results are presented 3

The dynamics of laminar unsteady hydrodynamic flows in a domain $\Omega$ is governed by the following non-dimensional incompressible Navier-Stokes equations

$$\mathcal{M}\frac{\partial\mathbf{q}}{\partial t}+\mathcal{N}(\mathbf{q})=\mathbf{0},$$

(1)

where $\mathbf{q}=(\mathbf{u},p)^{T}$ denotes the flow state composed of the flow velocity field $\mathbf{u}(\mathbf{x},t)$ and the pressure field $p(\mathbf{x},t)$ that depend on the space coordinates $\mathbf{x}$ and on the time variable $t$. The linear operator $\mathcal{M}$ is defined by its action on the flow state as $\mathcal{M}\mathbf{q}=(\mathbf{u},^{0})^{T}$ and the nonlinear residual is defined as

$$\mathcal{N}(\mathbf{q})=\begin{pmatrix}(\mathbf{u}\cdot\mathbf{\nabla})\mathbf{u}+\mathbf{\nabla}{p}-\textit{Re}^{-1}\mathbf{\Delta}\mathbf{u}\\[5.69054pt] \mathbf{\nabla}\cdot\mathbf{u}\end{pmatrix},$$

(2)

where the Reynolds number $\textit{Re}=U_{\infty}D/\nu$ is based on a reference velocity scale $U_{\infty}$ and a reference length scale $L$ used to make all the flow variables dimensionless. In the cylinder flow configuration investigated in the next section, they will correspond to the uniform inflow velocity and the cylinder diameter $L=D$.

Time-independent solutions of 1, often called base flows and denoted hereinafter $\mathbf{q}_{b}=(\mathbf{u}_{b},p_{b})^{T}$, satisfy the steady non-linear equations

$$\mathcal{N}(\mathbf{q}_{b})=\mathbf{0}.$$

(3)

Their temporal stability is investigated by superimposing an infinitesimal perturbation as

$$\mathbf{q}(\mathbf{x},t)=\mathbf{q}_{b}(\mathbf{x})+\epsilon\left(\hat{\mathbf{q}}(\mathbf{x})e^{\left(\sigma_{b}+\rm{i}\omega_{b}\right)t}+\mbox{ c.c.}\right)+O(\epsilon^{2})$$

(4)

where $\epsilon$ is an infinitesimally small number and c.c. is an abbreviation of complex conjugate. The linear perturbation is decomposed into a complex-valued spatial structure $\hat{\mathbf{q}}$ growing or decaying exponentially in time at the growth rate $\sigma_{b}$ and oscillating at the frequency $\omega_{b}$. Introducing the above decomposition into the non-linear 1, using 3 of the base-flow and neglecting higher-order terms $O(\epsilon^{2})$, yields the following generalized eigenvalue problem

$$\left[\left(\sigma_{b}+\rm{i}\omega_{b}\right)\mathcal{M}+\mathcal{L}(\mathbf{q}_{b})\right]\hat{\mathbf{q}}=\mathbf{0}$$

(5)

where $\mathcal{L}(\mathbf{q}_{b})$ is obtained by linearizing the non-linear residual $\mathcal{N}(\mathbf{q})$ around the base-flow $\mathbf{q}_{b}$. For the incompressible Navier-Stokes equations, the so-called Jacobian operator is defined as

$$\mathcal{L}(\mathbf{q}_{b})=\left(\begin{tabular}[]{lr}$(\mathbf{u}_{b}\cdot\mathbf{\nabla})+(\mathbf{\nabla}\mathbf{u}_{b})-\textit{Re}^{-1}\mathbf{\Delta}$&$\nabla$\\[5.69054pt] $\mathbf{\nabla}\cdot$&$0$\end{tabular}\right)$$

(6)

The stability of the base-flow is determined by computing the leading eigenvalue, i.e. the eigenvalue having the largest growth rate. When its growth rate is positive ($\sigma_{b}>0$), the base-flow is unstable. If the frequency of the unstable eigenvalue is different from zero ($\omega_{b}\neq 0$), the base-flow evolves towards a limit cycle which is an unsteady solution of 1 satisfying the periodicity condition

$$\mathbf{q}(\mathbf{x},t+T)=\mathbf{q}(\mathbf{x},t)$$

(7)

where $T$ is the period. The frequency of the limit cycle $\omega=2\pi/T$ is a priori different from the frequency of the unstable mode ($\omega\neq\omega_{b}$), that will be called theireafter the base-flow frequency. The limit-cycle solution is decomposed as

$\displaystyle\mathbf{q}(\mathbf{x},t)=\mathbf{q}_{0}(\mathbf{x})+\mathbf{q}^{\prime}(\mathbf{x},t)$

(8)

the sum of a time-independent component $\mathbf{q}_{0}(\mathbf{x})$, called the mean-flow, and a time-dependent component $\mathbf{q}^{\prime}(\mathbf{x},t)$ called the fluctuation. This fluctuation is time-periodic and can thus be expanded as a series of harmonics as

$\displaystyle\mathbf{q}^{\prime}(\mathbf{x},t)=\left(A\mathbf{q}_{1}(\mathbf{x})e^{\rm{i}\omega t}+A^{2}\mathbf{q}_{2}(\mathbf{x})e^{2\rm{i}\omega t}+\cdots\right)+\mbox{ c.c.},$

where $\mathbf{q}_{1}(\mathbf{x})$ is the (complex-valued) first harmonic oscillating at the fundamental frequency $\omega$ and $\mathbf{q}_{2}(\mathbf{x})$ is the second harmonic oscillating at twice the fundamental frequency. In the following, higher-harmonics are omitted unless stated otherwise. The first and second harmonics are scaled by the amplitude $A$ which is a complex number, implictely defined by the following normalization condition

$$\int_{\Omega_{s}}\mathbf{q}_{1}^{*}(\mathbf{x})\mathbf{q}_{1}(\mathbf{x})\,d\mathbf{x}=1\,.$$

(9)

The spatial domain $\Omega_{s}$ is enclosed in the computational domain $\Omega$ and allows a definition of the amplitude $A$ independent of the computational domain. The mean flow, first and second harmonics can be explicitely defined from the knowledge of the instantaneous solution $\mathbf{q}(t)$ over a period $[t_{0},t_{0}+T]$ as

$$A^{n}\mathbf{q}_{n}=\frac{1}{T}\int_{t_{0}}^{t_{0}+T}\mathbf{q}(t)e^{\rm{i}n\omega t}dt\;,$$

(10)

the instantaneous solution being computed by marching in time 1 and 2.

To better characterize the different non-linear mechanisms at play in a limit-cycle, it is interesting to formulate the governing non-linear equations in the frequency domain rather than in the time domain. To that aim, the truncated expansion 2 is introduced into 1. By balancing the time-independent terms and those oscillating at frequencies $\omega$ and $2\omega$, one obtains the so-called harmonic balanced formulation of the governing non-linear equation 1, truncated at the second order,

	$\displaystyle\mathcal{N}(\mathbf{q}_{0})=|A|^{2}\,\mathcal{C}(\mathbf{q}_{1}^{*})\mathbf{q}_{1}$	$\displaystyle+$	$\displaystyle|A|^{4}\,\mathcal{C}(\mathbf{q}_{2}^{*})\mathbf{q}_{2},$		(11)
	$\displaystyle\left[\rm{i}\omega\mathcal{M}+\mathcal{L}(\mathbf{q}_{0})\right]\mathbf{q}_{1}$	$\displaystyle=$	$\displaystyle|A|^{2}\,\mathcal{C}(\mathbf{q}_{2})\mathbf{q}_{1}^{*},$		(12)
	$\displaystyle\left[2\rm{i}\omega\mathcal{M}+\mathcal{L}(\mathbf{q}_{0})\right]\mathbf{q}_{2}$	$\displaystyle=$	$\displaystyle\frac{1}{2}\,\mathcal{C}(\mathbf{q}_{1})\mathbf{q}_{1},$		(13)

where $(\cdot)^{*}$ denotes the conjugation operation. The quadratic interaction between an harmonic $\mathbf{q}_{n}$ and an $\mathbf{q}_{m}$ is denoted as

$$\mathcal{C}(\mathbf{q}_{n})\mathbf{q}_{m}=\begin{pmatrix}-(\mathbf{u}_{n}\cdot\mathbf{\nabla})\mathbf{u}_{m}-(\mathbf{u}_{m}\cdot\mathbf{\nabla})\mathbf{u}_{n}\\ 0\end{pmatrix}$$

(14)

Note that the sign minus in the above definition. The mean flow is governed by 11 similar to 3 that governs the base flow, except for the right-hand side terms that are the quadratic interactions between the harmonics $\mathbf{q}_{n}$ and their complex conjugate $\mathbf{q}_{n}^{*}$, where $n=1,2$ because of the second-order expansion. The sum of all these terms is the divergence of the Reynolds-stress tensor. When the amplitude of the fluctuation is zero, i.e. $|A|=0$, it vanishes and the mean-flow equation is exactly the base-flow equation. 12 and 13 govern the first and second harmonics, respectively. The linear operators on the left-hand sides involve the Jacobian operator Refeq:Lns-def defined not with the base-flow $\mathbf{q}_{b}$ but with the mean-flow $\mathbf{q}_{0}$. The right-hand side terms couple the two equations.

Let us now focus on 12 that governs the first harmonic $\mathbf{q}_{1}$. As explained in [10], the right-hand side term $|A|^{2}\,\mathcal{C}(\mathbf{q}_{2})\mathbf{q}_{1}^{*}$ is negligible if the periodic flow is purely monochromatic ($\mathbf{q}_{2}=\mathbf{0}$) or if the amplitude of the first-harmonic is small ($|A|\ll 1$). In those two cases, 12 can be approximated by

$$\left[\rm{i}\omega\mathcal{M}+\mathcal{L}(\mathbf{q}_{0})\right]\mathbf{q}_{1}\approx\mathbf{0}$$

(15)

Therefore, when examining the eigenvalue spectrum of the Jacobian operator around the mean flow $\mathcal{L}(\mathbf{q}_{0})$ after solving the generalized eigenvalue problem

$$\left[\left(\sigma_{m}+\rm{i}\omega_{m}\right)\mathcal{M}+\mathcal{L}(\mathbf{q}_{0})\right]\hat{\mathbf{q}}_{m}=\mathbf{0}$$

(16)

one should find an eigenvalue whose growth rate $\sigma_{m}$ is close to zero, frequency $\omega_{m}$ is close to the fundamental frequency $\omega$ of the time-periodic flow and associated eigenmode $\hat{\mathbf{q}}_{m}$ approximates the spatial structure of the first harmonic $\mathbf{q}_{1}$. This property of limit cycles is called the Real Zero Imaginary Frequency (RZIF) in [10] and can be sum up as

$$\sigma_{m}\approx 0\;\;,\;\;\omega_{m}\approx\omega\;\;,\;\;\hat{\mathbf{q}}_{m}\approx\mathbf{q}_{1}\;.$$

(17)

It was first observed in [4] for the wake flow behind a circular cylinder flow in the range $50<\textit{Re}<100$. However, this property of limit cycles is not necessarily satisfied, as shown in [5] for an open-cavity flow configuration. In that case, the quadratic interaction term $(\mathcal{C}(\mathbf{q}_{2})\mathbf{q}_{1}^{*})$ in 12 is not negligible. Using a weakly non-linear expansion of the flow to compute the coefficients of the equation governing the temporal evolution of the amplitude $A$, Sipp & Lebedev [5] quantified the role of the mean-flow and second-harmonic in the nonlinear saturation of the amplitude. However, the validity of the weakly non-linear analysis is limited to small amplitudes that are reached close the bifurcation threshold.

To circumvent such limitation, we introduce a new eigenvalue analysis that allows to account for both effects, unlike the mean-flow analysis. The complex-valued first and second harmonics are first decomposed into their real and imaginary part as, $\mathbf{q}_{1}=\mathbf{q}_{1}^{r}+\rm{i}\mathbf{q}_{1}^{i}$ and $\mathbf{q}_{2}=\mathbf{q}_{2}^{r}+\rm{i}\mathbf{q}_{2}^{i}$. By introducing these decompositions in 12 and splitting the imaginary and real contributions, one obtains the coupled system of equations

	$\displaystyle\omega\mathcal{M}\mathbf{q}_{1}^{r}+\mathcal{L}\left[\mathbf{q}_{0}+|A|^{2}\,\mathbf{q}_{2}^{r}\right]\mathbf{q}_{1}^{i}-|A|^{2}\,\mathcal{C}(\mathbf{q}_{2}^{i})\mathbf{q}_{1}^{r}$	$\displaystyle=$	$\displaystyle 0$
	$\displaystyle-\omega\mathcal{M}\mathbf{q}_{1}^{i}+\mathcal{L}\left[\mathbf{q}_{0}-|A|^{2}\,\mathbf{q}_{2}^{r}\right]\mathbf{q}_{1}^{r}-|A|^{2}\,\mathcal{C}(\mathbf{q}_{2}^{i})\mathbf{q}_{1}^{i}$	$\displaystyle=$	$\displaystyle 0$

governing the real-valued variables $\mathbf{q}_{1}^{r}$ and $\mathbf{q}_{1}^{i}$. By gathering them into the (real-valued) extended vector $\mathbf{Q}_{1}=(\mathbf{q}_{1}^{r},\mathbf{q}_{1}^{i})^{T}$ and multiplying the second equation by $(-1)$, the above system is rewritten

$\displaystyle\left(\omega\mathcal{M}_{e}+\mathcal{L}_{e}(\mathbf{q}_{0},|A|^{2}\mathbf{q}_{2})\right)\mathbf{Q}_{1}=\mathbf{0}$

(18)

where $\mathcal{M}_{e}=\begin{bmatrix}\mathcal{M}&0\\[5.69054pt] 0&\mathcal{M}\end{bmatrix}$ and $\mathcal{L}_{e}$ is the extended mean-flow operator defined as

$\displaystyle\mathcal{L}_{e}(\mathbf{q}_{0},|A|^{2}\mathbf{q}_{2})=\begin{bmatrix}-|A|^{2}\,\mathcal{C}(\mathbf{q}_{2}^{i})&\mathcal{L}\left[\mathbf{q}_{0}+|A|^{2}\mathbf{q}_{2}^{r}\right]\\[5.69054pt] -\mathcal{L}\left[\mathbf{q}_{0}-|A|^{2}\mathbf{q}_{2}^{r}\right]&|A|^{2}\,\mathcal{C}(\mathbf{q}_{2}^{i})\end{bmatrix}$

that depends not only on the mean-flow $\mathbf{q}_{0}$ but also on the second harmonic $|A|^{2}\mathbf{q}_{2}$. Assuming they are both known (computed using 10 as explained before), the following generalized eigenvalue problem can be solved

$$\left[\left(\omega_{e}+\rm{i}\sigma_{e}\right)\mathcal{M}_{e}+\mathcal{L}_{e}(\mathbf{q}_{0},|A|^{2}\mathbf{q}_{2})\right]\mathbf{\hat{Q}}_{e}=\mathbf{0}$$

(19)

where $\mathbf{\hat{Q}}_{e}=(\mathbf{\hat{q}}_{e}^{1},\mathbf{\hat{q}}_{e}^{2})^{T}$ is a complex-valued eigenvector associated to a complex eigenvalue $\omega_{e}+\rm{i}\sigma_{e}$. In the following, they will be called the extended eigenmode and the extended eigenvalue, respectively. Note that the real part of the extended eigenvalue corresponds to the physical frequency while the imaginary part corresponds to the growth rate. Comparing 19 with 18, one expects the existence of an extended eigenvalue, whose imaginary part is equal to zero ($\sigma_{e}=0$) while its real part is equal to the fundamental frequency ($\omega_{e}=\omega$), and the real corresponding eigenvector approximate the first harmonic ($\mathbf{\hat{Q}}_{e}=\mathbf{Q}_{1}$). This property of limit cycles is called the Extended Real Zero Imaginary Frequency (ERZIF) and can be sum up as

$$\sigma_{e}=0\;\;,\;\;\omega_{e}=\omega\;\;,\;\;\mathbf{\hat{q}}_{e}^{1}=\mathbf{q}_{1}^{r}\;\;,\;\;\mathbf{\hat{q}}_{e}^{2}=\mathbf{q}_{1}^{i}$$

(20)

Before applying this extended mean flow analysis to the circular cylinder flow configuration, it is worthwhile to brieflty mention one property of the extended eigenvalue problem 19. It can be shown that if $\omega_{e}+\rm{i}\sigma_{e}$ is an eigenvalue associated to the eigenvector $\mathbf{\hat{Q}}_{e}$, then $\omega_{e}-\rm{i}\sigma_{e}$ is also an eigenvalue associated to the same eigenvector $\mathbf{\hat{Q}}_{e}$. The extended eigenvalue spectrum is therefore symmetric with respect to the axis $\sigma_{e}=0$.

The two-dimensional cylinder flow configuration is shown in 1 which indicates the lengths of the computational domain $\Omega$ in the streamwise $x$ and cross-stream $y$ directions. The nonlinear unsteady equations, nonlinear steady equations and eigenvalue problems introduced in the previous section are all solved using a finite-element method for the spatial discretization, as in [3]. To solve the nonlinear unsteady or steady equations, a unit streamwise velocity profile is imposed at the inlet boundary $\Gamma_{in}$ while the natural outflow boundary condition $-p+\textit{Re}^{-1}\partial_{x}u=0$, $\partial_{x}v=0$ is imposed at the outlet $\Gamma_{out}$, where $\mathbf{u}=(u,v)^{T}$ are the velocity components in the considered reference coordinate system. On the top and bottom boundaries, $\Gamma_{top}$ and $\Gamma_{bottom}$, symmetry conditions $\partial_{y}u=0$, $v=0$ are used. For solving the eigenvalue problems, these boundary conditions are applied in homogeneous form. The mesh is made of $47000$ triangles, with a refinement region close to the body wall and in the cylinder near-wake, resulting in a minimum mesh size of $\Delta x\approx\Delta y\approx 0.02$. Unsteady direct numerical simulations are performed using the non-dimensional time steps of $\Delta t=0.02$ for $\textit{Re}=60$, $\Delta t=0.015$ for $\textit{Re}=100$ and $\Delta t=0.01$ for $\textit{Re}=150$.

The two-dimensional base-flow computed for $\textit{Re}=100$ is depicted in 2(a) with the streamwise velocity. The stability analysis of this base flow is performed by computing the eigenvalues with largest growth rates in 5. They are depicted with open circles in 3(a). The large open circle, that lies in the gray plane, highlights the position of the unstable eigenvalue ($\sigma_{b}>0$). Initializing the unsteady nonlinear simulations with a superposition of the base flow and a small-amplitude of the corresponding eigenmode (not shown here), a linear growth followed by a nonlinear saturation of the fluctuation is observed. The flow then settles down to a fully developed vortex-shedding regime characterized by a period $T=6.01$ corresponding to the fundamental frequency $\omega=1.0446$. Once the limit cycle is established and the flow period determined, the solution is advanced in time over one period to compute the mean flow as well as the first and second harmonics using 10. The mean-flow shown in 2(b) displays a much shorter recirculation region than the base-flow. The first-harmonic, displayed in 2(c) with its real part, is normalized using 9 where the spatial domain $\Omega_{s}$ extends in the range $-15\leq x\leq 15$ and $-15\leq y\leq 15$. The amplitude of the first harmonic is $A=1.5472$ for that Reynolds number. 1 reports the amplitude and frequency of the limit cycles obtained for other values of the Reynolds number.

The first harmonic exhibits an oscillating pattern in the streamwise direction and its streamwise velocity breaks the symmetry of the mean flow with respect to the $x$-axis. The spatial pattern of the second harmonic, shown in 2(d), also oscillates in the streamwise direction, with a wavelength approximatively equal to half of the first-harmonic wavelength. The symmetry of the streamwise velocity is now similar to the one for the mean flow.

The mean-flow eigenvalues, computed by solving 16 and first reported by Barkley [4] for that flow configuration, are shown in 3(a) with black circles. The large black circle highlights the almost marginal eigenvalue. For $\textit{Re}=100$, its growth rate is $\sigma_{m}=0.0020$ and its frequency $\omega_{m}=1.0322$. Values obtained for other values of the Reynolds number are reported in 2. In all cases, the mean-flow frequency $\omega_{m}$ gives a much better approximation of the fundamental frequency $\omega$ than the base-flow frequency $\omega_{b}$.

The spatial structures of the first harmonic and marginal mean-flow mode and first harmonic are shown in 4(a) and (b), respectively. The magnitude of the velocity, equal to the square root of the mean kinetic energy, is displayed. In both cases, the largest amplitude is reached in the vicinity of the mean recirculation region. In the near-wake of the cylinder, the mean-flow eigenmode well approximates the first harmonic. In the far-wake of the cylinder, the amplitude of the mean-flow eigenmode does not decay, unlike what is observed for the first harmonic. We thus conclude that neglecting the quadratic interaction $(\mathcal{C}(\mathbf{q}_{2})\mathbf{q}_{1}^{*})$ in the first-harmonic equation does not impact neither the frequency nor the near-wake of the fluctuation, but changes the spatial decay in the far-wake.

The extended eigenvalue problem 19 is solved using the mean-flow and the second-harmonic previously shown for $\textit{Re}=100$. We are now interested by determining real eigenvalues. Indeed, complex eigenvalues of the extended problem write $\omega_{e}+\rm{i}\sigma_{e}$, where their real part corresponds to a physical frequency and their imaginary part to a physical growth rate. The extended eigenvalue spectrum is depicted in 3(b). First, we observe that the spectrum is symmetric with respect to the axis $\sigma_{e}=0$, meaning that when $\omega_{e}+\rm{i}\sigma_{e}$ is an eigenvalue, then $\omega_{e}-\rm{i}\sigma_{e}$ is also an eigenvalue. This is due to the structure of the extended operator, as noticed in the previous section. Secondly, the large white rectangle highlights the position of an eigenvalue close to the zero growth rate axis. The close-up view around that position displayed in 3(c) reveals that there are two eigenvalues lying exactly on the axis (zero growth rate) and having close but different frequencies. For comparison, the position of the marginally stable mean-flow eigenvalue is also reported in the figure with the black circle.

The exact values obtained for the growth rate and frequency of these two eigenvalues are reported in 2 for the three values of Reynolds number investigated. For the growth rates $\sigma_{e}$, we systematically report the value $0$ since they are lower than $10^{-11}$ in all cases. For the frequency, it is interesting to note that one extended eigenvalue gives a very good approximation ($\omega_{e}=1.0441$, right column) of the limit-cycle frequency ($\omega=1.0446$ for $\textit{Re}=100$), while the other one ($\omega_{e}=1.0164$, left column) gives a poorer approximation compared to the mean flow frequency ($\omega_{m}=1.0322$).

Let us now examine the spatial structure of the extended eigenmodes associated to these two eigenvalues, and called the high-frequency and low-frequency extended modes. They are respectively displayed in 5(b) and (c) by means of the magnitude of the velocity field. For comparison, the first harmonic is reproduced in 5(a). A substantial difference is observed in the shape of the two modes. Consistently with the above results, the high-frequency mode (5-b) is found to better approximate the shape of the first harmonic (5-a). In particular, the spatial decay of the velocity in the far-wake of the cylinder is well reproduced by the high-frequency mode, unlike the low-frequency mode (5-c) which even exhibits a spatial growth of the velocity in the far-wake. The same tendency has been observed for all values of the Reynolds number investigated. To select between the two extended modes without knowing the structure of the first harmonic, 13 governing the second-harmonic can be solved by replacing (in the right-hand side term) the first-harmonic with the spatial structure given by each extended mode. Results are not shown here but indicate that, consistently with the above results, the solution obtained with the high-frequency extended mode is close to the second-harmonic, unlike the one obtained with the low-frequency extended mode.

The extended mean-flow analysis has been introduced in this paper to analyze limit cycles. It is based on an eigenvalue analysis of an extended operator, that appears naturally when reformulating the equation of the first-harmonic of the time-periodic flow, and that depends not only on the mean-flow but also on the second-harmonic. Their effects are thus both taken into account, unlike the mean-flow analysis. This analysis has been applied to the circular cylinder flow configuration and revealed the existence of two eigenmodes whose growth rate is exactly equal to zero. One of this eigenmode oscillates at a frequency almost equal to the frequency of the time-periodic flow, and its spatial structure exhibits a decay of energy observed in the first-harmonic, but not in the eigenmode obtained with the mean-flow analysis.

This project has received funding from the European Research Council (ERC) under the European Union Horizon H2020 research and innovation programm (grant agreement number 638307).