# Time-accurate calculation and bifurcation analysis of the incompressible flow over a square cavity using variational multiscale modeling

Philippe Meliga, Elie Hachem

► To cite this version:

Philippe Meliga, Elie Hachem. Time-accurate calculation and bifurcation analysis of the incompressible flow over a square cavity using variational multiscale modeling. Journal of Computational Physics, 2019, 376, pp.952-972. 10.1016/j.jcp.2018.09.036. hal-01946893

HAL Id: hal-01946893

https://hal.science/hal-01946893

Submitted on 17 Feb 2020

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Time-accurate calculation and bifurcation analysis of the incompressible flow over a square cavity using variational multiscale modeling

P. Meliga<sup>a,b</sup>, E. Hachem<sup>b,*</sup>

<sup>a</sup>Aix-Marseille Univ., CNRS, Centrale Marseille, M2P2, Marseille, France

<sup>b</sup>MINES ParisTech, PSL Research University, Centre de mise en forme des matériaux (CEMEF), CNRS UMR 7635, 06904 Sophia Antipolis Cedex, France

## Abstract

A thorough variational multiscale (VMS) modeling of the Navier-Stokes equations is used to compute numerical solutions of the incompressible flow over an open cavity. This case features several competing instabilities, and is highly challenging for VMS methods with regard to frequency and pattern selection, because of the non-normality of the linearized Navier-Stokes operator. The relevance of the approach is thus carefully assessed by comparing to direct numerical simulation (DNS) data benchmarked at several Reynolds numbers, and highly accurate time advancing methods are shown to predict relevant evolutions of the transient and saturated solutions. The VMS reduces substantially the computational cost, by ~35% (resp. ~60%) in terms of CPU time using a semi-implicit discretization scheme based on backward differentiation formula (resp. the implicit Crank-Nicholson scheme), and by ~80% in terms of memory requirement. Eventually, the highly efficient semi-implicit VMS numerical framework is used to unravel the onset of the flow oscillations and the selection of the limit cycle frequency, that happens to involve a subcritical Neimark-Sacker bifurcation.

**Keywords:** Variational Multiscale modeling; Finite elements; Cavity flows; Limit cycle oscillations; Bifurcations.

<footnote>*Corresponding author
Email address: elie.hachem[@]mines-paristech.fr (E. Hachem)</footnote>
# 1. Introduction

The variational multiscale (VMS) modeling of the Navier-Stokes equations [1-5] has been widely developed for the numerical simulation of turbulent flows in several benchmarking and applicative context [6-8]. While the direct numerical simulation (DNS) requires a complete representation of the whole range of spatial and temporal turbulent scales at the discrete level (which may be untractable for several problems of practical interest, even at moderate Reynolds numbers), the VMS introduces an a priori decomposition of the solution into coarse and fine scale components, that correspond to different scales (or levels) of resolution. The general idea is that only the large scales of the flow field are fully represented and resolved at the discrete level, while the effect of the small unresolved scales is taken into account by means of consistently derived source terms proportional to the residual of the resolved scale solution, hence a more affordable computational cost.

In the context of finite elements (that remain widely used to simulate engineering relevant scenarios from computational fluid dynamics, because of the ability to handle complex geometries), the DNS solves boundary problems derived from a weak formulation of the Navier-Stokes equations (the so-called Galerkin method). The VMS add integrals over element interiors whose benefit is twofold. First, it allows relaxing the Babuska-Brezzi condition and thus accommodating equal-order interpolations for the velocity and the pressure [9, 10]. Second, it prevents the numerical instabilities developing in advection regimes at high Reynolds numbers, that induce numerical oscillations quickly spreading in the shear regions and polluting the entire solution domain. In short, the finite element VMS method provides a flexible and consistent framework for defining stable spatial approximation schemes, while also suitably representing the turbulence, which is now supported by an extensive body of work [11-15].

The question raised in this research is whether it can be similarly accurate at lower Reynolds numbers. This may seem uncanny at first glance, but such regimes can prove challenging because the linearized Navier-Stokes operator is very non-normal. In return, the stabilization terms can thus accidentally displace its eigenvalues by a substantial amount, even though their amplitude is small [16-18]. While this is not necessarily visible at high Reynolds num-
bers, where the induced variations remain small with respect to the leading growth rates, the effect can be dramatic close to the instability threshold. This is especially true under concurring linear instability mechanisms, since small eigenvalues variations lead to inaccurate instability thresholds, which may feed through the nonlinear interaction between unstable modes, and eventually affect the pattern and frequency selection. Just as with the DNS, we shall see that the key lies in highly efficient and accurate time advancing methods capable of predicting relevant evolutions of the transient and saturated solutions, a topic that has only recently become a focus of research among in the field of VMS modeling [19, 20].

The retained configuration is the open cavity flow documented in [21], that becomes unstable at Reynolds number $Re = 4140$ (based on the free-stream velocity and the cavity height), then exhibits a well-defined oscillatory behavior up to $Re \sim 7000$. This case is especially relevant for the intended purpose, insofar as there are three concurring instabilities, whose nonlinear interactions trigger a poorly understood change in the limit cycle frequency [22]. It also ensures that the stabilization is, by design, active, as the instability threshold happens to be much higher than in other classically benchmarked configurations, e.g., cylinder flows, whose critical Reynolds number $Re \sim 50$ (based on the free-stream velocity and the cylinder diameter) is smaller by two orders of magnitude [23–26], or backwards facing step flows, whose critical Reynolds number $Re \sim 800$ (based on the centerline velocity and the step height) is smaller by almost one order of magnitude [27, 28]. The objective is twofold: first, we assess the ability of the VMS to perform accurate calculations of the limit cycle oscillations while retaining the advantages of linear approximations ($P_1$ finite elements) for both the velocity and the pressure (which we recall breaks the Babuska–Brezzi condition). Then, we use efficient and accurate VMS methods to shed new light on the sequence of bifurcations triggering the onset of unsteadiness and the change in the limit cycle frequency.

The paper is organized as follows. The open cavity flow is described in Sec. 2. Exhaustive details regarding the numerical schemes used to perform the DNS and VMS computations are provided in Sec. 3. The statistics of the limit cycle oscillations are benchmarked in Sec. 4 using several DNS and VMS
<figure><img src="image_2.png" /><figcaption>(b)</figcaption></figure>

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

Figure 1: (a) Schematic of the cavity flow. (b) Structure of the computational mesh. The thick and thin solid lines have vertex densities $n_1 = 350$ and $n_2 = 300$. The dashed, dash-dotted and dash-dot-dotted lines have densities $n_3 = 200$, $n_4 = 150$ and $n_5 = 50$. The red, green and blue circles mark the position of the sensors S1, S2 and S3 used for benchmarking purposes.

discretizations. Physical interpretations for the obtained results are given in Sec. 5, together with a discussion on the numerical cost, as well as on the best compromise between numerical accuracy and time efficiency. Finally, the VMS is used in Sec. 6 to tackle the sequence of bifurcation triggering the onset of the limit cycle oscillations.

## 2. A numerical experiment: the open cavity flow

We consider the two-dimensional flow over the open square cavity sketched in Fig. 1. This case has received substantial attention in the past decade as a prototypal example of supercritical Hopf bifurcation leading to a periodic limit cycle [21, 22, 29, 30]. The general picture is as follows: a boundary layer develops upstream of the cavity from the position marked by the leftmost yellow circle in Fig. 1(a). It separates at the upstream cavity edge, and rolls up into large-scale vortices that plow into the downstream edge, after which a new boundary layer develops up to the position marked by the rightmost yellow circle. The break up of the large-scale vortices generates a pressure wave that travels upstream via the recirculating flow in the cavity, excites the shear layer at the upstream edge, and causes new perturbations to grow again into large-scale vortices via Kelvin-Helmholtz instability. For sufficiently large Reynolds numbers, this feedback loop leads to a linear global instability, whose nonlinear saturation yields periodic limit cycle oscillations.
In the sequel, we use Cartesian coordinates with $x$ horizontally along the cavity walls and $y$ vertically upwards. The governing equations to be solved numerically are the incompressible Navier-Stokes equations

$$ \nabla \cdot \mathbf{u} = 0, \quad \partial_t \mathbf{u} + \mathbf{u} \cdot \nabla \mathbf{u} + \nabla \cdot [p\mathbf{I} - \frac{2}{\mathrm{Re}}\boldsymbol{\varepsilon}(\mathbf{u})] = \mathbf{0}, \quad (1) $$

with $\mathbf{u} = (u, v)$ the velocity field of component $u$ (resp. $v$) in the $x$ (resp. $y$) direction, $p$ the pressure, $\mathbf{I}$ and $\boldsymbol{\varepsilon}(\mathbf{u})$ the identity and the rate of deformation tensors, and $\mathrm{Re}$ the Reynolds number. We use the free-stream velocity and the cavity height to make all quantities non-dimensional. The upstream and downstream edges of the cavity are thus at $(x, y) = (0, 0)$ and $(1, 0)$, and a uniform velocity $(u, v) = (1, 0)$ is prescribed at the inlet boundary, located at $x = -1.2$. A free-stress condition $-p\mathbf{n} + 2\mathrm{Re}^{-1}\boldsymbol{\varepsilon}(\mathbf{u}) \cdot \mathbf{n} = \mathbf{0}$ is applied at the outlet boundary, located at $x = 2.5$. A free-slip condition with zero tangential stress $\partial_y u = v = 0$ is prescribed on the cavity walls defined by $y = 0$ and $x \in [-1.2, -0.4] \cup [1.75, 2.5]$, as well as on the upper boundary located at $y = 0.5$. Finally, a no-slip condition $u = v = 0$ is imposed on the remainder of the cavity walls.

### 3. Numerical schemes

The above computational domain is denoted by $\Omega$, with $V \subset H_1(\Omega)^2$ and $Q = L^2(\Omega)$ proper functional spaces for the velocity and the pressure, respectively, chosen according to the above boundary conditions (we shall use the same notation for the continuous spaces and for their finite dimensional, discrete approximations, as the intending meaning is clear from the context).

#### 3.1. Finite elements discretization

