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

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 similar to 35% (resp. similar to 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 similar to 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.


Introduction
The variational multiscale (VMS) modeling of the Navier-Stokes equations [1][2][3][4][5] has been widely developed for the numerical simulation of turbulent flows in several benchmarking and applicative context [6][7][8]. While the direct numerical simulation (DNS) requires a complete representation of the whole range of 5 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 10 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 en- 15 gineering 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 ac- 20 comodating 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 25 defining stable spatial approximation schemes, while also suitably representing the turbulence, which is now supported by an extensive body of work [11][12][13][14][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 30 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][17][18]. While this is not necessarily visible at high Reynolds num-2 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 The retained configuration is the open cavity flow documented in [21], that becomes unstable at Reynolds number Re = 4140 (based on the free-stream ve-45 locity and the cavity height), then exhibits a well-defined oscillatory behavior up to Re ∼ 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 50 much higher than in other classically benchmarked configurations, e.g., cylinder flows, whose critical Reynolds number Re ∼ 50 (based on the free-stream velocity and the cylinder diameter) is smaller by two orders of magnitude [23][24][25][26], or backwards facing step flows, whose critical Reynolds number Re ∼ 800 (based on the centerline velocity and the step height) is smaller by almost one order 55 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 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 70 limit cycle oscillations.

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 75 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 80 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 largescale vortices via Kelvin-Helmholtz instability. For sufficiently large Reynolds numbers, this feedback loop leads to a linear global instability, whose nonlinear 85 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 with u =( u, v) the velocity field of component u (resp. v) in the x (resp. y) direction, p the pressure, I and ε(u) the identity and the rate of deformation tensors, and Re the Reynolds number. We use the free-stream velocity and the cavity height to make all quantities non-dimensional. The upstream and Finally, a no-slip condition u = v = 0 is imposed on the remainder of the cavity walls.

Numerical schemes
The above computational domain is denoted by Ω, with V ⊂ H 1 (Ω) 2 and 100 Q = L 2 (Ω) 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).

Finite elements discretization 105
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 with w (resp. q) an arbitrary test function in V (resp. Q), and ( • , • ) the L 2 inner product on Ω. The Babuska-Brezzi condition imposes to discretize V and 120 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 125 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 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].
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 where ( • , • ) K is the inner product on element K, R C,M being the residuals defined by and τ 1,2 are ad-hoc mesh-dependent stabilization parameters defined in [11,15] 130 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. is set to 10 −6 in terms of both the absolute and relative norm of the residual.

Time discretization
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]. 145 We consider first the implicit Crank-Nicholson (CN) scheme

Implicit Crank-Nicholson
where the i superscript refers to the solution at time t i = i∆t. A Newton root finding algorithm is used to handle the nonlinearity of (5) in u i : at the j-th All VMS stabilization terms in (3) are integrated explicitely with the forward Euler scheme. This yields the fully discretized stabilized weak form with residuals defined as and the weak form for the DNS deduces by forcing τ i,j−1 1,2 to zero.

Semi-implicit BDF
Alternatively, we consider the semi-implicit scheme where the time and convective derivatives are approximated with backward differentiation formulas (BDF) and backward Newton-Gregory polynomials (BNG) of same order σ, respectively. This yields where {β k ,γ k } are sequences defined recursively after [32,33] by hence β k>σ and γ k>σ 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 with residuals defined as and the weak form for the DNS deduces by forcing τ i,j−1 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., where the χ σ ′ coefficients reported in Table 1

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 x S2 =( 0 .75, 0.05) in the shear layer is marked by 160 the green circle in Fig. 1. For benchmarking purposes, we also report the values measured by two additional sensors S1 and S3, whose positions x S1 =(0.6, 0.05) and x S3 =(0.9, 0.05) are marked by the red and blue circles, respectively. The time step ∆t =1 0 −2 is the same as in [22], and allows sampling between ∼ 55 and 85 points over a period, depending on the Reynolds number. We consider a  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  ited discrepancies at each time step, that trigger a small phase shift in the limit cycle solution 2 , but the latter is erased in Fig. 2(c) to ease the presentation.

