1. Introduction
Optimal perturbations represent a fundamental concept of the linear stability theory of shear flows and one highly cited textbook (Schmid & Henningson Reference Schmid and Henningson2001) centres much of the discussion on this concept and on the associated idea of transient growth of disturbances for systems ruled by non-normal linear operators. It has been thought for quite some time that, under nominally stable conditions, the most dangerous initial disturbances computed in a linear setting might undergo a substantial amplification, eventually capable of bringing the fluid system outside of the domain of attraction of the laminar solution (and also outside of the domain of validity of the linearized equations), so that transition can occur because of nonlinear interactions. This has provided the rationale for studying linear optimal disturbances via a number of approaches including eigenfunction expansion (Butler & Farrell Reference Butler and Farrell1992), expansion in periodic functions (Criminale et al. Reference Criminale, Jackson, Lasseigne and Joslin1997) and adjoint equations (Andersson, Berggren & Henningson Reference Andersson, Berggren and Henningson1999; Luchini Reference Luchini2000; Corbett & Bottaro Reference Corbett and Bottaro2001).
Linear, optimal perturbations in simple parallel shear flows, such as Couette, Poiseuille or Hagen–Poiseuille flows, take the form of streamwise or quasi-streamwise vortices (the streamwise wavenumber
${\it\alpha}$
either vanishes or is very small, in the latter case decreasing like
$Re^{-1}$
, with
$Re$
the Reynolds number), and their downstream evolution turns them into spanwise-periodic streaks, elongated along the main flow direction,
$x$
. While the presence of vortices and streaks seems appropriate in view of experimental observations of transitional flow features, their very early presence in the flow field is unwarranted since
$x$
-invariant flow structures trigger nonlinear interactions very slowly, and only when their amplitude is large, through the onset and growth of a travelling wave instability (Andersson et al.
Reference Andersson, Brandt, Bottaro and Henningson2001; Gustavsson Reference Gustavsson2009).
The traditional concept of optimal disturbances (those which give rise to the largest disturbance energy growth over a given time window, or over a given spatial development length, should one consider flows developing in space) is thus misleading if one searches for the initial waveform most suited to instigate transition, at least for the simplest shear flows. This has been understood for a few years now, and has produced a rising number of contributions on nonlinear optimals (Zuccher, Bottaro & Luchini Reference Zuccher, Bottaro and Luchini2006; Pringle & Kerswell Reference Pringle and Kerswell2010; Cherubini et al. Reference Cherubini, De Palma, Robinet and Bottaro2011; Monokrousos et al. Reference Monokrousos, Bottaro, Brandt, Di Vita and Henningson2011; Pringle, Willis & Kerswell Reference Pringle, Willis and Kerswell2012; Rabin, Caulfield & Kerswell Reference Rabin, Caulfield and Kerswell2012; Cherubini & De Palma Reference Cherubini and De Palma2013; Duguet et al. Reference Duguet, Monokrousos, Brandt and Henningson2013) which have highlighted the profound differences between linear and nonlinear transient growth. In the abstract of their contribution, Rabin et al. (Reference Rabin, Caulfield and Kerswell2012) state that all of their results ‘highlight the irrelevance of the linear energy gain optimal perturbation for predicting or describing the lowest-energy flow structure which triggers turbulence’. Nonlinear optimals are typically found via adjoint looping (Luchini & Bottaro Reference Luchini and Bottaro2014) and are very demanding in terms of computer memory and CPU time, because of the need to carry out several direct simulations of the Navier–Stokes equations and to store (or recompute, via a checkpointing procedure) the whole flow field at all times for use in the adjoint equations. Thus, a thorough parametric search (in terms of wavenumbers, disturbance amplitude and target time of the optimization) is still out of reach.
A weakly nonlinear alternative, much cheaper computationally than the full nonlinear optimization approach, is explored in the present paper. The approach aims at identifying initial perturbations which distort efficiently, via primary nonlinear interactions, the main flow. It will be shown that such optimal solutions take the form of oblique structures. It has been known since Schmid & Henningson (Reference Schmid and Henningson1992) that the early appearance of oblique waves elicits faster transition than either optimal streamwise vortices or Tollmien–Schlichting waves (see also Schmid & Henningson Reference Schmid and Henningson2001). The mechanism described by Schmid & Henningson (Reference Schmid and Henningson1992) has been called oblique transition; they have expressed the generic disturbance field,
${\it\phi}$
, in Fourier series, e.g.
$$\begin{eqnarray}\displaystyle {\it\phi}(x,y,z,t)=\mathop{\sum }_{m=-M}^{M}\mathop{\sum }_{n=-N}^{N}\tilde{{\it\phi}}_{m,n}(y,t)\,\text{e}^{\text{i}(m{\it\alpha}x+n{\it\beta}z)}, & & \displaystyle\end{eqnarray}$$
where
$x,y$
and
$z$
indicate, respectively, the streamwise, wall-normal and spanwise directions;
${\it\alpha}$
and
${\it\beta}$
are the streamwise and spanwise wavenumbers;
$M$
and
$N$
are positive integers that define the truncation of their direct numerical simulation (DNS). By imposing at
$t=0$
a disturbance field in the form of a pair of oblique waves with
$(m,n)=(1,\pm 1)$
, Schmid & Henningson (Reference Schmid and Henningson1992) have observed a rapid amplification of the mean flow distortion, mode
$(0,0)$
, accompanied by a growth (equally rapid but less vigorous) of the streamwise-independent component,
$(0,\pm 2)$
; conversely, the
$(\pm 2,0)$
mode remained an order of magnitude smaller. The
$(0,\pm 2)$
component is ruled by a non-normal evolution operator very sensitive to external forcing, and can thus be enhanced linearly by a transient growth mechanism, yielding large amplitude streaks via the lift-up effect (Landahl Reference Landahl1975). Much attention has been paid to the instability of such streaks, with lower emphasis placed on the fact that the streaks develop on top of a base flow which is being strongly distorted. The predominance of the
$(0,0)$
nonlinearly generated component, even above the
$(1,\pm 1)$
mode, is, however, evident in figures 1 and 2 of the paper by Schmid & Henningson (Reference Schmid and Henningson1992). Thus there is a clear interest in studying the optimal linear (modal or non-modal) growth of an oblique disturbance developing over a base flow which is being deformed nonlinearly by the primary wave itself. This calls for an autonomous cycle, schematized in figure 1, with the first generation of nonlinear interactions (the Reynolds stresses), which feeds onto the mean flow. The present paper describes the simplest possible model of such a cycle; the creation of harmonics past the
$(0,0)$
mode which concentrates in itself the majority of growth, is neglected. We thus study the rapid transient growth of an oblique wave and search for the optimal initial disturbance capable of creating and maintaining the self-sustained, neutral cycle sketched in figure 1. The appropriateness of the present minimal model will be demonstrated by comparison with DNS.
Figure 1. Self-sustaining cycle centred on the optimal amplification of oblique waves.
2. Formulation of the weakly nonlinear problem
A self-sustained cycle can be created by decomposing velocity and pressure into a steady laminar parallel state, a travelling wave and a slowly varying time-dependent base flow distortion. The idea for this decomposition stems from Biau & Bottaro (Reference Biau and Bottaro2009); however, other papers in the recent literature (Hall & Sherwin Reference Hall and Sherwin2010; Blackburn, Hall & Sherwin Reference Blackburn, Hall and Sherwin2013; Beaume et al. Reference Beaume, Chini, Julien and Knobloch2015) have also split the disturbance field into streamwise-invariant (i.e. mean) and streamwise-varying (i.e. fluctuating) components, reaching a set of equations in which the Reynolds stresses force the growth of the mean flow. The decomposition employed here assumes no a priori scaling with the Reynolds number and reads:
$$\begin{eqnarray}\left[\begin{array}{@{}c@{}}U(y)\\ 0\\ 0\\ P(x)\end{array}\right]+{\it\epsilon}\left[\begin{array}{@{}c@{}}\tilde{u} (x,y,z,t)\\ \tilde{v}(x,y,z,t)\\ \tilde{w}(x,y,z,t)\\ \tilde{p}(x,y,z,t)\end{array}\right]+{\it\epsilon}^{2}\left[\begin{array}{@{}c@{}}u_{00}(y,t)\\ v_{00}(y,t)\\ w_{00}(y,t)\\ p_{00}(y,t)\end{array}\right],\end{eqnarray}$$
where
${\it\epsilon}\in \mathbb{R}$
denotes the wave amplitude. In the following, the simplest wall-bounded flow is considered, but the procedure is general and can easily be applied to semi-infinite or infinite domains. The
$O({\it\epsilon})$
terms correspond to the travelling wave which we wish to optimize; those of
$O({\it\epsilon}^{2})$
are the mean flow distortion terms which dominate the nonlinear interactions at early times. Clearly, the decomposition proposed is insufficient at large times when higher modes are created and interact; it is however adequate, as we will show, to approximate the lower bound of the amplitude
${\it\epsilon}$
beyond which the flow can undergo transition to turbulence. The flow decomposition (2.1) is introduced into the Navier–Stokes equations and like-order terms are collected. At order
$O(1)$
we obtain
$$\begin{eqnarray}P_{x}=\frac{1}{Re}U_{yy},\end{eqnarray}$$
with subscripts denoting partial differentiation, and no-slip boundary conditions for
$U$
. The Reynolds number,
$Re$
, is based on the half-channel thickness and the wall velocity. For the case of the Couette flow treated here,
$P$
is constant and
$U=y$
, for
$y\in [-1,1]$
. The disturbance at
$O({\it\epsilon})$
is taken to be a simple monochromatic wave which travels along a direction inclined by
$\text{tan}^{-1}({\it\beta}/{\it\alpha})$
with respect to the
$x$
-axis:
$$\begin{eqnarray}(\tilde{\boldsymbol{u}},\tilde{p})(x,y,z,t)=(\boldsymbol{u}_{11},p_{11})(y,t)\,\text{e}^{\text{i}({\it\alpha}x+{\it\beta}z)}+({\boldsymbol{u}_{11}}^{\ast },p_{11}^{\ast })(y,t)\,\text{e}^{-\text{i}({\it\alpha}x+{\it\beta}z)},\end{eqnarray}$$
where superscript
$^{\ast }$
denotes complex conjugation,
$\tilde{\boldsymbol{u}}=(\tilde{u} ,\tilde{v},\tilde{w})$
and
${\it\alpha},{\it\beta}\in \mathbb{R}$
are the streamwise and spanwise wave numbers, respectively. Using the expression (2.3), the governing equations, linearized around the perturbed and time-varying mean flow
$[U(y)+{\it\epsilon}^{2}\,u_{00}(y,t),0,{\it\epsilon}^{2}\,w_{00}(y,t)]$
, are given by
$$\begin{eqnarray}\displaystyle & \displaystyle \text{i}{\it\alpha}u_{11}+{v_{11}}_{y}+\text{i}{\it\beta}w_{11}=0,\quad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {u_{11}}_{t}+\text{i}{\it\alpha}(U+{\it\epsilon}^{2}\,u_{00})u_{11}+v_{11}(U+{\it\epsilon}^{2}\,u_{00})_{y}+\text{i}{\it\beta}({\it\epsilon}^{2}\,w_{00})u_{11}+\text{i}{\it\alpha}p_{11}=\frac{1}{Re}{\it\Delta}_{k}u_{11},\quad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {v_{11}}_{t}+\text{i}{\it\alpha}(U+{\it\epsilon}^{2}\,u_{00})v_{11}+\text{i}{\it\beta}({\it\epsilon}^{2}\,w_{00})v_{11}+{p_{11}}_{y}=\frac{1}{Re}{\it\Delta}_{k}v_{11},\quad & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {w_{11}}_{t}+\text{i}{\it\alpha}(U+{\it\epsilon}^{2}\,u_{00})w_{11}+v_{11}({\it\epsilon}^{2}\,{w_{00}}_{y})+\text{i}{\it\beta}({\it\epsilon}^{2}\,w_{00})w_{11}+\text{i}{\it\beta}p_{11}=\frac{1}{Re}{\it\Delta}_{k}w_{11},\quad & \displaystyle\end{eqnarray}$$
where
${\it\Delta}_{k}=\partial ^{2}/\partial y^{2}-k^{2}$
and
$k^{2}={\it\alpha}^{2}+{\it\beta}^{2}$
. Equations (2.4)–(2.7) are accompanied by no-slip conditions on both the upper and lower walls and an initial condition
$\boldsymbol{u}_{11}(y,t=0)$
which we will optimize for. At order
$O({\it\epsilon}^{2})$
, the equations read
$$\begin{eqnarray}\displaystyle & \displaystyle v_{00}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {u_{00}}_{t}-\frac{1}{Re}{u_{00}}_{yy}=-[v_{11}{u_{11}^{\ast }}_{y}+\text{i}{\it\beta}w_{11}^{\ast }u_{11}+\text{c.c.}], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {p_{00}}_{y}=-[\text{i}{\it\alpha}u_{11}^{\ast }v_{11}+v_{11}{v_{11}^{\ast }}_{y}+\text{i}{\it\beta}w_{11}^{\ast }v_{11}+\text{c.c.}], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle {w_{00}}_{t}-\frac{1}{Re}{w_{00}}_{yy}=-[\text{i}{\it\alpha}u_{11}^{\ast }w_{11}+v_{11}{w_{11}^{\ast }}_{y}+\text{c.c.}], & \displaystyle\end{eqnarray}$$
plus the usual no-slip conditions on both walls. The key point of the analysis is the fact of retaining the
$x$
- and
$z$
-component of the velocity at
$O({\it\epsilon}^{2})$
in the expression of the base state above which the mode
$(1,1)$
is superposed (note also that we have anticipated the fact that
$v_{00}=0$
in writing (2.5)–(2.7)). The fact of having lifted a higher-order term to equations of lower order is motivated by the need to allow for a distorted base state. From the mathematical point of view, this corresponds to replacing the original series (2.1) in powers of
${\it\epsilon}$
into a sequence of unknown (and decreasing) functions of
${\it\epsilon}$
(Van Dyke Reference Van Dyke1975; Cousteix & Mauss Reference Cousteix and Mauss2007). The Orr–Sommerfeld stability equation for the Blasius boundary-layer flow is the best-known example of such an instance.
The system (2.4)–(2.7) plus (2.9) and (2.11) are solved by maximising the energy of the oblique mode at the target time
$t=T$
, defined as
$$\begin{eqnarray}e(T)=\frac{{\it\epsilon}^{2}}{2}\left.\int _{-1}^{1}(u_{11}u_{11}^{\ast }+v_{11}v_{11}^{\ast }+w_{11}w_{11}^{\ast })\,\text{d}y\right|_{t=T}.\end{eqnarray}$$
Alternative measures of growth could have been considered, for example based on the right-hand side terms of (2.9)–(2.11). In light of previous experience (e.g. Bottaro, Soueid & Galletti Reference Bottaro, Soueid and Galletti2006; Monokrousos et al.
Reference Monokrousos, Bottaro, Brandt, Di Vita and Henningson2011; Rabin et al.
Reference Rabin, Caulfield and Kerswell2012), the conclusions reached here would probably be mildly affected; this remains however to be verified. For assigned
${\it\epsilon}$
, the nonlinear problem is solved iteratively in the following manner:
-
(1) Maximise the energy
$e(T)$
, constrained by (2.4)–(2.7), over a given time span
$T$
, to find
$u_{11},v_{11},w_{11}$
. The optimization procedure is performed via adjoint looping; in the first iteration set
$u_{00}=w_{00}=0$
$\forall y,t$
. -
(2) Compute the right-hand side of (2.9) and (2.11), and solve for
$u_{00}(y,t)$
and
$w_{00}(y,t)$
under the initial condition
$u_{00}(y,0)=w_{00}(y,0)=0$
; (2.9) and (2.11) are heat equations with Reynolds stress terms acting as source. Then, go back to (1). -
(3) Convergence is declared when the final wave energy
$e(T)$
converges to within a defined precision. The normalization employed is
$e(0)={\it\epsilon}^{2}$
.
In order to quantify the results, suitable quantities are defined which can also be used in the limit
${\it\epsilon}\rightarrow 0$
to compare with linear optimal results from the literature. The gain
$G$
is
$$\begin{eqnarray}G({\it\epsilon},Re,{\it\alpha},{\it\beta},T)=\frac{e(T)}{e(0)},\end{eqnarray}$$
and is obtained from the converged solution of the procedure outlined above for given values of
${\it\alpha},{\it\beta},Re$
and
${\it\epsilon}$
in a defined temporal interval
$[0,T]$
. Finally, the mean-flow distortion is quantified by evaluating its ‘energy’:
$$\begin{eqnarray}E_{U}(t)=\frac{{\it\epsilon}^{4}}{4}\int _{-1}^{1}u_{00}^{2}(y,t)\,\text{d}y\quad \text{and}\quad E_{W}(t)=\frac{{\it\epsilon}^{4}}{4}\int _{-1}^{1}w_{00}^{2}(y,t)\,\text{d}y.\end{eqnarray}$$
2.1. Numerical procedure
The problem outlined at step (1) in the solution procedure is a constrained maximization problem which can be solved using Lagrange multipliers. The governing equations are written in primitive variable form and the resulting adjoint equations are derived using a discrete approach (see Luchini & Bottaro Reference Luchini and Bottaro2014). The spatial derivatives are discretized using second-order finite differences and a semi-implicit second-order scheme is used to advance in time. A uniform grid is used in the
$y$
-direction and 300 discrete points are sufficient to obtain a converged solution. The code has been tested on several cases found in the literature for
${\it\epsilon}\rightarrow 0$
(i.e. the mean flow is not distorted); in particular, the value of the optimal gain
$G=0.00118Re^{2}$
and the corresponding time at which it is achieved,
$T=0.117Re$
, with
${\it\alpha}=35/Re$
and
${\it\beta}=1.60$
, are recovered with less than 0.1 % error (Schmid & Henningson Reference Schmid and Henningson2001). Results are obtained imposing that convergence is reached when the relative difference in gain between two consecutive iterations is less than
$10^{-8}$
.
3. Weakly nonlinear results
To illustrate results, the case
$Re=400$
is chosen. We take four values of
${\it\epsilon}$
and fix the wavenumbers
${\it\alpha}$
and
${\it\beta}$
as well as the final time
$T$
; all parameters selected are given in the caption of figure 2, which displays the curves of gain (2.13) together with the energies of the mean flow modification for the cases in which
${\it\epsilon}>0$
. We choose a rather small value of the final time in order to avoid situations in which the nonlinear effects become dominant, for which an expansion of the waves including higher harmonics would be necessary.
Figure 2. (a) Gain
$G$
as a function of the final time
$t$
for
${\it\epsilon}\rightarrow 0$
,
${\it\epsilon}=0.015,0.0153,0.016$
, for
${\it\alpha}_{opt}=0.93$
,
${\it\beta}_{opt}=1.44$
(obtained by choosing
$T=20$
as target time for the optimization) and
$Re=400$
. (b) Energy of the mean flow modification,
$E_{U}$
(——) and
$E_{W}$
(– – –), for
${\it\epsilon}=0.015,0.0153,0.016$
. The arrows display increasing values of
${\it\epsilon}$
.
For
${\it\epsilon}<0.0153$
, all solutions are damped for
$t>15$
. For
${\it\epsilon}\approx 0.0153$
a limit condition is reached, with the different energy measures neither increasing nor decreasing. This corresponds to a threshold between initial disturbance states which relaminarize and states which might lead the flow to turbulence. When
${\it\epsilon}>0.0153$
, the measure of energy in figure 2 shows a monotonic increase. We can note that the two components of the energy of the mean flow distortion differ by almost one order of magnitude. This fact might lead to the idea that including the
${\it\epsilon}^{2}w_{00}$
terms in the wave equations, (2.5)–(2.7), yields a negligible effect. This will be further discussed in § 3.1.
Figure 3 shows, for
${\it\epsilon}=0.0153$
, the mean flow
$U$
and its distortion at
$t=0$
, 10, 15, 20, highlighting the presence of inflection points in both the streamwise and the spanwise velocity components.
Figure 3. Components of the mean flow
$U+{\it\epsilon}^{2}u_{00}$
(a) and
${\it\epsilon}^{2}w_{00}$
(b) at
$t=0,10,15,20$
for the case
$Re=400$
,
${\it\epsilon}=0.0153$
,
${\it\alpha}_{opt}=0.93$
and
${\it\beta}_{opt}=1.44$
. The arrows display increasing values of
$t$
.
The optimal disturbance wave at
$t=0$
is displayed in modulus and phase in figure 4, and through isocontours of
$u_{11}$
and
$(v_{11},w_{11})$
vectors in figure 5. The modulus of the optimal disturbance displays a maximum, both at
$t=0$
and at the target time, in correspondence with the inflection point of the mean flow. The initial perturbation is tilted in the
$(y,z)$
plane at
$t=0$
, and grows in time, developing a slant in the opposite sense under the forcing of the mean flow distortion shown in figure 3.
Figure 4. Modulus of velocity components and phase of the optimal oblique wave at
$t=0$
(a,b) and its output at the target time
$T=20$
(c,d). The
$u_{11}$
,
$v_{11}$
and
$w_{11}$
components of the velocity are presented in (a,c) as solid, dashed and dash-dotted lines, respectively. The phase, in radians, of the
$u_{11}$
,
$v_{11}$
and
$w_{11}$
components are shown in (b,d) as
$\circ$
,
$+$
and
$\ast$
, respectively, and the parameters are those of figure 3.
Figure 5. Vectors and contours, in the
$(y,z)$
plane, of the optimal oblique wave at
$t=0$
(a) and at the target time
$T=20$
(b), shown over two spanwise periods. The shaded contours represent the positive and negative streamwise velocity components, whereas the vectors represent wall-normal and spanwise components. Same parameters as in figure 3.
Figure 6 displays the gain
$G$
at
$T=20$
, in the Fourier plane, for four different values of
${\it\epsilon}$
. For
${\it\epsilon}\rightarrow 0$
, the results are those from the classical optimal perturbation theory and the maximum value of
$G$
is found for
$({\it\alpha},{\it\beta})=(0.34,1.72)$
. As the value of
${\it\epsilon}$
exceeds
${\approx}0.015$
, a second region emerges for larger values of
${\it\alpha}$
and smaller values of
${\it\beta}$
. This is more evident in figure 6(c) with
${\it\epsilon}=0.0153$
, while in figure 6(d), this new peak, with
$({\it\alpha},{\it\beta})=(0.93,1.44)$
, is clearly dominant. If
${\it\epsilon}$
is further increased (not shown here) the values of
${\it\alpha}$
and
${\it\beta}$
, where the peak of
$G$
is found, remain at 0.93 and 1.44, respectively. We denote these weakly nonlinear optimal values of the wavenumbers as
${\it\alpha}_{opt}$
and
${\it\beta}_{opt}$
.
In figure 7(a) we report
${\it\alpha}_{opt}$
and
${\it\beta}_{opt}$
for different values of the final time
$T$
used in the optimization procedure. For
$T<15$
, the streamwise wavenumber is larger than the spanwise wavenumber, in particular at
$T=10$
${\it\alpha}\approx 2\,{\it\beta}$
. For
$T>15$
the opposite is true and both
${\it\alpha}_{opt}$
and
${\it\beta}_{opt}$
tend asymptotically to a constant value for large
$T$
. In particular, the asymptotic value of the spanwise wavenumber approaches its linear counterpart, i.e. 1.60. The value of
${\it\epsilon}$
which indicates a limit condition is here denoted
${\it\epsilon}_{crit}$
and its value as a function of the target time
$T$
of the optimization is shown in figure 7(b). As one would expect,
${\it\epsilon}_{crit}$
decreases with
$T$
, reaching the asymptotic value
${\it\epsilon}_{crit}\approx 0.015$
already for
$T$
as low as
$30$
.
Figure 6. Gain
$G$
in the
${\it\alpha}$
–
${\it\beta}$
plane for
$Re=400$
and
$T=20$
; (a)
${\it\epsilon}\rightarrow 0$
with contour values from 20 to 120, (b)
${\it\epsilon}=0.0145$
with contour values from 20 to 240, (c)
${\it\epsilon}=0.0153$
with contour values from 20 to 300, (d)
${\it\epsilon}=0.01603$
with contour values from 20 to 380. The interval among adjacent isolines is
${\rm\Delta}G=20$
in frames. The maximum value of
$G$
for each
${\it\epsilon}$
is denoted by a filled circle.
Figure 7. (a) Optimal streamwise (
${\it\alpha}_{opt}$
) and spanwise (
${\it\beta}_{opt}$
) wavenumbers as a function of
$T$
; (b) critical amplitude (
${\it\epsilon}_{crit}$
) as a function of
$T$
.
$Re=400$
.
3.1. Neglecting the correction to the transverse velocity
In figure 2 is was shown that the measure
$E_{W}$
was approximately one order of magnitude smaller than
$E_{U}$
. This seems to imply that the mean-flow component,
$w_{00}$
, can be safely omitted from the equations for the optimal oblique disturbance; such an approximation has been made, for example, by Farrell & Ioannou (Reference Farrell and Ioannou2012). This is, however, not true, at least in the present model, as the example here will show. Let us consider the case in which
$Re=400$
and the terms
${\it\epsilon}^{2}w_{00}$
are not included in (2.5)–(2.7). The corresponding optimal wavenumbers, optimized also with respect to
$T$
, are
${\it\alpha}_{opt}=0.27$
and
${\it\beta}_{opt}=1.083$
. In figure 8, we report the values of the gain
$G$
and energies
$E_{U}$
and
$E_{W}$
for three different values of
${\it\epsilon}$
. A limiting condition is found for
${\it\epsilon}\approx 0.0084$
, similar to what shown in figure 2. Here, the value of
$w_{00}$
is not included in the governing equations and the value of
$E_{W}$
is only shown for comparison with the value of
$E_{U}$
. Clearly,
$E_{W}$
is two orders of magnitude smaller than
$E_{U}$
. If we now include the
${\it\epsilon}^{2}w_{00}$
terms in (2.5)–(2.7) and compute the gain and energies due to the mean-flow distortion using the values of
${\it\alpha}_{opt}$
and
${\it\beta}_{opt}$
previously obtained, no limiting condition is found, and all solutions are damped for
$t>40$
(cf. figure 8
a, thick solid lines). This holds for all values of
${\it\epsilon}$
used and also for larger values of
${\it\epsilon}$
(not shown here).
The result is somewhat unexpected, and it confirms that acting on the spanwise mean flow is an efficient strategy to control fluid motion (Gatti & Quadrio Reference Gatti and Quadrio2013).
Figure 8. (a) Gain
$G$
as a function of time
$t$
for
${\it\epsilon}=0.0078,0.0084,0.0086$
,
${\it\alpha}_{opt}=0.27$
and
${\it\beta}_{opt}=1.083$
(obtained in the case excluding the
${\it\epsilon}^{2}w_{00}$
terms) and
$Re=400$
. Lines (——) denote the case including
${\it\epsilon}^{2}w_{00}$
-terms, (– – –) is the case without the
${\it\epsilon}^{2}w_{00}$
terms; the arrow displays increasing values of
${\it\epsilon}$
. (b) Energy of the mean flow modification when
${\it\epsilon}=0.0084$
for both
$E_{U}$
and
$E_{W}$
in the case (– – –) without
${\it\epsilon}^{2}w_{00}$
terms and (——) including
${\it\epsilon}^{2}w_{00}$
terms.
4. Nonlinear results and comparison
Several direct numerical simulations have been performed by initializing the perturbation field using the weakly nonlinear optimal solutions obtained for
$T=20$
,
${\it\alpha}_{opt}=0.93$
and
${\it\beta}_{opt}=1.44$
, and different values of
${\it\epsilon}$
ranging from
$0$
to
$0.016$
. The computational domain has length
$2{\rm\pi}/{\it\alpha}_{opt}$
in the streamwise direction and
$2{\rm\pi}/{\it\beta}_{opt}$
in the spanwise domain. The chosen domain is discretized using a staggered grid and the Navier–Stokes equations are solved by a fractional-step method with second-order accuracy in space and time (Verzicco & Orlandi Reference Verzicco and Orlandi1996). In order to provide a measure of the perturbation energy in the whole computational domain, we define two integral energies,
$E$
and
$E_{wavy}$
, the former based on the whole perturbation,
$\boldsymbol{u}$
, and the latter only on its wavy part,
$\boldsymbol{u}_{wavy}=\boldsymbol{u}-{\it\epsilon}^{2}\boldsymbol{u}_{00}$
. The perturbation and wavy integral energies are thus defined as:
$$\begin{eqnarray}E=\frac{1}{2V}\,\int _{V}u^{2}+v^{2}+w^{2}\,\text{d}V,\quad E_{wavy}=\frac{1}{2V}\int _{V}u_{wavy}^{2}+v_{wavy}^{2}+w_{wavy}^{2}\,\text{d}V,\end{eqnarray}$$
with
$V=(8{\rm\pi}^{2})/({\it\alpha}_{opt}{\it\beta}_{opt})$
the computational domain. As in § 3, gains are defined with respect to the energy at the initial time. In particular, the gain
$G$
predicted by the weakly nonlinear model will be compared with the wavy energy gain,
$E_{wavy}(t)/E_{wavy}(0)$
provided by the DNS, which, for brevity, will also be labelled
$G$
. Figure 9(a) shows the wavy energy gain obtained by DNS using the initial optimal solution for
${\it\epsilon}=0.014$
, compared with that predicted by the weakly nonlinear model for the same value of
${\it\epsilon}$
. The two curves show a remarkable agreement, indicating that the weakly nonlinear expansion succeeds in predicting the evolution of the perturbations, provided that the nonlinear effects are not too strong. However, for considerably larger values of
${\it\epsilon}$
(in particular, for
${\it\epsilon}>{\it\epsilon}_{crit}$
), the prediction provided by the weakly nonlinear model for
${\it\epsilon}=0.016$
, after a finite time (equal to approximately 10 time units in figure 9
b), begins to differ from the DNS. This can be observed in figure 9(b) by comparing the red dotted and the solid black curves. This discrepancy is clearly due to the fact that, for large enough amplitudes, nonlinear effects become too strong to be correctly approximated by a weakly nonlinear expansion. However, despite the differences in energy amplification, the perturbation resulting from the weakly nonlinear model is able to induce a self-sustaining cycle leading to turbulence. The dashed line in figure 9(b) shows that, when the field obtained at
$t=12$
by the model is injected as an initial perturbation into the DNS, sudden transition to turbulence is observed. A similar behaviour is observed when using optimal solutions extracted at larger times, provided that the initial value of
${\it\epsilon}$
is sufficiently large.
To determine the critical value of
${\it\epsilon}$
for transition to turbulence, simulations are performed initializing the perturbation field with several optimal solutions extracted at
$t=15$
; the initial conditions providing the
$t=15$
solutions are the weakly nonlinear optimal disturbances found at different values of
${\it\epsilon}$
for the same target time,
$T=20$
. This choice is motivated by the fact that the energy gain predicted by the model reaches a maximum at approximately
$t\approx 15$
, as shown, e.g. in figure 9.
Figure 9. (a) Energy gain computed by the weakly nonlinear model (dashed line) and by the DNS (solid line) both initialized with the optimal solution obtained at
$t=0$
for
$T=20$
,
${\it\epsilon}=0.014$
. (b) Energy gain versus time obtained by DNS (black solid line) and by the weakly nonlinear model (red dotted line) initialized with the initial optimal solutions for
${\it\epsilon}=0.016$
. The black dashed line represent the wavy energy gain for a DNS initialized by the optimal solution extracted at
$t=12$
for the same value of
${\it\epsilon}$
.
Figure 10. Disturbance energy
$E$
versus time measured by a DNS initialized with the weakly nonlinear optimal states at
$t=15$
,
${\it\epsilon}=0.0151$
, 0.0152, 0.0153, 0.0154 (increasing from bottom to top).
Figure 11. Variation in time of the amplitude of the different modes extracted by Fourier transform of the velocity snapshots obtained by DNS initialized with the weakly nonlinear solutions at
$t=15$
,
${\it\epsilon}=0.0151$
(a) and
${\it\epsilon}=0.0152$
(b). The legend indicates the streamwise and spanwise wavenumbers, as multiples of the
${\it\alpha}$
and
${\it\beta}$
values chosen for the optimization.
Figure 10 shows the integral energy curves obtained by DNS for four values of
${\it\epsilon}$
ranging from
$0.0151$
to
$0.0154$
. For the lowest value of
${\it\epsilon}$
, the perturbation energy decreases in time, leading the flow to relaminarization. For
${\it\epsilon}\geqslant 0.0152$
, the energy increases, reaching a statistically steady value, while the flow is characterized by a chaotic motion. This critical value of
${\it\epsilon}$
for transition to turbulence, determined by DNS, is very close to the critical value predicted by the model,
${\it\epsilon}_{crit}=0.0153$
. This result might have been anticipated, considering that the strong increase of the energy gain observed by the model for
${\it\epsilon}>{\it\epsilon}_{crit}$
is associated with the development of strong nonlinear effects. In order to investigate the mechanisms leading the optimal perturbation towards a sudden self-sustained state, we have performed a spatial Fourier transform of the streamwise disturbance velocity on the
$y=-0.4$
plane, at each
$2.5$
time units (by selecting a different wall-normal position or variable the results change only marginally). We have tracked the development in time of modes having different wavenumbers, using their amplitudes
$A_{m,n}=\sqrt{\text{Re}(\tilde{u} _{m,n})^{2}+\text{Im}(\tilde{u} _{m,n})^{2}}$
. The Fourier amplitudes of the most relevant modes, normalized with respect to the initial amplitude of the mode
$A_{0,0}$
, are provided in figure 11 for
${\it\epsilon}=0.0151$
and
$0.0152$
. The mode
$(1,1)$
, corresponding to the oblique wave, is represented by a dashed red line, the mean flow deformation
$(0,0)$
is given by the solid line and streaky structures,
$(0,2)$
and
$(0,1)$
by the blue long-dashed and the green dashed-dotted lines, respectively. For the lower value of
${\it\epsilon}$
, the scenario appears similar to the early phases of oblique transition. At the initial time of the DNS (
$t=15$
), the flow is composed of the wavy disturbance and its corresponding mean flow deformation. While the energy of the former decreases in time, the latter increases and saturates towards an almost steady state; at the same time, streaky modes begin to be generated from the initial numerical noise, and slowly saturate in time. However, the amplitude of the induced mean flow distortion is not sufficiently high to sustain a regeneration cycle, and after a long saturation phase, the streaks decay towards the laminar state (figure 11
a). However, when
${\it\epsilon}$
is slightly larger, the self-sustained cycle is triggered very suddenly, roughly after
$10$
time units. As one can observe in figure 11(b), the streaky modes generated by the initial noise do not experience the slow phase of saturation and secondary instability observed in the classical bypass transition scenario, but they rapidly reach a chaotic dynamics. This behaviour is not surprising since the perturbation has been optimized by taking into account nonlinear effects; thus, it should be able to skip some of the linear phases of growth of the perturbations, similarly to fully nonlinear optimals (cf. for instance, Cherubini et al.
Reference Cherubini, De Palma, Robinet and Bottaro2011).
Figure 12. Isosurfaces of the streamwise perturbation to the base flow (yellow and blue for
$u=\pm 0.3$
, respectively), extracted from a DNS initialized with the optimal solution at
$t=15$
for
${\it\epsilon}=0.0152$
. The snapshots are extracted at
$t=25~(a),30~(b),35~(c),45~(d)$
.
Figure 13. Isosurfaces of the streamwise perturbation to the base flow (yellow and blue for
$u=\pm 0.2$
, respectively), extracted from a DNS initialized with the optimal solution at
$t=15$
for
${\it\epsilon}=0.0151$
. The snapshots are extracted at
$t=25$
(a), 30 (b), 45 (c), 75 (d).
The different phases of the weakly nonlinear optimal transition are displayed in figure 12 where isosurfaces of the streamwise perturbation are plotted together with velocity vectors in a plane parallel to the direction of the oblique wave. At
$t=25$
(a), the mean flow distortion has induced an inflectional region at
$y\approx 0$
. The inflectional instability triggered in this region produces an array of vorticity rolls normal to the direction of the oblique wave (b). The rolls increase in amplitude, transporting the base flow momentum, creating inclined zones of strong high and slow streamwise velocity (c). It is noteworthy that the presence of inclined rolls and low- and high-momentum regions is a key feature found in fully nonlinear optimal perturbations in Couette flow (see Monokrousos et al.
Reference Monokrousos, Bottaro, Brandt, Di Vita and Henningson2011; Rabin et al.
Reference Rabin, Caulfield and Kerswell2012; Cherubini & De Palma Reference Cherubini and De Palma2013). Finally, after
$55$
time units, these chaotic structures lose memory of the initial inclination provided by the vortices, and a self-sustained chaotic motion is observed. However, for a lower value of
${\it\epsilon}$
, although similar structures are created at early times, they are not able to self-sustain. Figure 13 shows that after
$30$
time units, rolls of vorticity normal to the direction of the oblique wave begin to be generated. However, their amplitude being too weak, they cannot be sustained in time by nonlinear effects; thus, they begin to decay. During this slow decay phase, the streamwise vorticity component of these rolls slowly begin to generate streamwise streaks, which can be observed in figure 13(d). However, since the regeneration cycle is not sustained due to the weak effect of nonlinearities, the streaks saturate and then slowly decay.
The trajectory of the flow in these two cases can be observed in a projection of the phase space given by the energy input rate,
$I$
, and the integral dissipation rate,
$D$
, normalized by their respective laminar values and defined as:
$$\begin{eqnarray}\displaystyle & \displaystyle I=\frac{{\it\alpha}_{opt}{\it\beta}_{opt}}{4{\rm\pi}^{2}}\int _{0}^{2{\rm\pi}/{\it\alpha}_{opt}}\int _{0}^{2{\rm\pi}/{\it\beta}_{opt}}\frac{1}{2}\left[\left(\frac{\partial u}{\partial y}\right)_{y=-1}+\left(\frac{\partial u}{\partial y}\right)_{y=1}\right]\,\text{d}x\,\text{d}z, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle D=\frac{1}{V}\int _{0}^{2{\rm\pi}/{\it\alpha}_{opt}}\int _{-1}^{+1}\int _{0}^{2{\rm\pi}/{\it\beta}_{opt}}|{\bf\omega}|^{2}\,\text{d}x\,\text{d}y\,\text{d}z, & \displaystyle\end{eqnarray}$$
where
${\bf\omega}$
represents the vorticity vector. This projection of the phase space (see, e.g. Kawahara & Kida Reference Kawahara and Kida2001) allows visualization of the balance between the energy injected in the system through the moving walls, and that dissipated by viscosity in the whole domain. As discussed by Kawahara & Kida (Reference Kawahara and Kida2001), the laminar attractor sits in the point
$I=D=1$
; the turbulent statistically stationary state is close to
$I=D\approx 3$
, whereas all other equilibrium solutions are situated on the
$I=D$
diagonal, indicated in the two panels of figure 14 by dashed lines. As shown in figure 14(a), the trajectories start from the upper left part of the plane (the empty triangles indicating the starting point of the trajectory and its direction), corresponding to an energy input very close to the laminar one and a rather high dissipation value. This is not surprising since the initial flow structures are not placed close to the wall, where energy enters the system, but rather around the centre of the domain in
$y=0$
. The trajectory shown by the thinner line indicates that, for
${\it\epsilon}=0.0151$
, after a rapid decrease of the dissipation, the trajectory bends towards the
$I=D$
bisector, locus of equilibrium solutions. However, the perturbation amplitudes reached are not sufficiently high to escape from the laminar attractor. However, for
${\it\epsilon}=0.0153$
, the trajectory spends much time close to the bisector, and then leaves towards higher values of the energy input and dissipation rate, which are characteristic of the turbulent dynamics. The slow-down of the trajectory in correspondence with the bisector, in a region characterized by
$I$
and
$D$
values intermediate between the laminar and the turbulent ones, might be the signature of the passage of the perturbation close to a relative attractor on the laminar–turbulent boundary (the edge state, see Schneider & Eckhardt Reference Schneider and Eckhardt2006; Skufca, Yorke & Eckhardt Reference Skufca, Yorke and Eckhardt2006; Schneider et al.
Reference Schneider, Gibson, Lagha, De Lillo and Eckhardt2008). Indeed, for an initial amplitude far from the laminar–turbulent boundary,
${\it\epsilon}=0.0157$
(green line), the trajectory does not pass close to the bisector before reaching turbulence, but displays strong values of
$D$
typical of a bursting phenomenon (see also van Veen & Kawahara Reference van Veen and Kawahara2011; Cherubini & De Palma Reference Cherubini and De Palma2013).
Finally, in order to evaluate the efficiency of the weakly nonlinear optimal perturbations in inducing turbulence, we have compared these trajectories with those followed by a linear optimal perturbation, computed for the same target time. The linear optimal perturbation has been rescaled in amplitude in order to have, at
$t=15$
, the same integral energy of the weakly nonlinear solutions for
${\it\epsilon}=0.0151$
and
$0.0157$
. The two thinner curves in figure 14(b) show that, for these amplitudes, the linear optimal perturbation is not able to induce transition, but rather leads to very fast relaminarization. Increasing the initial energy of these perturbations by
$5$
times still does not lead the flow to transition (as shown by the thicker line in figure 14
b), but rather induces a strong bursting event which slowly decays in time. Thus, we can conclude that weakly nonlinear optimal perturbations are indeed more effective in inducing transition than linear ones.
Figure 14. Trajectories on a
$I$
–
$D$
projection of the phase space of the (a) weakly nonlinear optimal solutions obtained for
${\it\epsilon}=0.0151,0.0153,0.0157$
and (b) of the corresponding linear optimal solutions rescaled with increasing amplitudes. The starting points of the trajectories, as well as their direction, are indicated by the white triangles, whereas the full triangles indicate the ending points. The value of
${\it\epsilon}$
is increasing with the thickness of the curves.
5. Conclusions
A formulation of the optimal-perturbation problem in shear flows is proposed, which takes into account the distortion of the base flow by the action of the Reynolds stresses, modelled by averaging the products of single-mode disturbance waves. The technique proposed permits consideration of the effect of the amplitude
${\it\epsilon}$
of the initial disturbance field, and isolation of a strongly receptive region of wavenumber space, for
$T$
sufficiently small for the theory to remain tenable, which differs significantly from the most selective range of linear optimal disturbances.
The simplicity of the weakly nonlinear model permits computation of the thresholds separating the initial states which relaminarize from those which grow without bound, in a fraction of the time compared to full DNS and thereby the model represents a valuable tool for parametric studies of the early stages of the laminar–turbulent transition. The importance of including the spanwise mean-flow distortion terms in the equations describing the transient amplification of oblique flow structures has been demonstrated; although the
${\it\epsilon}^{2}w_{00}$
terms are small in magnitude, their influence on the cross-stream terms in the wave equations is not negligible.
It would, however, be illusory to expect for the present highly truncated model to yield accurate transition thresholds, since transition relies on the growth and interaction of several modes. We have thus performed DNS to unravel the transition scenarios induced by the weakly nonlinear optimal disturbances. Concerning the energy amplification, a remarkable agreement is found between the DNS and the weakly nonlinear model results, for initial amplitudes below threshold. For larger amplitudes, nonlinear effects become too strong to be correctly captured by a weakly nonlinear expansion. However, when
${\it\epsilon}\geqslant {\it\epsilon}_{crit}$
, the optimal perturbation resulting at finite time from the weakly nonlinear model is able to induce a self-sustaining cycle leading to turbulence, whereas a linear optimal disturbance computed for the same target time leads to relaminarization, even when the initial energy is rescaled to very large values. The self-sustained cycle identified is triggered by an inflectional instability at the centre of the domain, resulting from the mean flow distortion induced by Reynolds stresses. This inflectional instability generates an array of vorticity rolls in the direction normal to the oblique wave, which, while increasing in amplitude, transport the base flow momentum, creating inclined zones of high and low streamwise velocity. Finally, these chaotic structures lose memory of the initial inclination provided by the vortices, and sustained chaotic motion is observed.
Acknowledgement
The financial support to the PRIN 2012 project, no. D38C1300061000, from the Italian Ministry of Education is gratefully acknowledged.