The Delaunay-Voronoi algorithm is used to generate a triangulation of the computational domain. The mesh refinement is controlled by vertex densities, imposed on external and internal boundaries, as depicted in Fig. 1(b). The highest vertex density $n_1 = 350$ is used along the thick solid line (i.e., the size of each triangle side is $1/n_1 = 2.85 \times 10^{-3}$) to capture the small scales typical of the boundary and shear layers. For the range of Reynolds numbers from 4000 to 6500 considered herein, this allows distributing at least 15 triangles across
the boundary layer at the upstream cavity edge.<fnref n="1" /> The mesh is progressively stretched along the other fine lines arranged on either side of the shear region. Namely, the thin solid line has the second highest density $n_2 = 200$, while the dashed, dash-dotted and dash-dot-dotted lines have densities $n_3 = 150$, $n_4 = 100$ and $n_5 = 50$, respectively. These values are identical to those in [21, 22], and yield a mesh made of 194,447 triangles.

The relevance of the VMS for the cavity flow is established in Sec. 4 by cross-analyzing VMS and DNS data benchmarked at several Reynolds numbers. The DNS solves the weak form of (1), as obtained by integrating by parts the pressure and viscous terms, to give

$$ (\partial_t \mathbf{u} + \mathbf{u} \cdot \nabla \mathbf{u}, \mathbf{w}) + \frac{2}{\mathrm{Re}} (\varepsilon(\mathbf{u}), \varepsilon(\mathbf{w})) - (p, \nabla \cdot \mathbf{w}) + (\nabla \cdot \mathbf{u}, q) = 0, \quad (2) $$

with $\mathbf{w}$ (resp. $q$) an arbitrary test function in $V$ (resp. $Q$), and $(\cdot, \cdot)$ the $L^2$ inner product on $\Omega$. The Babuska-Brezzi condition imposes to discretize $V$ and $Q$ with different interpolation orders, which is done here with Taylor-Hood $P_2-P_1$ elements, i.e. continuous piecewise linear $P_1$ elements for $Q$ and continuous piecewise quadratic $P_2$ elements for $V$. This yields three degrees of freedom per triangle for the pressure (at the edges) and six for each velocity component (at the edges and at the middle of each side), hence 879,037 degrees of freedom for the present two-dimensional case. For proof of spatial accuracy, the reader is referred to [21], where the effect of refining the mesh up to 418,330 triangles (hence 1,888,003 degrees of freedom) on the flow linear and weakly nonlinear bifurcation parameters is reported to affect only the third decimal place.

The VMS conversely solves a stabilized formulation of (2) meant to circumvent the Babuska-Brezzi condition. We use here the simplest $P_1-P_1$ elements (i.e., $P_1$ elements for both the velocity and the pressure) yielding 294,087 degrees of freedom, hence a reduction by nearly 70% with respect to the DNS. We shall not go into the extensive details about the derivation of the stabilized formulation itself, for which the reader is referred to [31]. Suffice it to say here

<footnote n="1">¹We define the boundary layer thickness as the distance from the cavity wall to the point across the boundary layer where the $x$-velocity has reached 99% of its maximum velocity, which is because the velocity profile exhibits an overshoot in the near-wall region on behalf of a pressure gradient in the outer flow [21].</footnote>
that the flow quantities are split into coarse and fine scale components, that correspond to different scales (or levels) of resolution. The fine scales are solved in an approximate manner to allow modeling their effect into the large-scale equations, which gives rise to additional terms in the right-hand side of (2). In practice, this requires overlooking the self advection of the small scale momentum. Also, the advection time scale of the small scale is taken to be much smaller than that of the large scale, resulting in steady small scale equations using a quasi-static argument [31]. Ultimately, this yields a weak form for the large scale

$$ (\partial_t \mathbf{u} + \mathbf{u} \cdot \nabla \mathbf{u}, \mathbf{w}) + \frac{2}{\mathrm{Re}} (\varepsilon(\mathbf{u}), \varepsilon(\mathbf{w})) - (p, \nabla \cdot \mathbf{w}) + (\nabla \cdot \mathbf{u}, q) = \sum_{K \in \mathcal{T}_h} [(\tau_1 \mathcal{R}_M, \mathbf{u} \cdot \nabla \mathbf{w})_K + (\tau_1 \mathcal{R}_M, \nabla q)_K + (\tau_2 \mathcal{R}_C, \nabla \cdot \mathbf{w})_K], \quad (3) $$

where $(\cdot, \cdot)_K$ is the inner product on element $K$, $\mathcal{R}_{C,M}$ being the residuals defined by

$$ -\mathcal{R}_C = \nabla \cdot \mathbf{u}, \quad -\mathcal{R}_M = \partial_t \mathbf{u} + \mathbf{u} \cdot \nabla \mathbf{u} + \nabla p, \quad (4) $$

and $\tau_{1,2}$ are ad-hoc mesh-dependent stabilization parameters defined in [11, 15] from a Fourier analysis of the small scale problem; see also [9] for a review on the various ways to calculate these parameters in terms of the ratios of the norms of relevant matrices or vectors, as well as of local length scales and element-level Reynolds numbers.

## 3.2. Time discretization

The viscous, pressure and divergence terms are integrated implicitly with the backward Euler scheme. The discretization of the material time derivative and VMS stabilization terms is addressed in Secs. 3.2.1–3.2.2. All linear systems are preconditioned with a block Jacobi method supplemented by an incomplete LU factorization, and solved with the GMRES algorithm. The tolerance threshold is set to $10^{-6}$ in terms of both the absolute and relative norm of the residual. Decreasing the tolerance to $10^{-8}$ has been found to lead to a significant increase of the computational cost without any noticeable improvement of the numerical results, as also reported in [20].
### 3.2.1. Implicit Crank-Nicholson

We consider first the implicit Crank-Nicholson (CN) scheme

$$ \partial_t \mathbf{u} + \mathbf{u} \cdot \nabla \mathbf{u} = \frac{\mathbf{u}^i - \mathbf{u}^{i-1}}{\Delta t} + \frac{1}{2}(\mathbf{u}^i \cdot \nabla \mathbf{u}^i + \mathbf{u}^{i-1} \cdot \nabla \mathbf{u}^{i-1}), \quad (5) $$

where the i superscript refers to the solution at time $t_i = i\Delta t$. A Newton root finding algorithm is used to handle the nonlinearity of (5) in $\mathbf{u}^i$: at the j-th Newton iterate, the implicit component of the convective derivative is linearized about the last computed solution $\mathbf{u}^{i,j-1}$ (hence, $\mathbf{u}^{i,0} = \mathbf{u}^{i-1}$ and $\mathbf{u}^{i,n_i} = \mathbf{u}^i$ with $n_i$ the number of iterations needed to achieve convergence), to give

$$ \mathbf{u}^i \cdot \nabla \mathbf{u}^i = \mathbf{u}^{i,j-1} \cdot \nabla \mathbf{u}^{i,j} + \mathbf{u}^{i,j} \cdot \nabla \mathbf{u}^{i,j-1} - \mathbf{u}^{i,j-1} \cdot \nabla \mathbf{u}^{i,j-1}. \quad (6) $$

All VMS stabilization terms in (3) are integrated explicitely with the forward Euler scheme. This yields the fully discretized stabilized weak form

$$ \begin{aligned} & \left(\frac{1}{\Delta t}\mathbf{u}^{i,j} + \frac{1}{2}(\mathbf{u}^{i,j} \cdot \nabla \mathbf{u}^{i,j-1} + \mathbf{u}^{i,j-1} \cdot \nabla \mathbf{u}^{i,j}), \mathbf{w}\right) + \frac{2}{\mathrm{Re}}(\varepsilon(\mathbf{u}^{i,j}), \varepsilon(\mathbf{w})) \\ & - (\mathbf{p}^{i,j}, \nabla \cdot \mathbf{w}) + (\nabla \cdot \mathbf{u}^{i,j}, \mathbf{q}) \\ & - \left(\frac{1}{\Delta t}\mathbf{u}^{i-1} + \frac{1}{2}(\mathbf{u}^{i,j-1} \cdot \nabla \mathbf{u}^{i,j-1} - \mathbf{u}^{i-1} \cdot \nabla \mathbf{u}^{i-1}), \mathbf{w}\right) = \\ & \sum_{K \in \mathcal{T}_h} \left[ (\tau_1^{i,j-1} \mathcal{R}_M^{i,j}, \mathbf{u}^{i,j-1} \cdot \nabla \mathbf{w})_K + (\tau_1^{i,j-1} \mathcal{R}_M^{i,j}, \nabla q)_K + (\tau_2^{i,j-1} \mathcal{R}_C^{i,j}, \nabla \cdot \mathbf{w})_K \right] \\ & - \sum_{K \in \mathcal{T}_h} (\tau_1^{i,j-1} \mathcal{R}_M^{i,j-1}, \mathbf{u}^{i,j-1} \cdot \nabla \mathbf{w} + \nabla q)_K, \end{aligned} \quad (7) $$

with residuals defined as

$$ -\mathcal{R}_C^{i,j} = \nabla \cdot \mathbf{u}^{i,j}, \quad (8) $$

$$ -\mathcal{R}_M^{i,j} = \frac{\mathbf{u}^{i,j}}{\Delta t} + \frac{1}{2}(\mathbf{u}^{i,j} \cdot \nabla \mathbf{u}^{i,j-1} + \mathbf{u}^{i,j-1} \cdot \nabla \mathbf{u}^{i,j}) + \nabla p^{i,j}, \quad (9) $$

$$ -\mathcal{R}_M^{i,j-1} = \frac{\mathbf{u}^{i-1}}{\Delta t} + \frac{1}{2}(\mathbf{u}^{i,j-1} \cdot \nabla \mathbf{u}^{i,j-1} - \mathbf{u}^{i-1} \cdot \nabla \mathbf{u}^{i-1}), \quad (10) $$

and the weak form for the DNS deduces by forcing $\tau_{1,2}^{i,j-1}$ to zero.

### 3.2.2. Semi-implicit BDF

Alternatively, we consider the semi-implicit scheme

$$ \partial_t \mathbf{u} + \mathbf{u} \cdot \nabla \mathbf{u} = \frac{\alpha_{\sigma} \mathbf{u}^i - \mathbf{u}_{\mathrm{BDF}\sigma}^{i-1}}{\Delta t} + \mathbf{u}_{\mathrm{BNG}\sigma}^{i-1} \cdot \nabla \mathbf{u}^i, \quad (11) $$

where the time and convective derivatives are approximated with backward differentiation formulas (BDF) and backward Newton-Gregory polynomials (BNG)
<table><caption>Table 1: Properties of the optimized, multi-time-level BDF3† and BDF4† schemes.</caption><thead><tr><th>σ<sup>†</sup></th><th>χ<sub>2</sub></th><th>χ<sub>3</sub></th><th>χ<sub>4</sub></th></tr></thead><tbody><tr><td>3</td><td>1</td><td>1</td><td>0</td></tr><tr><td>4</td><td>5 - 3√2</td><td>-5 + 4√2</td><td>2 - √2</td></tr></tbody></table>

