Comptes Rendus Mécanique

. This paper revisits oblique wave and streamwise vortex scenarios in a plane Couette flow using restricted nonlinear simulations, where only a single Fourier mode for perturbation is retained. It is shown that this restriction of full dynamics gives a good approximation of the two subcritical paths. In particular, critical energy thresholds and edge states compare favorably with results obtained using direct numerical simulations by Duguet et al. ( Phys. Rev. E 82 (2010), 026316).


Introduction
Laminar/turbulent transition of wall-bounded shear flows is responsible for a dramatic increase in skin-friction drag for a wide variety of industrial applications.In this regard, the determination of the initial perturbation that leads to transition is fundamental to designing efficient control strategies.
The transition to turbulence often occurs even in the absence of an unstable linear mode associated with the laminar basic state.In this subcritical scenario, both the initial amplitude and the spatial shape of the initial perturbation play prominent roles in triggering turbulence.To gain an improved understanding of the subcritical transition process, researchers have mainly focused on a small number of incompressible canonical flows.
For instance, Kreiss et al. [1] and Reddy et al. [2] first investigated the critical energy thresholds in both plane Couette flow (pCf) and channel flow for the streamwise vortex (SV) and oblique wave (OW) scenarios using direct numerical simulations (DNSs).For the SV scenario, the transition is initiated by streamwise vortices that generate streaks due to the lift-up effect [3].When secondary perturbations are superimposed, the streak experiences secondary instability and breaks down to turbulence.The second route to turbulence begins with a pair of OWs that interact nonlinearly to create streamwise vortices.This leads to the formation of streaks that break down to turbulence due to secondary instability.For the pCf, Reddy et al. [2] found that the threshold energy density associated with the initial perturbation scales as Re −2.2 and Re −2.5 , where Re is the Reynolds number, for the SV and OW scenarios, respectively.Later, Duguet et al. [4] corrected the values given by Reddy et al. [2].In particular, the authors found that the critical energy threshold scales as Re −2 for both OW and SV paths for the pCf.These scaling factors are consistent with the theoretical analysis by Chapman [5].More recently, Karp et al. [6,7] showed that the optimal set of parameters for transition in the pCf is associated with initial disturbances, which generate the strongest inflection point in the mean flow profile when nonlinearities come into play.In the same vein, Cossu et al. [8] investigated the streak breakdown in a two-dimensional parameter space with respect to the amplitudes of the low-speed streak and its secondary perturbation.The authors showed that the transition is mainly associated with finite-amplitude disturbances that push the flow state along an unstable direction of the edge of chaos (Skufca [9]).In the same spirit, Duguet et al. [4] carried out optimizations to seek a minimal initial perturbation based on a linear combination of a finite number of linear optimal modes.When retaining three linear optimal modes, the authors found that the energy threshold of such a minimal perturbation is only 2% less than that for a pair of symmetric OWs.For both SV and OW scenarios, they found that the trajectories closely approach two different steady states, the so-called edge states, with special symmetries.
Over the past decade, a great deal of effort has been made in finding the weakest disturbance that can touch the edge of chaos.This perturbation is referred to as the "minimal seed" and has the lowest energy on the edge.Using nonlinear optimizations, Monokrousos et al. [10], Rabin et al. [11], and Duguet et al. [12] identified a localized minimal disturbance in the pCf.In particular, Duguet et al. [12] proposed a new scaling law for the minimal energy perturbation O(Re −2.7 ) well below the energy threshold for the OW scenario O(Re −2 ).The temporal evolution of the minimal perturbation exhibits well-known linear mechanisms (Orr, lift-up, and secondary streak instability), which the nonlinearity of the Navier-Stokes equations couples together [13].Therefore, several attempts have been made to develop simplified models that are able to describe the main steps of the subcritical transition.A common aspect of all simplified models relies on the decomposition of the instantaneous flow fields into a streamwise-averaged flow and a fluctuation.Nonlinearities are removed from the fluctuation equation, and a nonlinear feedback term is present in the streamwise-averaged part.For instance, using a single streamwise Fourier mode for the fluctuation and within the framework of the vortex-wave interaction (VWI) theory [14], Deguchi et al. and Deguchi [15,16] established a Reynolds number scaling for transition in the pCf: Re −2.5 , which is close to that found using DNS (Re −2.7 ).One may recall that the VWI theory is consistent with the self-sustained process for asymptotically large Reynolds numbers [17,18].Biau and Bottaro [19] carried out nonlinear optimizations based on a similar decomposition, the so-called restricted nonlinear (RNL) system [20], to identify an optimal path to transition in a linearly stable duct flow.Similarly, Pralits et al. [21] showed that nonlinear optimizations based on an RNL system with a single streamwise Fourier mode and two spanwise Fourier modes are accurate for finding an optimal path according to the OW scenario.Nevertheless, none of the studies were devoted to clearly establishing the ability of such RNL models to describe the OW and SV paths.In particular, are we able to recover the scaling laws found by Duguet et al. [4] using DNS as suggested by the theoretical approach of Chapman [5]?In particular, how are the edge states associated with OW and SV paths captured by RNL models?
The aim of this paper is to shed some light on these two points.For this purpose, the pCf case studied by Duguet et al. [4] is investigated.A single Fourier mode for fluctuation will be retained hereafter.The paper is organized as follows.In Section 2, we recall the equations of the RNL model.After presenting numerical methods and the decomposition of the initial perturbation as a finite number of linear optimal modes, the results are given in Section 3. Specific attention is devoted to providing the Reynolds number scaling for both scenarios within an RNL framework and revealing the corresponding edge states.Finally, Section 4 is devoted to conclusions and prospects.