Discretization ∆t
Other than that, the limit cycle statistics in Table 2 clearly assess the accuracy 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. 185 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 190 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 195 the optimized BDF schemes improve significantly the accuracy, namely BDF3 † All results pertain to the same phase in the limit cycle, corresponding to a peak velocity at the S2 sensor.
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.  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 230 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 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 240 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 overes-  All results pertain to the same phase in the limit cycle, corresponding to a peak velocity at the S2 sensor.
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 250 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 255 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 ver-260 tical cuts at x =0 .9 in Fig. 7(d) for which the shear layer velocity is strongly distorted by the inner vortices.  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.

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 265 Re = 5000. As shown in Figure 9 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   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 squareroot of the departure from criticality, i.e.,

Discretization ∆t
with Re c the critical Reynolds number and ξ an amplitude coefficient accounting    (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.
into the technical details, those are free from time-discretization errors, P 2 -305 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 is smaller than the error on the amplitude by one order of magnitude) and an amplitude coefficient ξ =4.24 overestimating the theoretical value by as few as

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 10 −3 . The present study conversely uses a value 10 −6 . and Re cB = 4349. Interestingly, the first eigenmode to bifurcate, termed A, 320 remains more linearly unstable up to Re cAB = 4553, after which mode B takes over. Those thresholds (just as all eigenfrequencies ω/2π 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 bechmark results can thus be interpreted as follows:

325
• 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 eigenfrequency ω A /2π =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 330 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.

335
• 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 ω B /2π =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 340 matches well the eigenfrequency ω A /2π =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 345 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 ω B /2π =1.70 agrees well with the limit cycle frequency St B =1 .78. Note, the statistics of the BDF1 limit  375 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 380 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 385 (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 390 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σ and u BNGσ 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 395 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 400 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. 405 In the course of the discussion, we have reported additional VMS data obtained with a refined time step ∆t =5× 10 −3 . Tables 2-4 shows that they remain identical within ∼ 1% to those computed with ∆t =1 0 −2 , which is a needed proof of accuracy as the time step is featured in the VMS stabilization parameter τ 1 [15] (for the record, the DNS results exhibit the same level of 410 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 415 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  27 cavity flow at Re = 5000, regardless of the time discretization, which makes the cost of a BDF2 run with ∆t =5× 10 −3 twice as large as that of a BDF4 † run with ∆t =1 0 −2 . Even worse, it is larger by 20% than that of using CN with 420 ∆t =1 0 −2 , so the advantage of the semi-explicit discretization is essentially lost, and we eventually retain BDF4 † with ∆t =1 0 −2 as the best compromise between numerical accuracy and time efficiency.

Limit cycle selection
In this section, it is proposed to use the VMS -BDF4 † numerical framework 425 to shed some light on the sequence of bifurcation responsible for the limit cycle selection. We report in Fig. 14 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.

435
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 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,û A andû B are the first two (complex) eigenmode to bifurcate, computed at the current Reynolds number and normalized to (û,û)=1 , { } denotes the real part of a complex quantity, and the α coefficients are (real) 440 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  Fig. 15(b). The difference in the initial condition seemingly triggers different transient in Fig. 15(c), as it takes about 75 450 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 455 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 460 of the frequency, as evidenced by the close-up in Fig. 17(d) emphasizing the 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 475 initial disturbance with A =5× 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  Fig. 19(d), hence the need for a highly efficient numerical framework.

485
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 B =5× 10 −3 as above, that eventually selects LCB. At t = 150, we substitute for the limit cycle velocity, whereũ is a random velocity field with controlled amplitude , hence normalized to (ũ,ũ) = 1, and we time march the system until it settles again on a limit cycle.  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 495 of the saturated mode B (resp. mode A).

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 ve-500 locity 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 505 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 is alleviated by ∼ 80%, regardless of the time-discretization scheme. This is all 510 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 515 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 520 mechanism underlying the selection of the limit cycle frequency. Using con- be Re cA→B = 4455 ± 5a n dR e cB→A = 4635 ± 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.