of same order $\sigma$, respectively. This yields

$$ \alpha_{\sigma} = \sum_{k=1}^{\sigma} \frac{1}{k}, \quad \mathbf{u}_{\text{BDF}\sigma}^{i-1} = \sum_{k=1}^{\sigma} \beta_k \mathbf{u}^{i-k}, \quad \mathbf{u}_{\text{BNG}\sigma}^{i-1} = \sum_{k=1}^{\sigma} \gamma_k \mathbf{u}^{i-k}, \quad (12) $$

where $\{\beta_k, \gamma_k\}$ are sequences defined recursively after [32, 33] by $\beta_1 = \gamma_1 = \sigma$ and

$$ \beta_{k \ge 2} = \beta_{k-1} \frac{k(k-\sigma)}{k^2}, \quad \gamma_{k \ge 2} = \frac{k-\sigma}{k+1} \gamma_{k-1}, \quad (13) $$

hence $\beta_{k>\sigma}$ and $\gamma_{k>\sigma}$ are trivially zero. Again, the VMS stabilization terms are integrated explicitly with the forward Euler scheme. This gives rise to the fully discretized stabilized weak form

$$ \begin{aligned} & \left(\frac{\alpha_{\sigma}}{\Delta t} \mathbf{u}^i + \mathbf{u}_{\text{BNG}\sigma}^{i-1} \cdot \nabla \mathbf{u}^i, \mathbf{w}\right) + \frac{2}{\text{Re}}(\varepsilon(\mathbf{u}^i), \varepsilon(\mathbf{w})) \\ & - (p^i, \nabla \cdot \mathbf{w}) + (\nabla \cdot \mathbf{u}^i, q) \\ & - \left(\frac{1}{\Delta t} \mathbf{u}_{\text{BDF}\sigma}^{i-1}, \mathbf{w}\right) = \\ & \sum_{K \in \mathcal{T}_h} \left[ \left(\tau_{1\sigma}^{i-1} \mathcal{R}_M^i, \mathbf{u}_{\text{BNG}\sigma}^{i-1} \cdot \nabla \mathbf{w}\right)_K + (\tau_{1\sigma}^{i-1} \mathcal{R}_M^i, \nabla q)_K + (\tau_{2\sigma}^{i-1} \mathcal{R}_C^i, \nabla \cdot \mathbf{w})_K \right] \\ & - \sum_{K \in \mathcal{T}_h} (\tau_{1\sigma}^{i-1} \mathcal{R}_M^i, \mathbf{u}_{\text{BNG}\sigma}^{i-1} \cdot \nabla \mathbf{w} + \nabla q)_K, \end{aligned} \quad (14) $$

with residuals defined as

$$ -\mathcal{R}_C^i = \nabla \cdot \mathbf{u}^i, \quad (15) $$

$$ -\mathcal{R}_M^i = \frac{\alpha_{\sigma}}{\Delta t} \mathbf{u}^i + \mathbf{u}_{\text{BNG}\sigma}^{i-1} \cdot \nabla \mathbf{u}^i + \nabla p^i, \quad (16) $$

$$ -\mathcal{R}_M^{i-1} = \frac{1}{\Delta t} \mathbf{u}_{\text{BDF}\sigma}^{i-1}, \quad (17) $$

and the weak form for the DNS deduces by forcing $\tau_{1,2}^{i,j-1}$ to zero.

It is worth mentioning that only BDF1 and BDF2 are unconditionally stable, a result known as the second Dahlquist barrier. Without anticipating the results of the upcoming sections, this is an issue because both schemes perform quite
poorly on the cavity flow, and we did run into numerical instabilities with BDF3 and BDF4; see also [34]. We thus circumvent the difficulty using modified schemes combining linearly BDF2, BDF3 and BDF4, i.e.,

$$ \alpha_{\sigma\dagger} = \frac{1}{2} \sum_{\sigma'=2}^{\sigma} \chi_{\sigma'} \alpha_{\sigma'}, \quad \mathbf{u}_{\text{BDF}\sigma\dagger}^{i} = \frac{1}{2} \sum_{\sigma'=2}^{\sigma} \chi_{\sigma'} \mathbf{u}_{\text{BDF}\sigma'}^{i}, \qquad (18) $$