Restricted nonlinear model
We introduce the Cartesian coordinate system (x, y, z) and consider the pCf of an incompressible fluid with kinematic viscosity ν between two parallel plates located at y = ±h.The two plates move in opposite directions with velocity (±U 0 , 0, 0).The Reynolds number is defined as Re = U 0 h/ν.We introduce the streamwise averaging operator over a streamwise distance L x for any quantity ψ: The velocity disturbance to the Couette laminar solution (U b = y, 0, 0) is decomposed into its streamwise-averaged part U = u, called mean flow hereafter, and its fluctuation u: u(x, y, z, t ) = U(y, z, t ) + u(x, y, z, t ). ( The components u = (u, v, w) are streamwise, wall-normal, and spanwise velocities, respectively.The pressure is similarly split as follows: p(x, y, z, t ) = P (y, z, t ) + p(x, y, z, t ).
Hereafter, we will consider a minimal RNL system where only a single streamwise component is retained.The fluctuations are thus expressed as where α is the streamwise wavenumber (α = 2π/L x ) and * is used to denote the complex conjugate.The system of equations governing the streamwise-averaged flow reads as which are associated with homogeneous boundary conditions U = V = W = 0 at the walls.The equations for the fluctuation linearized around the mean flow are together with û = v = ŵ = 0 at the walls.The coupled system (3), ( 4) is referred to as the RNL model.

Numerical methods and flow cases
For numerical integration of the RNL system, spectral approximations of the velocity fields are adopted using Fourier expansions in both streamwise and spanwise directions and Chebyshev polynomials in the wall-normal direction.The numerical methods were described by Alizard [22].The code was also validated by Alizard [23].

Initial perturbations
Following previous studies based on DNS [2,4], the initial perturbation is constructed by the summation of linear optimal modes.A linear optimal mode is defined as the initial perturbation that maximizes kinetic energy gain over all times when the linearized Navier-Stokes operator is advanced forward in time [24].For this, we have where each mode U i corresponds to an optimal mode with α = 0 and different harmonics in z.Modes u i correspond to optimal modes with α = 0.5 (L x = 4π) and different harmonics in z.The complex amplitudes of the different modes are denoted as A i and a i .
The kinetic energy threshold for transition is determined using a bisection algorithm.In this regard, we introduce the scalar function as in Duguet et al. [4], where N and M modes are considered for the streamwise-averaged quantity and fluctuation, respectively.
Hereafter, Fourier modes corresponding to the wavevector components (nα, mβ) with (α, β) = (2π/L x , 2π/L z ) are labeled by (m, n).In the present study, we restrict our analysis to two different cases.For the oblique path, a pair of OWs with equal amplitude (1, ±1) is chosen as the initial perturbation.For the OW scenario, the only unknown is the amplitude of the perturbation, which is fixed through the bisection algorithm.For the SV scenario, the initial subspace for the perturbation is reduced to two optimal modes (0, 2) and (1, 2).Considering the symmetries in x and z, the amplitude coefficients ℜ(A 1 ), ℑ(A 1 ), ℜ(a 1 ), ℑ(a 1 ), . . .are reduced to one parameter denoted as λ.The integrated kinetic energies for the perturbation and the mean flow are denoted as E p (t ) and E M (t ), respectively.They are defined as Without loss of generality, λ is fixed as E p (t = 0)/E M (t = 0).In this case, the expression of E c is reduced to E c (λ).In Duguet et al. [4], the minimization of E c is carried out using a Newton algorithm.Here, we investigate the distribution of E c with λ.The stopping criterion of the bisection algorithm is fixed as for a given iteration i .The critical energy threshold includes both contributions for the SV scenario (E M (t = 0) + E p (t = 0)) and only the kinetic energy of the fluctuation for the OW case (E p (t = 0)).

Kinetic energy thresholds
The results of RNL simulations corresponding to the pure OW and SV scenarios are now presented in parallel.In Figure 1, we show the integrated kinetic energy of the perturbation (both fluctuation and mean flow contributions) as a function of time when the initial amplitude is varied.For the purpose of illustration, the Reynolds number is fixed as Re = 1500.The figure shows that for the two cases, the bisection algorithm converges to an almost steady state, where E p (t ) and E M (t ) reach a plateau before breakdown.At convergence, the critical energy threshold is E c = E M (t = 0) + E p (t = 0).One may also note that for an initial kinetic energy slightly above E c , the integration forward in time of the initial perturbation leads to a self-sustained state for long times.
The variations in threshold energy E c with the Reynolds number are shown in Figure 2 for the pure OW scenario and the SV scenario.The minimization of E c is used to fix λ for the SV case. Figure 2 shows that the exponent predicted by the RNL model is in perfect agreement with the full DNS for both scenarios (i.e., −2).In addition, fitted curves of the kinetic energy thresholds obtained using either the RNL system or the full nonlinear system are nearly overlapping.This finding is in agreement with the theoretical development by Chapman [5].

Figure 2. Energy threshold versus the Reynolds number
Re for Re = 750, 1500, 2000, 3000, and 3500 obtained by using RNL simulations (symbols).The fitted curves for oblique waves and streamwise vortex scenarios using RNL simulations (full lines) are compared with values provided by [4] using fully nonlinear simulations (dashed lines).The Re-scaling law for the OW scenario is 4.5Re −2.03 .For the SV scenario, the Re-scaling law is 23.3Re −1.99 .To better appreciate the relevance of the strong restriction associated with the RNL model in describing the OW path, the variation in the threshold kinetic energy with the streamwise computational domain size is shown in Figure 3 for Re = 400.Results provided by Duguet et al. [4] are also illustrated.We recall that for this case, we consider (N y , N z ) = (33, 48).The OW structure leading to transition with the lowest energy density has a streamwise wavelength L x ≈ 6π when considering the RNL system.The full DNS gives L x ≈ 5.5π [4].This corresponds to an angle with respect to the streamwise direction of 20 • for the DNS and 18.5 • for the RNL model.This illustrates that the RNL model is quite accurate in providing the direction of the most efficient OW.Furthermore, the kinetic energy threshold for the optimal wavelength obtained using RNL simulations (2.09 × 10 −5 ) is also quite close to that obtained by DNS (2.92 × 10 −5 ).

Edge states
In the previous section, it was shown that for both OW and SV paths, when the initial kinetic energy is fixed as E c , the trajectory passes close to a steady state.It is suggested that such a state corresponds to an edge state.To better confirm this assumption, we consider hereafter the same flow case as that in Duguet et al. [4] (Re = 400).To find the corresponding edge state, the stopping criterion associated with the bisection algorithm is fixed as r = 10 −6 .The edge state found using a pair of OWs with equal amplitudes is shown in Figure 4.This steady state exhibits the same features that the edge state referenced as E 1 by Duguet et al. [4].In particular, E 1 displays lowspeed streaks with a varicose structure in both simulations.Both states lie in the same spatial symmetry subspace associated with the following transformations (the notation from Duguet et al. [4] is used), as well as any combination of these symmetries.
For the SV scenario, the initial perturbation is driven by modes (0, 2) and (1, 2).In Figure 5, the distribution of the critical initial kinetic energy as a function of λ is displayed.The figure shows that an optimal value of λ is approximately 0.14.For the full DNS, Duguet et al. [4] found λ = 0.0965, which is quite close to that provided by the RNL approximation.Furthermore, the figure also shows that the improvement in terms of critical energy is weak (only 5% for the range of λ that is investigated).This remark is consistent with the observations by Duguet et al. [4].The edge state is shown in Figure 6.This edge state shares strong similarities with the state referred to as E 2 in [4].In particular, E 2 exhibits a sinuous motion for a pair of low-and high-speed streaks  and has the following symmetries (the notation from Duguet et al. [4] is used): Hence, phase space trajectories are seen to be also well reproduced for both OW and SV scenarios.

Conclusions
This study revisits the OW and SV scenarios for subcritical transition in a pCf within an RNL framework.It is shown that the strong restriction of the full dynamics associated with the RNL system using a single streamwise Fourier mode for perturbation allows us to describe well both subcritical paths.Mainly, the critical energy threshold is seen to scale as Re −2 for the two subcritical scenarios.This is in agreement with the theoretical study by Chapman [5] and DNSs carried out by Duguet et al. [4].More interestingly, the fitted curves for the variations of the critical energy thresholds with the Reynolds number obtained using either RNL simulations or DNS are nearly overlapping.Furthermore, the corresponding phase space trajectories are also well described by the RNL approximation for the two scenarios.In particular, the edge states found by Duguet et al. [4] are well recovered using a simplified set of equations.
We believe that such a confirmation is a preliminary step toward approximating the minimal seed in parallel shear flows using RNL models that are classically developed using extensive nonlinear optimizations based on full DNS [13,25].From such a perspective, the optimization procedure derived by Biau and Bottaro [19] could be an attractive prospect.This work is currently in progress.In particular, since the RNL system prevents any streamwise localization, this analysis should also indicate the role played by spanwise localization in the Re-scaling law proposed by Duguet et al. [12].
It should also be interesting to carry out optimizations as those done by Duguet et al. [4] using a larger subspace for the linear optimal mode basis and including the target times as an optimization parameter.Although within a DNS framework, such optimizations require a huge numerical effort, this could be more easily achieved using RNL simplification.In this regard, it would be of interest to evaluate the possibility of spanwise localization for the initial perturbation when the number of linear optimal modes is increased.
Finally, since subcritical paths leading to turbulence (OW and SV) are well described within the RNL framework, one may suggest using the RNL system to design controllers that can be utilized in full nonlinear dynamics.

Figure 1 .
Figure 1.Illustration of the bisection algorithm for Re = 1500.Integrated kinetic energy of the perturbation velocity versus time when the initial amplitude is varied.(a) Pure oblique wave scenario.(b) Streamwise vortex scenario.For case (b), the amplitude λ is associated with the minimum of E c (λ).

Figure 3 .
Figure 3. Re = 400.Threshold kinetic energy for the transition initiated by a pair of symmetric oblique waves as a function of streamwise length.The DNS results are extracted from Duguet et al. [4].

Figure 4 .
Figure 4. Re = 400.Edge state E 1 approached by an initial perturbation composed of a pair of symmetric oblique waves.Isocontours of streamwise velocity fields extracted at planes y = 0 (a) and y = 0.4 (c).Vectors for cross-stream components extracted at planes x = 3.2 (b) and x = 0.4 (d).