where the $\chi_{\sigma'}$ coefficients reported in Table 1 for $\sigma^\dagger = 3$ and 4 are computed analytically in [35] to minimize the leading-order truncation error at each time step. These so-called optimized schemes, denoted by BDF3$^\dagger$ and BDF4$^\dagger$ in the following, are shown in the aforementioned reference to formally maintain the second-order accuracy of BDF2 but decrease the truncation error by a factor $2/\chi_2$ (hence, 2 for BDF3$^\dagger$ and $\sim 2.65$ for BDF4$^\dagger$) while maintaining the numerical stability.

## 4. Benchmark results

This section reports VMS and DNS data pertaining to several Reynolds numbers in the periodic regime. For a given numerical framework, it assesses the accuracy from the statistics of the limit cycle oscillations measured by a sensor labeled S2, whose position $\mathbf{x}_{S2} = (0.75, 0.05)$ in the shear layer is marked by the green circle in Fig. 1. For benchmarking purposes, we also report the values measured by two additional sensors S1 and S3, whose positions $\mathbf{x}_{S1} = (0.6, 0.05)$ and $\mathbf{x}_{S3} = (0.9, 0.05)$ are marked by the red and blue circles, respectively. The time step $\Delta t = 10^{-2}$ is the same as in [22], and allows sampling between $\sim 55$ and 85 points over a period, depending on the Reynolds number. We consider a time span of 25 oscillation periods, which allows determining the mean velocity $\bar{u}$, the root mean square (rms) of the velocity fluctuations $u'$ and the fundamental Strouhal number $\text{St}$ (checked to be independent of the sensor position) to a numerical accuracy of 0.5%, as has been assessed comparing to results obtained from a doubled time span.

### 4.1. Case I - Re = 5000

First, we examine the results obtained at Re = 5000 with the Crank-Nicholson scheme described in Sec. 3.2.1, for VMS and DNS computations starting from the exact same (albeit generated randomly) initial condition. As
<table><caption>Table 2: Statistics of the limit-cycle oscillation at Re = 5000.</caption><thead><tr><th rowspan="2">Discretization</th><th rowspan="2">Δt</th><th colspan="3">ū (×10<sup>-1</sup>)</th><th colspan="3">u' (×10<sup>-2</sup>)</th><th rowspan="2">St</th></tr><tr><th>S1</th><th>S2</th><th>S3</th><th>S1</th><th>S2</th><th>S3</th></tr></thead><tbody><tr><td rowspan="3">CN</td><td>DNS</td><td>10<sup>-2</sup></td><td>9.965</td><td>9.730</td><td>9.315</td><td>2.380</td><td>4.065</td><td>4.443</td><td>1.692</td></tr><tr><td>VMS</td><td>10<sup>-2</sup></td><td>9.959</td><td>9.724</td><td>9.312</td><td>2.350</td><td>4.018</td><td>4.404</td><td>1.692</td></tr><tr><td>VMS</td><td>5×10<sup>-3</sup></td><td>9.960</td><td>9.726</td><td>9.316</td><td>2.345</td><td>4.007</td><td>4.394</td><td>1.694</td></tr><tr><td rowspan="3">BDF1</td><td>DNS</td><td>10<sup>-2</sup></td><td>9.984</td><td>9.848</td><td>9.634</td><td>0.375</td><td>0.726</td><td>1.016</td><td>1.177</td></tr><tr><td>VMS</td><td>10<sup>-2</sup></td><td>9.978</td><td>9.841</td><td>9.627</td><td>0.403</td><td>0.778</td><td>1.084</td><td>1.178</td></tr><tr><td>VMS</td><td>5×10<sup>-3</sup></td><td>9.978</td><td>9.821</td><td>9.560</td><td>0.932</td><td>1.812</td><td>2.509</td><td>1.204</td></tr><tr><td rowspan="3">BDF2</td><td>DNS</td><td>10<sup>-2</sup></td><td>9.955</td><td>9.683</td><td>9.225</td><td>2.751</td><td>4.550</td><td>4.659</td><td>1.703</td></tr><tr><td>VMS</td><td>10<sup>-2</sup></td><td>9.948</td><td>9.677</td><td>9.224</td><td>2.728</td><td>4.516</td><td>4.628</td><td>1.703</td></tr><tr><td>VMS</td><td>5×10<sup>-3</sup></td><td>9.958</td><td>9.716</td><td>9.297</td><td>2.428</td><td>4.117</td><td>4.450</td><td>1.698</td></tr><tr><td rowspan="3">BDF3<sup>†</sup></td><td>DNS</td><td>10<sup>-2</sup></td><td>9.967</td><td>9.736</td><td>9.329</td><td>2.322</td><td>3.984</td><td>4.401</td><td>1.694</td></tr><tr><td>VMS</td><td>10<sup>-2</sup></td><td>9.960</td><td>9.730</td><td>9.327</td><td>2.299</td><td>3.947</td><td>4.362</td><td>1.694</td></tr><tr><td>VMS</td><td>5×10<sup>-3</sup></td><td>9.961</td><td>9.729</td><td>9.322</td><td>2.325</td><td>3.979</td><td>4.377</td><td>1.695</td></tr><tr><td rowspan="3">BDF4<sup>†</sup></td><td>DNS</td><td>10<sup>-2</sup></td><td>9.967</td><td>9.735</td><td>9.326</td><td>2.333</td><td>3.999</td><td>4.404</td><td>1.694</td></tr><tr><td>VMS</td><td>10<sup>-2</sup></td><td>9.960</td><td>9.729</td><td>9.325</td><td>2.309</td><td>3.961</td><td>4.370</td><td>1.694</td></tr><tr><td>VMS</td><td>5×10<sup>-3</sup></td><td>9.961</td><td>9.728</td><td>9.322</td><td>2.328</td><td>3.981</td><td>4.377</td><td>1.695</td></tr></tbody></table>

evidenced in Fig. 2(a) showing the time evolution of $u_{S2}$ (i.e., the x-component of the velocity measured by the S2 sensor, here computed by VMS), well-defined oscillations emerge after roughly 25 time units, and the solution subsequently slowly converges to its limit cycle within 50 time units. The VMS predicts especially well the transient evolution in Fig. 2(b) (CN data are shown as the grey lines, with solid and dashed patterns for VMS and DNS results). There are limited discrepancies at each time step, that trigger a small phase shift in the limit cycle solution<fnref n="2" />, but the latter is erased in Fig. 2(c) to ease the presentation. Other than that, the limit cycle statistics in Table 2 clearly assess the accuracy

<footnote n="2">²The latter is arbitrary in the sense that a slight phase advance or delay of the order of one tenth of the oscillation period is observed depending on the scheme and the Reynolds number.</footnote>
<figure><img src="image_4.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(c)</figcaption></figure>

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

Figure 2: (a) Time evolution of the x-velocity at the S2 sensor, computed at $Re = 5000$ by VMS with the Crank-Nicholson scheme. (b,c) Focus on the (b) transient and (c) periodic regimes, as computed by VMS (solid lines) and DNS (dashed lines) with the Crank-Nicholson (grey lines) and BDF4† (black lines) scheme. The considered time intervals are marked in (a) as the grey boxes.

of the VMS solution. The oscillation amplitude is especially well predicted, the discrepancy by 0.2% being of the same order as the statistical error.

We now take the above CN data as reference, and compare against the values obtained with the BDF schemes described in Sec. 3.2.2. We compare VMS-BDF against VMS-CN and DNS-BDF against DNS-CN, considering that all computations start from the same initial condition. We quickly pass over BDF1, found in Table 2 to yield irrelevant limit cycle oscillations, whether it be in terms of frequency or amplitude. This is true regardless of the numerical framework, and is ascribed to an excessive dissipativity, as will be further discussed in Sec. 5. BDF2 restores limit cycle oscillations of the proper frequency but performs poorly overall, as the oscillation amplitude is widely overestimated, by 12.4% using VMS and 13.3% using DNS. The data in Table 2 indicate that the optimized BDF schemes improve significantly the accuracy, namely BDF3†
<figure><img src="image_7.png" /></figure>

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

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

<figure><img src="image_10.png" /><figcaption>Figure 3: (a,b) Isocontours of the fluctuating x-velocity, computed at $Re = 5000$ by (a) VMS and (b) DNS with the Crank-Nicholson scheme. (c,d) Same as (a,b) using the BDF4† scheme. All results pertain to the same phase in the limit cycle, corresponding to a peak velocity at the S2 sensor.</figcaption></figure>

cuts the VMS error to 1.7% and BDF4† cuts it further down to 1.3% (which incidentally assesses the convergence in terms of the approximation order). The VMS therefore has a residual (albeit small) error associated with it, whereas the DNS results obtained with the same schemes are accurate to 0.5%. Again, this is a matter to which we return later. If we now pay a closer attention to BDF4†, a relevant evolution of the transient sensor velocity is computed by VMS; see the solid and dashed lines in Fig. 2(b). Note, there are subtle discrepancies with respect to the CN results (in grey), owing to the fact that we initially time-march the system with BDF1 and progressively upgrade to BDF4† during the first three time units.

The high accuracy of the VMS-CN and VMS-BDF4† numerical frameworks is further illustrated in Fig. 3 showing the spatial distribution of the x-velocity at the same phase in the limit cycle, corresponding to a peak velocity at the S2 sensor, together with the mean shear layer position (as marked by the superimposed grey line showing the time-averaged streamline issuing from the upstream cavity edge). The roll-up mechanism yields counter-rotating vortical structures
<figure><img src="image_11.png" /><figcaption>Figure 4: Cuts of the *x*-velocity at (a) *y* = 0.05, (b) *x* = 0.6, (c) *x* = 0.75 and (d) *x* = 0.9, computed at *Re* = 5000 by VMS (solid lines) and DNS (dashed lines) using the Crank-Nicholson (grey lines) and BDF4† (black lines) schemes, respectively. The sensor positions considered in Table 2 are marked by the coloured circles. All results pertain to the same phase in the limit cycle, corresponding to a peak velocity at the S2 sensor.</figcaption></figure>

traveling and strengthening across the shear layer. The VMS captures especially well the number of vortices, the peak vorticity at the downstream cavity edge, but also the finest structures developing in the inner cavity and the downstream boundary layer. A more quantitative indication of the agreement is displayed in Fig. 4(a) showing the related horizontal cuts at *y* = 0.05, where the footprint of the vortices is visible through the succession of peaks and troughs resembling a wave distribution (again, CN and BDF4† data are shown as the grey and black lines, with solid and dashed patterns for VMS and DNS results). The same level of agreement is reported in Figs. 4(c) and 4(d) showing the evolution of the shear layer velocity through vertical cuts at *x* = 0.6, 0.75 and 0.9. The latter case unveils a severe distortion corresponding to the trail left by the vortices traveling upstream via the recirculating flow.
<table><caption>Table 3: Statistics of the limit-cycle oscillation at Re = 6000.</caption><thead><tr><th rowspan="2">Discretization</th><th rowspan="2">Δt</th><th colspan="3">ū (×10<sup>-1</sup>)</th><th colspan="3">u' (×10<sup>-2</sup>)</th><th rowspan="2">St</th></tr><tr><th>S1</th><th>S2</th><th>S3</th><th>S1</th><th>S2</th><th>S3</th></tr></thead><tbody><tr><td rowspan="3">CN</td><td>DNS</td><td>10<sup>-2</sup></td><td>10.04</td><td>9.600</td><td>9.012</td><td>4.184</td><td>6.314</td><td>5.530</td><td>1.781</td></tr><tr><td>VMS</td><td>10<sup>-2</sup></td><td>10.04</td><td>9.595</td><td>9.011</td><td>4.162</td><td>6.296</td><td>5.517</td><td>1.781</td></tr><tr><td>VMS</td><td>5×10<sup>-3</sup></td><td>10.04</td><td>9.601</td><td>9.013</td><td>4.150</td><td>6.247</td><td>5.504</td><td>1.782</td></tr><tr><td rowspan="3">BDF1</td><td>DNS</td><td>10<sup>-2</sup></td><td>10.12</td><td>9.934</td><td>9.593</td><td>2.047</td><td>3.494</td><td>4.049</td><td>1.655</td></tr><tr><td>VMS</td><td>10<sup>-2</sup></td><td>10.11</td><td>9.928</td><td>9.591</td><td>2.051</td><td>3.499</td><td>4.042</td><td>1.647</td></tr><tr><td>VMS</td><td>5×10<sup>-3</sup></td><td>10.20</td><td>9.811</td><td>9.309</td><td>3.043</td><td>4.818</td><td>4.921</td><td>1.712</td></tr><tr><td rowspan="3">BDF2</td><td>DNS</td><td>10<sup>-2</sup></td><td>10.00</td><td>9.503</td><td>8.948</td><td>4.764</td><td>6.919</td><td>5.313</td><td>1.792</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>2.095</td></tr><tr><td>VMS</td><td>10<sup>-2</sup></td><td>10.00</td><td>9.502</td><td>8.947</td><td>4.705</td><td>6.869</td><td>5.312</td><td>1.791</td></tr><tr><td rowspan="3">BDF3†</td><td>VMS</td><td>10<sup>-2</sup></td><td>10.03</td><td>9.580</td><td>8.997</td><td>4.285</td><td>6.384</td><td>5.459</td><td>1.786</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>VMS</td><td>5×10<sup>-3</sup></td><td>10.05</td><td>9.615</td><td>9.026</td><td>4.123</td><td>6.214</td><td>5.515</td><td>1.781</td></tr><tr><td rowspan="3">BDF4†</td><td>VMS</td><td>10<sup>-2</sup></td><td>10.04</td><td>9.610</td><td>9.024</td><td>4.116</td><td>6.212</td><td>5.498</td><td>1.781</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>VMS</td><td>5×10<sup>-3</sup></td><td>10.04</td><td>9.605</td><td>9.018</td><td>4.154</td><td>6.242</td><td>5.499</td><td>1.783</td></tr><tr><td rowspan="3">BDF4‡</td><td>VMS</td><td>10<sup>-2</sup></td><td>10.05</td><td>9.612</td><td>9.020</td><td>4.141</td><td>6.250</td><td>5.534</td><td>1.781</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>VMS</td><td>5×10<sup>-3</sup></td><td>10.04</td><td>9.606</td><td>9.019</td><td>4.133</td><td>6.247</td><td>5.516</td><td>1.781</td></tr><tr><td rowspan="3">BDF4‡</td><td>VMS</td><td>10<sup>-2</sup></td><td>10.04</td><td>9.604</td><td>9.017</td><td>4.156</td><td>6.244</td><td>5.500</td><td>1.783</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>VMS</td><td>5×10<sup>-3</sup></td><td>10.04</td><td>9.604</td><td>9.017</td><td>4.156</td><td>6.244</td><td>5.500</td><td>1.783</td></tr></tbody></table>

## 4.2. Case II - Re = 6000

We examine now the results obtained at Re = 6000 with the Crank-Nicholson scheme, starting from the exact same initial condition as at Re = 5000. As evidenced in Fig. 5(a) showing the time evolution of the $u_{S2}$ velocity (here computed by VMS), well-defined oscillations emerge after only 10 time units, but it takes about 350 times units for the solution to converge to its limit cycle, as the transient exhibits persistent modulations resulting from a non-normal competition between instabilities (this point will be addressed later). Also, the periodic velocity in Fig. 5(c) significantly departs from a pure wave because the nonlinearity now produces a substantial second harmonic responsible for a complex frequency selection mechanism [22]. As already noted at Re = 5000, the VMS
<figure><img src="image_12.png" /><figcaption>Figure 5: (a) Time evolution of the x-velocity at the S2 sensor, computed at $Re = 6000$ by VMS with the Crank–Nicholson scheme. (b,c) Focus on the (b) transient and (c) periodic regimes, as computed by VMS (solid lines) and DNS (dashed lines) with the Crank–Nicholson (grey lines) and BDF4† (black lines) scheme. The considered time intervals are marked in (a) as the grey boxes.</figcaption></figure>

predicts well the transient and periodic evolutions of the sensor velocity, up to a small phase factor deliberately erased in Fig. 5(c), and the VMS oscillation amplitude is accurate to 0.3%, which is of the same order as the statistical error.

If we now take the CN data as reference and compare against the BDF results, we observe the same trends as for $Re = 5000$, namely BDF1 keeps yielding off-topic limit cycle oscillations (although this is to a far lesser extent than at $Re = 5000$, for instance the oscillation amplitude is widely underestimated by 25% but the frequency is somehow accurate to 7.5%), and BDF2 keeps overestimating the amplitude (9% by VMS, vs. 9.6% by DNS). The accuracy again improves dramatically with the optimized BDF† schemes, namely BDF3† cuts the error to 1.5%, and BDF4† further cuts it down to 0.7% (as a comparison, the DNS is perfectly accurate to the statistical error). Note, for this case, BDF2 has another drawback, that has to do with the very nature of the computed solution.
<figure><img src="image_13.png" /></figure>

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

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

<figure><img src="image_16.png" /><figcaption>Figure 6: (a,b) Isocontours of the fluctuating $x$-velocity, computed at $Re = 6000$ by (a) VMS and (b) DNS with the Crank-Nicholson scheme. (c,d) Same as (a,b) using the BDF4$^{†}$ scheme. All results pertain to the same phase in the limit cycle, corresponding to a peak velocity at the S2 sensor.</figcaption></figure>

On the one hand, the VMS-BDF2 solution is periodic, as evidenced by the velocity fluctuations in Fig. 8(a), whose single-sided amplitude spectrum is shown in Fig. 8(c). The DNS however settles on a quasi-periodic orbit because the transient modulations never die out, hence the two frequencies featured in the amplitude spectrum (Fig. 8(d)) and reported in Table 3. While the dominant peak is identical to its VMS-BDF2 counterpart, there is a small, secondary peak at a slightly larger frequency that yields the parasitic, low-frequency modulation of the sensor velocity in Fig. 8(b). This inconsistency (to which we come back in Sec. 5) stresses again the need for highly accurate time-integration schemes, namely CN and BDF4$^{†}$, whose accuracy is further illustrated in Fig. 6 depicting isocontours of the fluctuating velocity. The same level of agreement is reported in Fig. 7 showing several cuts of the shear layer velocity; see especially the vertical cuts at $x = 0.9$ in Fig. 7(d) for which the shear layer velocity is strongly distorted by the inner vortices.
<figure><img src="image_17.png" /></figure>

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

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

<figure><img src="image_20.png" /><figcaption>Figure 7: Cuts of the x-velocity at (a) $y = 0.05$, (b) $x = 0.6$, (c) $x = 0.75$ and (d) $x = 0.9$, computed at $Re = 6000$ by VMS (solid lines) and DNS (dashed lines) using the Crank-Nicholson (grey lines) and BDF4† (black lines) schemes, respectively. The sensor positions considered in Table 3 are marked by the coloured circles. All results pertain to the same phase in the limit cycle, corresponding to a peak velocity at the S2 sensor.</figcaption></figure>

### 4.3. Case III - Re = 4350

Finally, we consider the results obtained at $Re = 4350$ with the Crank-Nicholson scheme, starting again from the exact same initial condition as at $Re = 5000$. As shown in Figure 9(a) showing the time evolution of the $u_{S2}$ velocity (here computed by VMS), well-defined oscillations emerge after 50 time units (which is because the Reynolds number is close to the instability threshold, and the linear growth rate of the disturbances is accordingly smaller), and the solution settles on its limit cycle after 100 time units. As already noted at $Re = 5000$ and $6000$, the VMS predicts well the transient and periodic evolutions of the sensor velocity, up to a small phase factor again deliberately erased in Fig. 9 (c). Other than that, the limit cycle statistics in Table 4 assess the accuracy of the VMS solution, whose only noticeable shortcoming is an overestimation of the
<figure><img src="image_21.png" /><figcaption>Figure 8: (a,b) Time evolution of the *x*-velocity at the S2 sensor, computed at $Re = 6000$ by (a) VMS and (b) DNS with the BDF2 scheme. (c,d) Single-sided amplitude spectrum of the corresponding (c) VMS and (d) DNS *x*-velocity fluctuations.</figcaption></figure>

rms by 1.9%. Regarding the BDF results, we shall again overlook BDF1, whose excessive dissipativity now yields an erroneous steady solution. BDF2 restores limit cycle oscillations of the proper frequency but keeps performing poorly overall, as the oscillation amplitude remains widely overestimated (13.8% by VMS, vs. 15.2% by DNS). The accuracy again improves significantly using the optimized BDF† schemes, namely BDF3† cuts the error to 4.1%, and BDF4† further cuts it down to 3.5% (which is comparable to the accuracy of the DNS, found to be to 4.3%). Although the residual error is substantially larger than at $Re = 5000$ and 6000, the VMS-CN and VMS-BDF4† numerical frameworks still predict the unsteady cavity flow to a high accuracy, as an especially close agreement is emphasized Fig. 10 displaying the isocontours of the fluctuating velocity. Again, the VMS captures well the number of vortices, that happens to be substantially smaller than at $Re = 5000$ and 6000, hence the reduced oscillation frequency. A more quantitative indication of the agreement is also
<table><caption>Table 4: Same as Table 2 at Re = 4350.</caption><thead><tr><th rowspan="2">Discretization</th><th rowspan="2">Δt</th><th colspan="3">ū (×10<sup>-1</sup>)</th><th colspan="3">u' (×10<sup>-2</sup>)</th><th rowspan="2">St</th></tr><tr><th>S1</th><th>S2</th><th>S3</th><th>S1</th><th>S2</th><th>S3</th></tr></thead><tbody><tr><td rowspan="3">CN</td><td>DNS</td><td>10<sup>-2</sup></td><td>9.844</td><td>9.686</td><td>9.433</td><td>0.773</td><td>1.503</td><td>1.998</td><td>1.202</td></tr><tr><td>VMS</td><td>10<sup>-2</sup></td><td>9.836</td><td>9.677</td><td>9.424</td><td>0.789</td><td>1.532</td><td>2.032</td><td>1.202</td></tr><tr><td>VMS</td><td>5×10<sup>-3</sup></td><td>9.837</td><td>9.678</td><td>9.426</td><td>0.780</td><td>1.512</td><td>2.007</td><td>1.203</td></tr><tr><td rowspan="3">BDF1</td><td>DNS</td><td>10<sup>-2</sup></td><td>9.847</td><td>9.703</td><td>9.484</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>VMS</td><td>10<sup>-2</sup></td><td>9.839</td><td>9.695</td><td>9.477</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>VMS</td><td>5×10<sup>-3</sup></td><td>9.840</td><td>9.696</td><td>9.478</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td rowspan="3">BDF2</td><td>DNS</td><td>10<sup>-2</sup></td><td>9.843</td><td>9.680</td><td>9.416</td><td>0.890</td><td>1.732</td><td>2.301</td><td>1.204</td></tr><tr><td>VMS</td><td>10<sup>-2</sup></td><td>9.835</td><td>9.671</td><td>9.406</td><td>0.915</td><td>1.777</td><td>2.355</td><td>1.205</td></tr><tr><td>VMS</td><td>5×10<sup>-3</sup></td><td>9.837</td><td>9.677</td><td>9.423</td><td>0.807</td><td>1.563</td><td>2.075</td><td>1.204</td></tr><tr><td rowspan="3">BDF3<sup>†</sup></td><td>DNS</td><td>10<sup>-2</sup></td><td>9.844</td><td>9.688</td><td>9.438</td><td>0.735</td><td>1.428</td><td>1.904</td><td>1.202</td></tr><tr><td>VMS</td><td>10<sup>-2</sup></td><td>9.836</td><td>9.678</td><td>9.428</td><td>0.757</td><td>1.469</td><td>1.953</td><td>1.203</td></tr><tr><td>VMS</td><td>5×10<sup>-3</sup></td><td>9.837</td><td>9.679</td><td>9.428</td><td>0.766</td><td>1.485</td><td>1.972</td><td>1.203</td></tr><tr><td rowspan="3">BDF4<sup>†</sup></td><td>DNS</td><td>10<sup>-2</sup></td><td>9.844</td><td>9.687</td><td>9.437</td><td>0.739</td><td>1.437</td><td>1.914</td><td>1.202</td></tr><tr><td>VMS</td><td>10<sup>-2</sup></td><td>9.836</td><td>9.678</td><td>9.428</td><td>0.762</td><td>1.478</td><td>1.964</td><td>1.203</td></tr><tr><td>VMS</td><td>5×10<sup>-3</sup></td><td>9.837</td><td>9.679</td><td>9.428</td><td>0.766</td><td>1.486</td><td>1.973</td><td>1.203</td></tr></tbody></table>

unveiled in the various cuts shown in Fig. 11.

The larger residual error is ascribed to the fact that the Reynolds number is only closely above the instability threshold [21, 30]. For this case of supercritical Hopf bifurcation [21], the limit cycle oscillation amplitude varies as the square-root of the departure from criticality, i.e.,

$$u' = \xi \sqrt{\frac{1}{\text{Re}_c} - \frac{1}{\text{Re}}}, \qquad (19)$$

<footnote n="290">290 with Re_c the critical Reynolds number and ξ an amplitude coefficient accounting for the saturation of the instability. It is thus very nonlinear in Re, which in turn magnifies the small inaccuracies inherent to the VMS (and also those inherent to</footnote>
<figure><img src="image_22.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure 9: Same as Figs. 2 and 5 at Re = 4350.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>Figure 10: Same as Figs. 3 and 6 at Re = 4350.</figcaption></figure>
<figure><img src="image_29.png" /></figure>

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

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

<figure><img src="image_32.png" /><figcaption>Figure 11: Same as Figs. 4 and 7 at Re = 4350.</figcaption></figure>

the time integration scheme, as evidenced by the increased variance between the DNS-CN and DNS-BDF4† data). For all that, the effect on the near-critical unsteady dynamics remains remarkably small, as has been assessed carefully from the energy of the *x*-velocity fluctuations at the sensor position, i.e., $u_{S2}^2$. As illustrated in Fig. 12 showing the results obtained for several values of Re in the range [4100; 4200], the latter exhibits the expected linear increase with Re. This allows extrapolating the instability threshold $Re_c$ from the intersection of the linear regression (performed here on the 6 most near-critical values) with the horizontal axis, and the amplitude coefficient $\xi$ from the square root of the slope (after the results have been recast in terms of the departure from threshold $Re_c^{-1} - Re^{-1}$). For the record, the theoretical values shown as the grey lines/symbols in Fig. 12 are $Re_c = 4130^3$ and $\xi = 4.18$. Without going

295

300

<footnote n="3">³Note, there is a slight difference with the value $Re_c = 4140$ reported in [21], which is because the aforementioned reference assumes the instability threshold to be determined to a sufficient accuracy when the growth rate of the most unstable disturbances is smaller than</footnote>
<figure><img src="image_33.png" /><figcaption>Figure 12: Energy of the $x$-velocity fluctuations at the S2 sensor, computed against Re with (a) VMS-CN, and (b) VMS-BDF4†. The circles are the actual data points, and the solid line is the linear regression of the 6 most near-critical values. The black vertical dots mark the so-determined instability threshold. The solid grey line is the energy computed from (19), using the theoretical instability threshold (grey vertical dots) and amplitude coefficient stemming from a multiple time scale analysis, using the same solver as in [36]. The dashed grey line is the energy corrected for the discrepancy in the instability thresholds.</figcaption></figure>

into the technical details, those are free from time-discretization errors, $P_2-P_1$ values computed by multiple time-scale analysis for the exact same mesh as in the present study, using the same solver as in [36]. By way of comparison, VMS-CN yields an instability threshold $Re_c = 4108$ accurate to 0.5% (which is smaller than the error on the amplitude by one order of magnitude) and an amplitude coefficient $\xi = 4.24$ overestimating the theoretical value by as few as 1.43%. VMS-BDF4† conversely yields a threshold $Re_c = 4128$ perfectly accurate to the statistical error, and an amplitude coefficient $\xi = 4.27$ overestimated by 2.15%. Despite the overall high accuracy, this suggests that, BDF4† predicts more accurately the linear features of the instability, while CN captures more accurately the nonlinear saturation mechanism.

# 5. Discussion

## 5.1. *Physical interpretations*

The base cavity flow (i.e., the solution to the steady Navier-Stokes equation) undergoes two Hopf bifurcations in a row, at Reynolds numbers $Re_{cA} = 4130$

<footnote>$10^{-3}$. The present study conversely uses a value $10^{-6}$.</footnote>
and $Re_{CB} = 4349$. Interestingly, the first eigenmode to bifurcate, termed A, remains more linearly unstable up to $Re_{CAB} = 4553$, after which mode B takes over. Those thresholds (just as all eigenfrequencies $\omega/2\pi$ to follow) are again free from time-discretization errors, $P_2-P_1$ values computed by linear stability analysis for the exact same mesh as in the present study, using the same solver as in [36]. The above benchmark results can thus be interpreted as follows:

* At $Re = 4350$, both modes are unstable but mode B is almost exactly neutrally stable. The limit cycle oscillations therefore proceed from the nonlinear saturation of mode A (hence, the LCA moniker), whose eigen-frequency $\omega_A/2\pi = 1.20$ is indistinguishable from the limit cycle frequency $St_A$. This is because the growth rate of mode A is small enough for the nonlinear dynamics to follow the linear theory. Somehow, it is too small to sustain the high dissipativity of BDF1, whose effect is to artificially reduce the Reynolds number (which amounts to assuming that all unstable modes are somehow damped to the same extent), hence the erroneous steady solution computed with this scheme.

* At $Re = 5000$, the limit cycle oscillations now proceed from the nonlinear saturation of mode B (LCB). The latter has taken over as the dominant instability mode, and its eigenfrequency $\omega_B/2\pi = 1.67$ is almost identical to the limit cycle frequency $St_B = 1.69$. The BDF1 solution is now periodic but erroneously settles on LCA (whose frequency $St_A = 1.18$ matches well the eigenfrequency $\omega_A/2\pi = 1.21$ of mode A computed at this Reynolds number) because the excessive dissipation requires a larger Reynolds number for LCB to take over. This illustrates how a lack of accuracy in the numerical integration can alter the way competing linear instability mechanisms interact with one another, which in turn alters the outcome of the limit cycle selection.

* At $Re = 6000$, the limit cycle oscillations continue to proceed from the nonlinear saturation of mode B (LCB). The latter remains the dominant instability mode, and its eigenfrequency $\omega_B/2\pi = 1.70$ agrees well with the limit cycle frequency $St_B = 1.78$. Note, the statistics of the BDF1 limit cycle are close to those computed at $Re = 5000$ with CN/BDF4$^\dagger$, which
supports the idea that the excessive dissipation of this scheme acts by artificially reducing the Reynolds number (only the latter is now large enough for the solution to settle on LCB). As for the quasi-periodic oscillations computed with DNS-BDF2, they are related to the occurrence of a third Hopf bifurcation at $Re_{cc} = 5685$, as the eigenfrequency $\omega_c/2\pi = 2.15$ computed at this Reynolds number is very close to the secondary fundamental $St_c = 2.10$. Even though a rigorous proof is lacking, this likely yields a Neimark-Sacker bifurcation, i.e. a Hopf bifurcation for the limit cycles, and thus a quasi-periodic T2 torus on which modes B and C oscillate with finite amplitudes. The mechanism itself is relevant, in the sense that the same quasi-periodic behavior has been observed at a larger Reynolds number $Re = 7500$ [37]. Still, the presence of such oscillations at a Reynolds number as low as 6000 is spurious, and results from the lack of accuracy of BDF2, that widely underestimates the onset of this secondary bifurcation. This again illustrates the dramatic effect of an insufficiently accurate numerical integration in the presence of competing linear instability mechanisms. As for the fact that VMS-BDF2 recovers the proper periodic regime, it is fortuitous and likely originates from the VMS stabilization terms delaying the miscalculated instability threshold, so both errors somehow cancel in the end. This further stresses the quality of the VMS-CN and VMS-BDF4† numerical frameworks, that both accurately predict the main features of the flow oscillations despite the sensitivity of the selection mechanism.

## 5.2. CPU cost

We have mentioned in Sec. 3 that using $P_1-P_1$ elements instead of $P_2-P_1$ elements scales down the number of degrees of freedom (by nearly 70% for this case, from 879,037 with $P_2-P_1$ to 294,087 with $P_1-P_1$). It also considerably cuts the numerical cost, as has been assessed by benchmarking the total memory requirement and CPU time needed to compute 1 time unit of the periodic cavity flow at $Re = 5000$. The required memory is essentially independent of the time discretization and drops by 80%, from 2.6 Gb (DNS) to 0.5 Gb (VMS). In order to smooth over performance differences, the CPU time has been measured from the average over 15 independent runs, with Table 5 (also Fig. 13) providing all
results in arbitrary units to erase the dependency on the hardware resources (e.g., type and number of CPUs). For the DNS, the cost of BDF4† is smaller than that of CN by 38%. This is because the semi-implicit nature of the BDF schemes allows (i) assembling the Jacobian matrix, (ii) preconditioning, and (iii) solving the linear system only once per time step, while Crank-Nicholson requires all operations to be performed at each Newton iterate (which amounts to about 4 times per time step for this case). Moreover, switching from BDF2 to BDF4† comes with little to no extra cost, as the detrimental effect of computing the velocities $u_{BDF\sigma}$ and $u_{BNG\sigma}$ at each time step is essentially paid once BDF1 has been dismissed. Here, the variation is by 4%, which is of the same order as the statistical error (see the vertical bars in Fig. 13 representing the coefficient of variation). The same trends are observed for the VMS, for which the cost of BDF4† is smaller by 39% than that of CN, and barely larger (by 4%) than that of BDF2. More interestingly, if we now compare side-by-side, the cost of VMS-BDF4† is smaller than that of DNS-BDF4† by 32%, and the cost of VMS-CN is smaller than that of DNS-CN by a tremendous 57%. In closing this section, we draw the reader's attention to the fact that the above values are relevant to the present 2D case only. The cost reduction triggered by the VMS in large-scale, 3D problems depends on code-dependent, numerical scalability considerations, and should thus be assessed carefully on a case by a case basis.

### 5.3. Effect of the time step

In the course of the discussion, we have reported additional VMS data obtained with a refined time step $\Delta t = 5 \times 10^{-3}$. Tables 2–4 shows that they remain identical within ~ 1% to those computed with $\Delta t = 10^{-2}$, which is a needed proof of accuracy as the time step is featured in the VMS stabilization parameter $\tau_1$ [15] (for the record, the DNS results exhibit the same level of convergence, but are not reported to ease the reading). The only noticeable difference is with BDF2, whose accuracy improves significantly with a reduced time step, for instance the discrepancy with the VMS-CN results drops to 2.2% at $Re = 5000$ and 6000, and 3.3% at $Re = 4350$. This raises the question about the best interplay between time step and time integration scheme, for which CPU cost is a main point of concern. As seen in Table 5, halving the time step more than doubles the time needed to compute 1 time unit of the periodic
<table><caption>Table 5: CPU time (in arbitrary units) needed to compute 1 time unit of the periodic cavity flow on 1 processor by VMS at Re = 5000. All values proceed from the average over 15 independent runs.</caption><thead><tr><th></th><th>Δt</th><th>CN</th><th>BDF1</th><th>BDF2</th><th>BDF3†</th><th>BDF4†</th></tr></thead><tbody><tr><td>DNS</td><td>10<sup>-2</sup></td><td>1.00</td><td>0.28</td><td>0.38</td><td>0.36</td><td>0.38</td></tr><tr><td>VMS</td><td>10<sup>-2</sup></td><td>0.42</td><td>0.17</td><td>0.25</td><td>0.25</td><td>0.26</td></tr><tr><td>VMS</td><td>5×10<sup>-3</sup></td><td>0.85</td><td>0.40</td><td>0.51</td><td>0.54</td><td>0.55</td></tr></tbody></table>

<caption>Figure 13: Same as Table 5. The black and dark grey bars denote DNS and VMS results computed with time step $10^{-2}$. The light grey bars conversely pertain to VMS results computed with a refined time step $5 \times 10^{-3}$. The vertical error bars stand for the coefficient of variation computed over the 15 independent runs.</caption>
cavity flow at Re = 5000, regardless of the time discretization, which makes the cost of a BDF2 run with $\Delta t = 5 \times 10^{-3}$ twice as large as that of a BDF4† run with $\Delta t = 10^{-2}$. Even worse, it is larger by 20% than that of using CN with $\Delta t = 10^{-2}$, so the advantage of the semi-explicit discretization is essentially lost, and we eventually retain BDF4† with $\Delta t = 10^{-2}$ as the best compromise between numerical accuracy and time efficiency.

## 6. Limit cycle selection

In this section, it is proposed to use the VMS-BDF4† numerical framework to shed some light on the sequence of bifurcation responsible for the limit cycle selection. We report in Fig. 14 the amplitude and frequency of the periodic cavity oscillations, that have been determined by monotonically increasing the Reynolds number from the instability threshold, up to Re = 6500. The change in the dominant mode leads to a quick increase of the oscillation amplitude in Fig. 14(a), but the transition is even more visible in Fig. 14(b), where the oscillation frequency undergoes a discontinuity at Re = 4525 (±25). The latter value is almost identical to that Re$_{CAB}$ = 4553 at which both modes exchange linear dominance (vertical dots in Fig. 14), which suggests that the mode selection essentially proceeds from a 'largest growth rate' criterion.

Matters are not quite so simple, however, because the transition from LCA to mode LCB happens to be sensitive to the initial condition (we recall that all results proceed so far from the same, randomly generated condition). This is further examined using from now on controlled initial conditions

$$u_0 + \Re\{\epsilon_A \hat{u}_A + \epsilon_B \hat{u}_B\}, \quad (20)$$

where $u_0$ is the base cavity flow, i.e., the steady solution to the steady Navier-Stokes equations whose linear instability gives rise to the observed limit cycle oscillations, $\hat{u}_A$ and $\hat{u}_B$ are the first two (complex) eigenmode to bifurcate, computed at the current Reynolds number and normalized to $(\hat{u}, \hat{u}) = 1$, $\Re\{\}$ denotes the real part of a complex quantity, and the $\alpha$ coefficients are (real) perturbation amplitudes chosen small enough for a phase of linear growth to take place before the nonlinear competition between modes sets in. In practice, we use the same solver as in [36] to compute the P$_2$-P$_1$ initial condition, and
<figure><img src="image_34.png" /><figcaption>Figure 14: (a) Amplitude of the x-velocity fluctuations at the S2 sensor, computed against Re with VMS-BDF4†. The symbols (circles and squares for LCA and LCB oscillations, respectively) are the actual data points, and the solid line is the value obtained from (19) after the instability threshold and the amplitude coefficient have been computed by linear regression of the 6 most near-critical squared amplitudes. (b) Oscillation frequency. The vertical dots mark the theoretical threshold value $Re_{cAB} = 4553$ at which modes A and B exchange linear dominance [22].</figcaption></figure>

deduce the $P_1-P_1$ approximation without any interpolation, by retaining the velocity and the pressure at each triangle edge.

Figure 15(a) shows the time evolution of the sensor velocity computed at $Re = 4400$ for a pure A initial perturbation of amplitude $\epsilon_A = 5 \times 10^{-3}$ (hence, $\epsilon_B = 0$). Its counterpart for a pure B initial perturbation with $\epsilon_B = 5 \times 10^{-3}$ and $\epsilon_A = 0$ is provided for comparison in Fig. 15(b). The difference in the initial condition seemingly triggers different transient in Fig. 15(c), as it takes about 75 times units for the periodic regime to emerge starting from mode A, but nearly twice as much starting from mode B. Eventually, both solutions settle on the same limit cycle (up to a phase factor erased in Fig. 15(d)) that proceeds from the saturation of mode A, as can be seen from the detailed statistics in Table 6. The results obtained at $Re = 4700$ are shown in Fig. 16, where both solutions again settle on the same limit cycle, only it now proceeds from the saturation of mode B. Those results obtained at $Re = 4550$ are shown in Fig. 17, where we notice that a pure A initial perturbation selects LCA oscillations, while a pure B initial perturbation selects LCB oscillations. The difference, albeit hardly visible in terms of the oscillation amplitude, is unambiguous in terms of the frequency, as evidenced by the close-up in Fig. 17(d) emphasizing the
<table><caption>Table 6: Statistics of the limit cycle oscillations at Re = 4400, 4550 and 4700.</caption><thead><tr><th>Re</th><th>ε<sub>A</sub></th><th>ε<sub>B</sub></th><th>u&#773;<sub>S2</sub> (×10<sup>-1</sup>)</th><th>u'<sub>S2</sub> (×10<sup>-2</sup>)</th><th>St</th></tr></thead><tbody><tr><td rowspan="3">4400</td><td>5×10<sup>-3</sup></td><td>0</td><td>9.688</td><td>1.618</td><td>1.205</td></tr><tr><td>5×10<sup>-3</sup></td><td>5×10<sup>-3</sup></td><td>9.688</td><td>1.618</td><td>1.205</td></tr><tr><td>0</td><td>5×10<sup>-3</sup></td><td>9.688</td><td>1.618</td><td>1.205</td></tr><tr><td rowspan="3">4550</td><td>5×10<sup>-3</sup></td><td>0</td><td>9.716</td><td>1.948</td><td>1.212</td></tr><tr><td>5×10<sup>-3</sup></td><td>5×10<sup>-3</sup></td><td>9.721</td><td>1.989</td><td>1.658</td></tr><tr><td>0</td><td>5×10<sup>-3</sup></td><td>9.721</td><td>1.989</td><td>1.658</td></tr><tr><td rowspan="3">4700</td><td>5×10<sup>-3</sup></td><td>0</td><td>9.730</td><td>2.781</td><td>1.668</td></tr><tr><td>5×10<sup>-3</sup></td><td>5×10<sup>-3</sup></td><td>9.730</td><td>2.781</td><td>1.668</td></tr><tr><td>0</td><td>5×10<sup>-3</sup></td><td>9.730</td><td>2.781</td><td>1.668</td></tr></tbody></table>

related periods $\tau = 1/\text{St}$. The selection proceeds from both linear and nonlinear considerations, which is best seen in Fig. 18 showing the results pertaining to various combinations of A/B initial disturbances. The general picture is that LCA is selected if the initial amplitude of mode B is smaller than that of mode A by at least an order of magnitude, otherwise LCB is selected even though the initial amplitude of mode B is as small as $5 \times 10^{-5}$ and mode A is more linearly unstable at this Reynolds number.

The lower and upper bounds of the hysteresis domain have been determined numerically to be $Re_{cA→B} = 4455 \pm 5$ and $Re_{cB→A} = 4635 \pm 5$. These values proceed from about 50 runs performed at several Reynolds in the range [4400; 4700]. A particular attention has been paid to assessing systematically that the system ultimately selects LCA or LCB. No other oscillation pattern has been found, which is best exemplified in Fig. 19(a) showing the time evolution of the sensor velocity computed at Re = 4650, starting from a pure A initial disturbance with $\epsilon_A = 5 \times 10^{-3}$. The system exhibits a long transient reflecting a complex competition between modes, and we show in Fig. 19(b) that the single-sided amplitude spectrum, computed at Re = 4650 over the first 300 time units does feature both modes oscillating with finite amplitude at their incommensurate frequencies. This quasi-periodic solution is not stable,
<figure><img src="image_35.png" /><figcaption>Figure 15: (a,b) Time evolution of the x-velocity at the S2 sensor, computed at Re = 4400 with VMS-BDF4†, starting from (a) a pure A initial disturbance with $\epsilon_A = 5 \times 10^{-3}$ ($\epsilon_B = 0$), and (b) a pure B initial disturbance with $\epsilon_B = 5 \times 10^{-3}$ ($\epsilon_A = 0$). (c,d) Close-ups on the (c) transient and (d) periodic regimes. The results pertaining to modes A and B are reported from (a,b) as the black and grey lines, respectively. The considered time intervals are marked in (a,b) by the grey boxes.</figcaption></figure>

however, as the spectrum computed over the next 300 time units are pure LCB oscillations; see Fig. 19(c). Such transients can last close to 1000 time units if the Reynolds number approaches the cusp of the hysteresis domain; see the results obtained at Re = 4640 in Fig. 19(d), hence the need for a highly efficient numerical framework.

In the terminology of the normal form theory, the observed hysteretic behavior is typical of a subcritical Neimark-Sacker bifurcation. The subcritical nature has especially been confirmed by perturbing numerically the limit cycle solution. In Fig. 20(a), we start from the same pure B initial disturbance with $\epsilon_B = 5 \times 10^{-3}$ as above, that eventually selects LCB. At $t = 150$, we substitute

$$u_{\text{LCB}} + \epsilon \tilde{u}, \qquad (21)$$
<figure><img src="image_36.png" /><figcaption>Figure 16: Same as Fig. 15 at Re = 4700.</figcaption></figure>

for the limit cycle velocity, where $\tilde{u}$ is a random velocity field with controlled amplitude $\epsilon$, hence normalized to $(\tilde{u}, \tilde{u}) = 1$, and we time march the system until it settles again on a limit cycle. The nominal, unperturbed results are displayed in black, and the perturbed results in grey. For a small value $\epsilon = 5 \times 10^{-3}$ of the amplitude, the LCB oscillations are promptly restored. More surprisingly, the same behavior is observed with an amplitude as large as $\epsilon = 10^{-1}$, as we show in Figs. 20(c) and 20(d) that LCB is restored after a short transient of about 20 time units. We have also considered perturbing with a pure A disturbance of controlled amplitude $\epsilon_A$, as achieved substituting instead

$$u_{\text{LCB}} + \epsilon_A \hat{u}_A, \qquad (22)$$

for the limit cycle velocity. For small perturbation amplitudes, the LCB oscillations are again promptly restored, as evidenced in Figs. 21(a) and 21(b) for $\epsilon_A = 5 \times 10^{-3}$. A different behavior is observed using a large amplitude $\epsilon_A = 10^{-1}$, as Figs. 21(c) and 21(d) exhibit a long transient of about 150 time units during which both modes compete against one another, before the even-
<figure><img src="image_37.png" /></figure>

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

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

<figure><img src="image_40.png" /><figcaption>Figure 17: Same as Fig. 15 at Re = 4550.</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>Figure 18: Limit cycle computed at Re = 4550 with VMS-BDF4†, starting from various blends of A/B initial disturbances. The circle and square symbols denote LCA and LCB oscillations, respectively, with the grey shade to visually enhance the LCB domain.</figcaption></figure>

tual outbreak of LCA. Similar results have been obtained perturbing LCA with either random or pure B disturbances, meaning that LCB (resp. LCA) is destabilized only if mode A (resp. mode B) is applied to the system with a sufficiently large amplitude, otherwise the small amount of mode A (resp. mode B) present in a random disturbance is quickly damped nonlinearly by the Reynolds stress of the saturated mode B (resp. mode A).

495
<figure><img src="image_42.png" /><figcaption>Figure 19: (a) Time evolution of the x-velocity at the S2 sensor, computed at Re = 4650 with VMS-BDF4†, starting from a pure A initial disturbance with $\epsilon_A = 5 \times 10^{-3}$. (c,d) Single-sided amplitude spectrum of the velocity fluctuations, computed over the time intervals (c) $t \in [0; 300]$, and (d) $t \in [300; 600]$. (e) Same as (a) at Re = 4640.</figcaption></figure>

7. Conclusion

This research uses the VMS modeling of the Navier-Stokes equations to compute numerical solutions of the periodic flow over an open cavity. The space discretization uses linear approximations ($P_1$ finite elements) for both the velocity and the pressure, which breaks the Babuska-Brezzi condition and is thus bound to fail in a DNS. The solution is marched in time with two discretization schemes, a semi-implicit BDF scheme based on backward differentiation formulas, and the implicit Crank-Nicholson scheme.

A detailed comparison with DNS data shows that the VMS-BDF4† and VMS-CN numerical frameworks accurately compute the limit cycle oscillations at several Reynolds numbers $Re = 4350, 5000$ and $6000$ above the instability threshold, as well as the near-critical dynamics. By doing so, the CPU cost is reduced by ~35% using BDF4† (resp. ~60% using CN), while the memory load
<figure><img src="image_43.png" /><figcaption>Figure 20: (a) Time evolution of the x-velocity at the S2 sensor, computed at $Re = 4550$ with VMS-BDF4†. The black line is for the pure B initial disturbance with $\epsilon_B = 5 \times 10^{-3}$, found in Fig. 17(a) to yield LCB oscillations. The grey line is for perturbed results obtained by superimposing at $t = 150$ (vertical dots) a random velocity field of controlled amplitude $\epsilon = 5 \times 10^{-3}$. (b) Close-up on the periodic regime. The considered time interval is marked in (a) as the grey box. (c,d) Same as (a,b) with $\epsilon = 10^{-1}$.</figcaption></figure>

is alleviated by ~80%, regardless of the time-discretization scheme. This is all the more remarkable given that the cavity flow undergoes two Hopf bifurcations close one to the other, at $Re_{cA} = 4130$ and $Re_{cB} = 4349$, in a way such that there exist two competing instabilities (termed here A and B) in this range of Reynolds numbers. In these circumstances, even a small miscalculation of the related eigenvalues (whether it be because the time discretization scheme lacks accuracy or because the VMS stabilization terms can accidentally shift the eigenspectrum of the linearized Navier-Stokes operator) can affect the linear instability thresholds, the nonlinear interaction between unstable modes, and eventually, the frequency selection.

Finally, the highly efficient VMS-BDF4† framework is used to analyze the mechanism underlying the selection of the limit cycle frequency. Using con-
<figure><img src="image_44.png" /><figcaption>Figure 21: Same as Fig. 20 superposing at $t = 150$ pure A disturbances of controlled amplitude (a,b) $\epsilon_A = 5 \times 10^{-3}$ and (c,d) $\epsilon_A = 10^{-1}$.</figcaption></figure>

trolled initial conditions built from the unstable modes, we show that the latter involves a hysteresis (and thus a discontinuity of the frequency curve against Re) reflecting the occurrence of a subcritical Neimark-Sacker bifurcation. The lower and upper bounds of the hysteresis domain are determined numerically to be $Re_{cA \to B} = 4455 \pm 5$ and $Re_{cB \to A} = 4635 \pm 5$. The system is shown to switch from one limit cycle to another (say from LCB to LCA) only if the competing mode (hence, mode A) is applied with a sufficiently large amplitude.

Acknowledgements

This work is supported by the «Investissements d'Avenir» French Government program, 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). It was granted access to the HPC resources of Aix-Marseille University financed by the project Equip@Meso (ANR-10- EQPX-29-01).
References

[1] T. J. R. Hughes, G. R. Feijóo, L. Mazzei, J.-B. Quincy, The variational multiscale method - a paradigm for computational mechanics, Comput. Methods Appl. Mech. Engrg. 166 (1998) 3-24.

[2] R. Codina, Stabilization of incompressibility and convection through orthogonal sub-scales in finite element methods, Comput. Methods Appl. Mech. Engrg. 190 (2000) 1579-1599.

[3] Y. Bazilevs, V. M. Calo, J. A. Cottrell, T. J. R. Hughes, A. Reali, G. Scovazzi, Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows, Comput. Methods Appl. Mech. Engrg. 197 (2007) 173-201.

[4] R. Codina, J. Principe, Dynamic subscales in the finite element approximation of thermally coupled incompressible flows, Int. J. Numer. Meth. Fl. 54 (2007) 707-730.

[5] G. Scovazzi, Lagrangian shock hydrodynamics on tetrahedral meshes: a stable and accurate variational multiscale approach, J. Comput. Phys. 231 (2012) 8029-8069.

[6] U. Rasthofer, V. Gravemeier, Multifractal subgrid-scale modeling within a variational multiscale method for large-eddy simulation of turbulent flow, J. Comput. Phys. 234 (2013) 79-107.

[7] T. Coupez, E. Hachem, Solution of high-reynolds incompressible flow with stabilized finite element and adaptive anisotropic meshing, Comput. Methods Appl. Mech. Engrg. 267 (2013) 65-85.

[8] R. Calderer, A. Masud, Residual-based variational multiscale turbulence models for unstructured tetrahedral meshes, Comput. Methods Appl. Mech. Engrg. 254 (2013) 238-253.

[9] T. Tezduyar, Y. Osawa, Finite element stabilization parameters computed from element matrices and vectors, Comput. Methods Appl. Mech. Engrg. 190 (3-4) (2000) 411-430.
565

[10] T. J. R. Hughes, G. Scovazzi, P. B. Bochev, A. Buffa, A multiscale discontinuous galerkin method with the computational structure of a continuous galerkin method, Comput. Methods Appl. Mech. Engrg. 195 (2006) 2761–2787.

570

[11] R. Codina, Stabilized finite element approximation of transient incompressible flows using orthogonal subscales, Comput. Methods Appl. Mech. Engrg. 191 (2002) 4295–4321.

[12] C. Farhat, A. Rajasekharana, B. Koobus, A dynamic variational multiscale method for large eddy simulations on unstructured meshes, Comput. Methods Appl. Mech. Engrg. 195 (2006) 1667–1691.

575

[13] K. B. Nakshatrala, D. Z. Turner, K. D. Hjelmstad, A. Masud, A stabilized mixed finite element method for darcy flow based on the multiscale decomposition of the solution, Comput. Methods Appl. Mech. Engrg. 196 (2006) 120–142.

580

[14] I. Akkerman, Y. Bazilevs, V. M. Calo, T. J. R. Hughes, S. Hulshoff, The role of continuity in residual-based variational multiscale modeling of turbulence, Comput. Mech. 41 (2008) 371–378.

[15] E. Hachem, S. Feghali, R. Codina, T. Coupez, Immersed stress method for fluid-structure interaction using anisotropic mesh adaptation, Int. J. Numer. Meth. Eng. 94 (2013) 805–825.

585

[16] S. C. Reddy, P. J. Schmid, D. S. Henningson, Pseudospectra of the orrisommerfeld operator, SIAM J. Appl. Math. 53 (1) (1993) 15–47.

[17] L. N. Trefethen, A. E. Trefethen, S. C. Reddy, T. A. Driscoll, Hydrodynamic stability without eigenvalues, Science 261 (5121) (1993) 578–584.

590

[18] J.-M. Chomaz, Global instabilities in spatially developing flows: Nonnormality and nonlinearity, Annu. Rev. Fluid Mech. 37 (2005) 357–392.

[19] O. Colomès, S. Badia, R. Codina, J. Principe, Assessment of variational multiscale models for the large eddy simulation of turbulent incompressible flows, Comput. Methods Appl. Mech. Engrg. 285 (2015) 32–63.
595

[20] D. Forti, L. Dedè, Semi-implicit BDF time discretization of the Navier-Stokes equations with VMS-LES modeling in a high performance computing framework, *Comp. Fluids* 117 (2015) 168-182.

600

[21] D. Sipp, A. Lebedev, Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows, *J. Fluid Mech.* 593 (2007) 333-358.

605

[22] P. Meliga, Harmonics generation and the mechanics of saturation in flow over an open cavity: a second-order, self-consistent description, *J. Fluid Mech.* 826 (2017) 503-521.

[23] C. Mathis, M. Provansal, L. Boyer, Bénard-von kármán instability: an experimental study near the threshold, *J. Physique Lett.* 45 (1984) 483-491.

[24] M. Provansal, C. Mathis, L. Boyer, Bénard-von kármán instability: transient and forced regimes, *J. Fluid Mech.* 182 (1987) 1-22.

[25] C. P. Jackson, A finite-element study of the onset of vortex shedding in flow past variously shaped bodies, *J. Fluid Mech.* 270 (1987) 23-45.

610

[26] A. Zebib, Stability of viscous flow past a circular cylinder, *J. Engrg. Maths* 21 (1987) 155-165.

[27] D. Barkley, M. G. M. Gomes, R. D. Henderson, Three-dimensional instability in flow over a backward-facing step, *J. Fluid Mech.* 473 (2002) 167-190.

615

[28] O. Marquet, M. Lombardi, J.-M. Chomaz, D. Sipp, L. Jacquin, Direct and adjoint global modes of a recirculation bubble: lift-up and convective non-normalities, *J. Fluid Mech.* 622 (2009) 1-21.

[29] E. Åkervik, U. Ehrenstein, F. Gallaire, D. S. Henningson, Global two-dimensional stability measures of the flat plate boundary-layer flow, *Eur. J. Mech. B-Fluid* 27 (5) (2008) 501-513.

620

[30] A. Barbagallo, D. Sipp, P. J. Schmid, Closed-loop control of an open cavity flow using reduced-order models, *J. Fluid Mech.* 641 (2009) 1-50.
[31] E. Hachem, B. Rivaux, T. Kloczko, H. Digonnet, T. Coupez, Stabilized finite element method for incompressible flows with high Reynolds number, J. Comput. Phys. 229 (23) (2010) 8643–8665.

625

[32] F. E. Cellier, E. Kofman, Continuous system simulation, New York: Springer-Verlag, 2006.

[33] G. S. Rao, Numerical analysis., New Delhi: New Age International Publishers, 2009.

630

[34] V. N. Vatsa, M. H. Carpenter, Higher order temporal schemes with error controllers for unsteady Navier–Stokes equations, AIAA 2005-5245.

[35] V. N. Vatsa, M. H. Carpenter, D. P. Lockhard, Re-evaluation of an optimized second order backward difference (BDF2OPT) scheme for unsteady flow applications, AIAA 2010-0122.

635

[36] P. Meliga, J.-M. Chomaz, An asymptotic expansion for the vortex-induced vibrations of a circular cylinder, J. Fluid Mech. 671 (2011) 137–167.

[37] F. Alizard, X. Merle, J.-C. Robinet, X. Gloerfelt, Linear and non-linear dynamics of a cavity flow, presented at the 8th SIG 33 ERCOFTAC workshop on Global Instabilities of open flows, Nice, France (2010).