# Turbulent transition in a truncated one-dimensional model for shear flow

J. H. P. Dawes, W. J. Giles

## Abstract

We present a reduced model for the transition to turbulence in shear flow that is simple enough to admit a thorough numerical investigation, while allowing spatio-temporal dynamics that are substantially more complex than those allowed in previous modal truncations. Our model allows a comparison of the dynamics resulting from initial perturbations that are localized in the spanwise direction with those resulting from sinusoidal perturbations. For spanwise-localized initial conditions, the subcritical transition to a ‘turbulent’ state (i) takes place more abruptly, with a boundary between laminar and turbulent flows that appears to be much less ‘structured’ and (ii) results in a spatio-temporally chaotic regime within which the lifetimes of spatio-temporally complicated transients are longer, and are even more sensitive to initial conditions. The minimum initial energy E0 required for a spanwise-localized initial perturbation to excite a chaotic transient has a power-law scaling with the Reynolds number E0Rep with p≈−4.3. The exponent p depends only weakly on the width of the localized perturbation and is lower than that commonly observed in previous low-dimensional models where typically p≈−2. The distributions of lifetimes of chaotic transients at the fixed Reynolds number are found to be consistent with exponential distributions.

## 1. Introduction

The transition from laminar to turbulent states has been a central problem in fluid mechanics for many decades. Since the 1960s, promising lines of attack have been opened up through the use of ideas from dynamical systems theory coupled to the thorough investigation of reduced versions of the Navier–Stokes equations. Perhaps the best-known example of such a reduction is the derivation of the Lorenz equations as a model of the onset of thermal convection in a layer of fluid heated from below (Lorenz 1963). Although the set of three nonlinear ordinary differential equations (ODEs) that comprise the (Lorenz 1963) model displays a wealth of interesting dynamical behaviour, little of this is relevant to the original fluid-mechanical problem. However, similar approaches yield convincing agreement over wide ranges of parameter values in other fluid-mechanical situations, for example, thermal convection in the presence of a magnetic field, and the onset of Taylor vortices in the flow between rotating coaxial concentric cylinders. In situations such as these, the flow undergoes a series of bifurcations before becoming chaotic, in a sense that can be given a clear meaning in terms of the behaviour of the reduced model comprising nonlinear ODEs.

In contrast, the transition to turbulence in shear flows appears not to proceed through a sequence of bifurcations, but to be linked to the appearance, at a critical Reynolds number Rec, of a chaotic saddle in the phase space: a complicated collection of unstable equilibria and time-periodic orbits which results in ever longer transients before the flow relaxes to a purely laminar state. For a review, and substantial numbers of references to the literature (although with an emphasis on the pipe flow), see Kerswell (2005). The significance of Rec is that for Re<Rec, the laminar state is a global attractor and trajectories appear to evolve rapidly towards it, while for Re>Rec other (usually unstable) invariant sets exist in the phase space. These new invariant sets cause increasingly long transient excursions to take place before relaminarization occurs, for initial conditions that are sufficiently far from the laminar profile. Small perturbations to the laminar state still decay rapidly towards it, and there appears to be a distinct boundary separating the behaviour of trajectories into those that relax rapidly and those which undergo long transient excursions. This laminar-turbulent boundary is sometimes referred to as the ‘edge of chaos’. Despite these differences in phenomenology and dynamics, reduced models constructed along the same lines as the Lorenz model have provided substantial insight into both the physical origin of this transition to self-sustaining complicated motion, and the mathematical organization of equilibria and periodic orbits inside the chaotic saddle.

The study by Waleffe (1997) provides the direct inspiration for the present work. Waleffe showed that a low-order reduced model could be constructed that elucidated the different elements of a self-sustaining process (SSP) that allowed sufficiently large deviations from the laminar flow profile to persist indefinitely. The SSP can be briefly described as follows. Weak streamwise vortices (i.e. vortical rolls whose axes are aligned with the primary flow direction) distort the streamwise velocity profile by moving high- and low-speed fluid around. This distortion generates streamwise streaks of fluid that are moving faster and slower than the fluid around them. The streaks are unstable to modes which create vortical eddies oriented in the wall-normal direction, orthogonal to the streamwise vortices, and the resulting three-dimensional re-organization of the flow reinforces the streamwise vortices.

Waleffe's model simplified the Navier–Stokes equations first by considering, not plane Couette flow between rigid boundaries, but a modified problem in which the boundaries are stress-free and the laminar profile is sinusoidal, sustained by an artificially applied pressure term. It appears that the physics of the SSP is rather insensitive to these modifications. The second set of simplifications concern the Galerkin expansion of the velocity field in all three directions: streamwise (x), wall normal (y) and spanwise (z). Thus, Waleffe considered a model comprising eight of the lowest-wavenumber Fourier modes in this Galerkin expansion. These eight modes were chosen in order to capture the central elements of the SSP, and to be self-consistent in the sense that the nonlinear (quadratic) interactions between these eight modes preserve energy just as the full advective nonlinearity in the Navier–Stokes equations does. Having projected out the spatial dependence of the dynamics onto these modes, the problem reduces to a far simpler set of ODEs for the time-dependent mode amplitudes. Subsequent work, in particular by Eckhardt & Mersmann (1999) and Moehlis et al. (2004, 2005), extended Waleffe's model to include an additional physical insight: that the basic laminar profile of the shear flow will itself be modified by the nonlinear interactions with streamwise vortices and streaks. Moehlis et al. developed a nine-mode ODE model that was amenable to investigation in substantial detail. In particular, they discussed the lifetimes of perturbations as the Reynolds number increased, locating the onset of complicated dynamics, and they also discussed the probabilistic distribution of lifetimes at the fixed Reynolds number from randomized initial conditions of equal energy.

It is clear that the assumptions made by these authors concerning spatial periodicity seems much easier to argue for in the streamwise and wall-normal directions than in the spanwise direction; indeed, recent work by Schneider et al. (2009, 2010a,b) (see also Duguet et al. 2009) has concentrated in understanding the formation of structures which are spatially localized in z: very far from periodic in this direction.

Bearing this in mind, we propose, in this paper, an extension of Waleffe's model which is a reduction of the Navier–Stokes equations to a collection of partial differential equations (PDEs) in z and t: we adopt the Fourier mode truncations that Waleffe used in order to remove the dependence on the streamwise x and wall-normal y coordinates since numerical work shows that, at least for some of these equilibrium and periodic orbit states, the flow structure in these coordinates can be well approximated by a small number of Fourier modes.

By retaining full dependence on the spanwise coordinate z we admit both the spatially periodic solutions of Waleffe and the formation of localized states. In addition, a wealth of spatio-temporal complexity is allowed. Our study is very similar in spirit to work by Manneville & Locher (2000), Manneville (2004) and Lagha & Manneville (2007), who preserve full resolution in two directions (streamwise and spanwise) and use the Galerkin truncation only in the third (‘wall-normal’) direction. This enables these authors to consider the dynamics around turbulent spots that are localized in both the streamwise and spanwise directions, at the cost of more intensive numerical computations. Their work is therefore complementary to that presented here. We also remark that reduced models, of different kinds, have been used both in pipe flow (Willis & Kerswell 2009) and in understanding the formation of turbulent-laminar bands in plane Couette flow (Barkley & Tuckerman 2007).

The structure of the paper is as follows. In §2 we discuss the derivation of our PDE extension of Waleffe's ODE model. In §3 we present the results of our numerical investigations into the transition to turbulence described by the model. We conclude in §4.

## 2. Derivation of the PDE model

In this section we define sinusoidal shear flow and summarize the Galerkin truncation that we use to derive our simplified model for turbulent transition.

Following Waleffe (1997) and Moehlis et al. (2004), we use the usual Cartesian coordinate conventions that x is the downstream direction (‘streamwise’), y is the direction of the shear gradient, i.e. normal to the sidewalls and z is the spanwise direction. We write the Navier–Stokes equations for the incompressible flow in the non-dimensionalized form: 2.1where we have scaled lengths by h/2, where h is the width of the channel, velocities by U0 the velocity of the laminar profile at a distance h/4 from the upper boundary and pressure by where ρ is the density of the fluid. The Reynolds number Re is therefore given by Re=U0h/(2ν), ν being the kinematic viscosity. The evolution of equation (2.1) is subject to the usual incompressibility condition ∇⋅u=0 and the conditions for impermeable and stress-free upper and lower boundaries 2.2The flow is assumed to be periodic in the x- and z-directions, with periodicities Lx and Lz, respectively. We take the body force term F(y) to be 2.3where ex is the unit vector in the x-direction and it is convenient to define the constant β=π/2. The laminar profile is then the steady solution 2.4of equations (2.1)–(2.3), as shown in figure 1. We note that although the ‘sinusoidal shear flow’ profile (2.4) has an inflection point, the flow is linearly stable for all Re (Drazin & Reid 1981).

Figure 1.

Geometry of the domain, illustrating the basic sinusoidal shear profile which is sustained by the applied body force term.

We now turn to our solution ansatz. We write the velocity field in the form 2.5where the subscripts denote mean, toroidal and poloidal components which are, respectively, expressed as sums of the first few Fourier modes in each case: and We define the streamwise and spanwise wavenumbers α=2π/Lx and γ=2π/Lz for notational convenience. The amplitudes A1,…,A8 are functions of z and t whose evolution can be obtained by substituting equation (2.5) into equation (2.1). Note that the form of the ansatz implies that incompressibility and the boundary conditions (2.2) are automatically satisfied.

The amplitudes A1,…,A8 correspond exactly, in terms of their Fourier dependence in x and y, to the modes selected by Waleffe (1997). For reference, table 1 summarizes the correspondence.

View this table:
Table 1.

Comparison of notation.

We now briefly describe the role that each mode plays in the dynamics. A1 is the amplitude of the sinusoidal shear profile; the laminar state corresponds to , A2=⋯=A8=0. A2 describes variations in z of the streamwise velocity, i.e. the formation of streamwise streaks. A3 describes the formation of x-independent streamwise vortices that redistribute the shear profile. Modes A4,…,A7 describe, as in Waleffe (1997), the development of x-dependent distortions of the streaks, and, in particular, the linear instability of the x-independent streaks described by A2. These modes have no velocity component in the vertical (i.e. y) direction. A8 describes ‘oblique rolls’ and, in contrast to modes A4,…,A7, has a non-zero vertical velocity component but no vertical vorticity component.

To derive evolution equations for the modes A1,…,A8, we use Fourier orthogonality in the x- and y-directions combined with projections onto individual components of either the velocity field u, the vorticity field ω=∇×u or (for A8) the curl of the vorticity field. For later convenience, we introduce the differential operators , and which correspond to the action of −∇2 on different Fourier modes: We denote ∂Aj/∂z by Aj for j=1,…,8. Considering first the ex component of u in equation (2.1), we obtain the following PDEs in z and t for A1 and A2: 2.6and 2.7Now we turn to the vorticity ω. For A3, A4 and A6 we find that only one component of ω contains a contribution from each of these; taking the curl of equation (2.5) the terms involving A3, A4 and A6 are Therefore, it is straightforward to consider the x and y components of the vorticity equation obtained by applying the operators ex⋅∇× and ey⋅∇× to equation (2.1) in order to obtain evolution equations for A3, A4 and A6: 2.8 2.9 and 2.10

Evolution equations for A5, A7 and A8 are derived similarly, using (for A5 and A7) the projection operator ey⋅∇× and (for A8) the projection operator ey⋅∇×∇× applied to equation (2.1). The resulting evolution equations are 2.11 2.12 and 2.13 Equations (2.6)–(2.13) are a closed set of nonlinear PDEs in z and t which form a truncated model of the dynamics of sinusoidal shear flow. Crucially these equations satisfy two consistency checks: the nonlinear terms in equations (2.6)–(2.13) conserve energy and so reflect the conservative nature of the full u⋅∇u nonlinearity in equation (2.1). Secondly, this model reduces to the model of Waleffe (1997) in the special case in which the amplitudes A1,…,A8 are taken to be periodic in the z-direction, after applying appropriate projection onto orthogonal Fourier modes in z. We discuss each of these important issues in more detail in the following subsections.

### (a) Nonlinear terms and energy conservation

In this subsection we show that the nonlinear terms in equations (2.6)–(2.13) do not contribute to the energy budget for the dynamics. This is a physically crucial property in constructing any reasonable reduced model for shear flows: the only energy source term must be the body force F(y) that drives the laminar flow, and the only source of dissipation must be (linear) viscous diffusion.

We define the (dimensionless) total kinetic energy of the flow to be 2.14where Ω=[0,π/α]×[−1,1]×[0,Lz] is (one half of) the domain (in x,y,z coordinates) occupied by the fluid. We substitute our ansatz (2.5) into equation (2.14) and carry out the x and y integrals, noting that Fourier orthogonality enables us to remove all cross-terms except those involving A7 and A8 since they have the same Fourier dependencies; their contributions to u are Hence the kinetic energy E is given by 2.15which defines the ‘local’ energy quantity . After integrating several terms by parts, and noting that the boundary contributions vanish (since we use a periodic boundary condition in z), and also after some manipulation of the cross-terms involving A7 and A8, we obtain the expression We can now compute the time evolution of the kinetic energy by differentiating with respect to time. After carrying out further integrations by parts we find 2.16We now substitute for the time derivatives using the PDEs (2.6)–(2.13). The interest in pursuing this calculation is that at this stage we find that the conservative nature of the quadratic nonlinearities in equations (2.6)–(2.13) becomes apparent since all the cubic terms cancel out (after appropriate integrations by parts). We are then left with a single linear source term and a collection of quadratic dissipation terms: The derivative of this equation for the evolution of the total kinetic energy, which describes the balance between the driving provided by the body force term F(y) and viscous dissipation, makes it clear that the nonlinear terms in equations (2.6)–(2.13) conserve energy. This justifies the self-consistency of the selection of just the eight modes A1,…,A8 in the reduced model: by including exactly these modes, no unphysical effects are introduced into the energy evolution.

### (b) Relation with the model of Waleffe (1997)

Waleffe (1997) considered a far simpler modal truncation in which each mode contains only a single Fourier mode in z. This introduces an additional parameter: the wavenumber γ in the z-direction, but the model comprises ODEs rather than PDEs for the eight amplitudes and is therefore much more amenable to analysis. The derivation of such a modal truncation is slightly simpler than that described above since one can use Fourier orthogonality directly in the z-direction as well as in x and y. In fact, the reduced model (2.6)–(2.13) reduces exactly (after rescalings of the amplitude variables) to that discussed by Waleffe if one inserts the corresponding trigonometric dependencies and then projects out unwanted modes that arise from some of the nonlinear terms. This projection step is self-consistent in the sense that in both equations (2.6)–(2.13) and the ODEs derived by Waleffe the nonlinear terms conserve energy. Table 1 lists the Fourier dependence in the z-direction that Waleffe (1997) assumed for each mode.

Waleffe observed that his ODE model had a number of deficiencies, for example, the streaks (described by A2) were unstable to x-dependent perturbations, with even and odd symmetry in y, only if the wavenumber γ was sufficiently large: γ2>α2 and γ2>α2+β2, respectively, to be precise.

More seriously, he discusses the inadequacy of his ODE model in capturing the interaction between the mean shear (M, or equivalently A1) and the x-dependent modes (A, C, B and D, or equivalently A4,…,A7) that arise from consideration of advection by the mean shear. We observe that in equations (2.6)–(2.13) these quadratic interactions take far more complex forms that in many cases vanish identically when the simple Fourier mode dependencies in z that are listed in table 1 are imposed. In particular, the mean shear A1 is influenced by new combinations of x-dependent modes which do not arise in the ODE model because the expressions A6A7′′−A6′′A7 and A4′′A5A4A5′′ which are present in equation (2.6) vanish identically in the ODE reduction. Terms with a structure identical to this (i.e. of the form AnAm′′−An′′Am) appear in several of the other amplitude equations. Compared with Waleffe's ODE model, the other qualitatively new couplings introduced in equations (2.6)–(2.13) are the term β(A3A6′′)′ in equation (2.11) and the term −2α2βA1A6 in equation (2.13).

## 3. Numerical results

In this section we present the results of time-stepping the system of PDEs (2.6)–(2.13) over the range 50≤Re≤200. Our numerical method is the pseudospectral exponential time-stepping scheme referred to as ‘ETD2’ by Cox & Matthews (2002), using 128 Fourier modes in z, and suitably small timesteps such that our results were insensitive to the timestep used.

As in Moehlis et al. (2004) our primary interest is in the emergence of a chaotic saddle in phase space. This is indicated by the lifetimes of chaotic transients as trajectories evolve towards the linearly stable laminar state, having started from initial conditions that are far from the laminar equilibrium.

### (a) Initial conditions and domain parameters

The lifetimes of chaotic transients are of course sensitive to the choice of initial condition, and in order to investigate the response of the flow to a spatially localized perturbation we took initial conditions corresponding to a Gaussian profile, scaled so that the four modes A3, A4, A5 and A6 gave equal contributions to the initial kinetic energy E0. Later, in §3(f), we compare these results with those obtained using a sinusoidal initial condition.

Specifically our Gaussian initial condition takes the form A1=A2=A7=A8=0 and 3.1for j=3,…,6, where the coefficient σ describes the width of the Gaussian, and the normalization constants cj are given by 3.2 3.3 and 3.4 The differences in the expressions for c3,…,c6 reflect the different contributions made by A3,…,A6 to the kinetic energy (2.15).

We present results for a domain size Lx=1.75π, Lz=1.2π corresponding to the ‘minimal flow unit’ identified by previous authors (Hamilton et al. 1995, and adopted by Moehlis et al. 2004) as the smallest domain in which sustained spatio-temporally chaotic dynamics have been found numerically. We consider values for σ in the range 0.2≤σ≤0.8, initially setting σ=0.2, so that the initial Gaussian disturbance is always spatially well localized in the z-direction.

### (b) Results at fixed Reynolds numbers

Our numerical integrations are carried out until either the dynamics approaches very close to the laminar state, or until 1000 dimensionless time units of h/(2U0) have elapsed: this is the maximum lifetime that transients are followed for in our computations. We follow transients for the range 50≤Re≤200, increasing Re in steps of unity, and increasing initial kinetic energy E0 in the range 0<E0<3.0 in steps of 0.01.

At low Reynolds numbers, Re<117, we find numerically that all initial conditions decay monotonically towards the stable laminar equilibrium, with a lifetime that increases slowly with E0 in the range 0<E0<1.7 approximately, and then decreases slowly as E0 increases further. For Re≥117 the dynamics changes abruptly and initial energies around E0=1.7 show far longer transients, as indicated in figure 2 which shows the lifetimes of trajectories for four fixed values of Re, as E0 varies. The overall impression of figure 2 is of a well-defined transition from laminar to spatio-temporally complex dynamics for a range of initial perturbations. It is clear that windows of rapid attraction to the laminar state remain, for example, for Re=143 and initial energies around E0=1.3.

Figure 2.

Lifetimes of transients started from the Gaussian initial condition, for initial kinetic energies 0<E0<3.0 for Re=108 (black line with plus symbols), Re=117 (red line with filled circles), Re=128 (blue line with filled diamonds) and Re=143 (magenta line with filled squares). Parameters are σ=0.2, Lx=1.75π and Lz=1.2π. (Online version in colour.)

### (c) Transition boundary and structure of the chaotic saddle

Figure 3.

Lifetimes of transients started from the Gaussian initial condition, over the range 0<E0<3.0 and 50≤Re≤200. Parameters are σ=0.2, Lx=1.75π and Lz=1.2π.

Figure 4.

Enlargements of figure 3 showing a lack of coherent organization to the lifetimes at higher resolution in Re and E0. (a) Transient lifetimes over the range 1.6<E0<1.7 and 130≤Re≤140. (b) Transient lifetimes over the range 1.65<E0<1.66 and 135≤Re≤136 as indicated by the black square in (a). For both figures, the other parameter values are σ=0.2, Lx=1.75π and Lz=1.2π. (Online version in colour.)

The lower boundary of the chaotic saddle is of particular physical interest, since it describes the smallest amplitude perturbation required to produce an extended chaotic transient. For the particular form of perturbation used here, in the case σ=0.2, we find numerically that the lower boundary in figure 3 is remarkably well fitted by a power-law curve as illustrated in the log–log plot in figure 5. Additional computations show that a power law continues to provide a very good fit as Re is increased, up to at least Re=103, as shown in figure 6a. The black curve in figure 6a is the best-fit curve 3.5Blue squares in figure 6 indicate that the flow returned to laminar before the computational time limit Tmax was reached, while red plus signs indicate that the flow did not relaminarize before t=Tmax. We set Tmax=5000 for 200<Re≤750 and Tmax=10 000 for 800≤Re≤950.

Figure 5.

Log–log plot of transient lifetimes over the range 0.1<E<3.0 and 100≤ Re≤200 to show power-law form of the lower boundary of the attractor. The solid red line indicates the power-law E0=(125/Re)4.65. Parameter values are σ=0.2, Lx=1.75π and Lz=1.2π.

Figure 6.

Estimates of the lower boundary of the chaotic saddle over the range 200<Re<1000. (a) σ=0.2; squares (blue) indicate a rapid return to the laminar state. Plus symbols denote a long-lived chaotic transient. The black line indicates the power-law scaling of the lower boundary. (b) Best-fit power laws for a range of values of σ, showing a broad insensitivity to the width parameter. Parameter values are Lx=1.75π and Lz=1.2π (solid black line, σ=0.2; solid blue line, σ=0.3; solid red line, σ=0.5 and dotted line, σ=0.8). (Online version in colour.)

In figure 6b, we show the best-fit power law for the boundary between laminar and spatio-temporally complicated flows for different widths of the Gaussian initial condition (3.1). All are excellent fits to the data over the range 200<Re<1000, and the results are identical for computations using 128 modes or 32 modes in z. Since the slopes of the power laws become more negative as the coefficient increases, the overall envelope of perturbation energies as Re increases, formed by taking the minimum value over all four curves, has contributions from each curve at sufficiently large Re. This envelope is an estimate of the true lower boundary for which we would ideally compute the minimum over all possible forms of initial perturbation.

For perturbations with a Gaussian profile (3.1), and for the range of σ considered here, it appears unlikely that the differences in exponents would be detectable experimentally over the range of relevant and accessible Re. Table 2 gives the details of the best-fit parameters of the power-laws shown in figure 6b, using a functional form E0=(c/Re)p.

View this table:
Table 2.

Dependence of power-law scalings on the width σ of initial condition.

Moehlis et al. also computed the distribution of lifetimes of turbulent transients at fixed initial energies and Reynolds numbers. They found that the survival probability P(T) that the solution had not decayed back to the laminar state after a time T was distributed exponentially, with a mean lifetime that increased rapidly with Re above the critical value at which the chaotic saddle appeared. Numerically, we construct a lifetime distribution by distributing the initial energy randomly across a subset of the modes. For our spatially extended model there are two possible approaches to randomizing the distribution of energy.

In the first approach, the total initial energy is distributed uniformly on a spherical energy shell. This can be achieved straightforwardly by modifying the expressions (3.2)–(3.4) for the normalization coefficients by replacing E0 with in the expression for coefficient cj, where j=3,…,6. The coordinates ξj are those of points distributed at random over the unit sphere in given by . We refer to this as randomizing over the amplitudes of the perturbations.

In the second approach, the central position of each Gaussian perturbation Aj is shifted randomly in z away from the centre of the domain Lz/2, i.e. the values of the coefficients c3,…,c6 are left unchanged but the centre of the Gaussian in equation (3.1) is modified for each amplitude Aj. Since we employ periodic boundary conditions, the total energy in the perturbation is preserved. We refer to this procedure as randomizing over the locations of the perturbations. Figure 7 shows the results of the first approach at two Reynolds numbers quite close to the transition boundary. At low lifetimes (T<200) the survival probability remains unity since any transient takes at least this finite amount of time to decay to be within the prescribed threshold of the laminar state. As T increases there is an initial sharp drop indicating that a substantial proportion of trajectories do decay rapidly, with lifetimes in the range 200<T<400. At larger T the distribution becomes close to a straight line on the linear-log plot, which is consistent with an exponential distribution for T≥400 (approx.), as indicated by the dashed line. For Re=130 the best-fit coefficient values are a1=3.93×10−1 and a2=8.65×10−4. For Re=140 the best-fit coefficient values are a1=4.57×10−1 and a2=5.11×10−4. Similarly, figure 8 shows the lifetime distribution computed from the second approach, randomizing over the locations of the perturbations. The distribution shows very similar characteristics to that in figure 7a and is again well described by an exponential distribution with best-fit parameters a1=4.24×10−1 and a2=8.32×10−4. We observe that the values for the coefficient a2 are very similar between the two randomization methods. In all cases, we expect that as Re increases, the distribution exhibits a systematic deviation from exponential and flattens out. We anticipate that, as in the ODE case, this is due to the appearance of stable attractors near the chaotic saddle, as Re increases. These attracting sets then absorb trajectories with positive probability, leading to a proportion of trajectories which never return to the laminar state (Moehlis et al. 2005).

Figure 7.

Distribution of lifetimes P(T) computed over an ensemble of 2000 initial conditions produced by randomizing over amplitudes of perturbations at E0=1.0. (a) Re=130. (b) Re=140. Parameter values are σ=0.2, Lx=1.75π and Lz=1.2π. (Online version in colour.)

Figure 8.

Distribution of lifetimes P(T) computed over an ensemble of 2000 initial conditions produced by randomizing over the locations of perturbations at E0=1.0 and Re=130. Parameter values are σ=0.2, Lx=1.75π and Lz=1.2π. (Online version in colour.)

### (e) Spanwise resolution

The very high resolution used in the z-direction (N=128 Fourier modes) is greater than required, for such a small domain, in order to capture the dynamics of the ‘active’ modes of the system. In order to probe the range of wavenumbers that make substantial contributions to the dynamics we investigated the effect of varying the truncation level of the numerical scheme in z. We define the ‘m-mode truncation’ of the PDEs (2.6)–(2.13) by keeping the Fourier modes ∼e±inγz for 0≤nm−2 for A1, A4 and A5, and keeping the modes for 0≤nm−1 for A2, A3, A6, A7 and A8. The m-mode truncation is thus a collection of (at most) 16m−14 real ODEs. This is an upper bound on the effective dimension of the ODE dynamics since not every real and imaginary part is coupled for every m.

By construction, this further reduction resembles the Waleffe model in the particular case m=2, but by varying m we are able to probe the influence of the higher-wavenumber modes on the location of the lower boundary for the onset of temporally complex dynamics. For each value of E0 and m, a single initial condition was used. This initial condition was the projection of the Gaussian profiles described in §3(a) onto the available Fourier modes in z. Results are shown in figure 9. Figure 9a, for the small domain Lz=1.2π, shows that computations with truncation levels m≥5 produce an identical indication of the boundary between the laminar and spatio-temporally complex behaviour. The case m=4 is qualitatively, but not quantitatively, correct. For the larger domain Lz=6π, figure 9b indicates that the numerical results are essentially unchanged for m≥17. The number of Fourier modes required to maintain accuracy in the larger domain is therefore broadly in line with, although slightly lower than, what one might naively expect.

Figure 9.

Lifetimes of transients T as a function of initial energy E0 for Re=200 in domains of widths (a) Lz=1.2π (filled black square and solid line, m=16; filled red circle and solid line, m=9; filled blue diamond and solid line, m=5; filled magenta small square and solid line, m=4; filled black square in a dashed line, m=3 and filled red circle in a dashed line, m=2). (b) Lz=6π, showing the effect of varying the truncation level m (filled black square and solid line, m=32; filled red circle and solid line, m=19; filled blue diamond and solid line, m=17; filled black square in a dashed line, m=15; filled red circle in a dashed line, m=13; filled blue diamond in a dashed line, m=11; filled magenta small square in a dashed line, m=9 and filled red circle in a dash-dotted line, m=5). Vertical axis is logarithmic for clarity, and shows the additional lifetime after subtracting a constant: T−365 in (a), T−414 in (b). In (a) the data for m=9 and m=5 lie exactly on the m=16 values: they are artificially offset vertically for clarity. In (b) the data for m=19, m=17 and m=15 lie exactly on the m=32 values and are similarly offset for clarity. Computations were terminated at Tmax=5000. Parameter values are σ=0.8 and Lx=1.75π. (Online version in colour.)

Finally we note that as m decreases further, the boundary appears consistently to move to lower E0. This behaviour is purely a function of the organization of invariant sets in the phase space; it is not clear that there is a physical reason why these should appear to move in one direction or the other.

### (f) Spatio-temporal dynamics and effect of initial conditions

In this final subsection, we comment on the spatio-temporal dynamics of the PDEs for initial energies near the transition boundary, and we compare the results of §3(c) for the lifetimes of transients with those in this section obtained using sinusoidal initial conditions. Figure 10 shows the spatial and temporal evolution of the PDEs for Re=200 and E0=0.11, just above the transition boundary. The quantity is plotted in the figure, so that the laminar state is at level zero, where is defined in equation (2.15). The initial monotonic decay towards the laminar state is interrupted by the growth of oscillatory disturbances. These disturbances generate a sharp spike in the kinetic energy, localized in both space and time, before the solution settles into a spatio-temporally chaotic state. For comparison, at E0=0.1 we observe only a monotonic decay to the laminar state.

Figure 10.

Space–time plot of the local energy quantity indicated by both surface height and colour, for the PDEs (2.6)–(2.13) for Re=200 in the minimal flow unit domain. is defined in equation (2.15). The additional square terms shift the laminar state to the zero level in the plot and serve to highlight the localized ‘edge’ state at t≈450. Parameter values are E0=0.11, σ=0.2, Lx=1.75π and Lz=1.2π. (Online version in colour.)

Analysis of the energy distribution across Fourier modes shows that the small-scale oscillations just prior to the localized spike involve many higher-wavenumber Fourier modes. These mode amplitudes grow very rapidly as we approach the spike state. Since the formation of the spike is, in almost every case, the ‘edge state’ that is the precursor to spatio-temporally complicated dynamics, one interpretation of the results in the previous subsection is that the spatial resolution required to correctly determine the boundary of the basin of attraction of the laminar state is indicated by the resolution required to describe accurately these spatially localized spikes.

More generally, the dynamics near the boundary between monotonic relaminarization and spatio-temporal complexity appear to depend on the evolution of high-wavenumber modes, seeded by the use of a Gaussian initial condition which injects energy into every available mode. This is in contrast with the sinusoidal initial conditions used in previous reduced models, where, by construction, such a sinusoidal initial condition was the only possible choice.

In figure 11, we show the analogous plot to figure 3 for the lifetimes of transients, but in this case using a sinusoidal initial condition instead of a Gaussian profile. For the sinusoidal initial condition, we set ; ; A4=c4; A5=c5; A1=A6=A7=A8=0 where the constants c2,…,c5 are given by so that the initial kinetic energy E0, defined in equation (2.15), is distributed equally between the four non-zero amplitudes.

Figure 11.

Lifetimes of transients started from the sinusoidal initial condition, over the range 0<E0<3.0 and 100≤Re≤200. Solid (red) line indicates the boundary E0=(75.7/Re)4.3 corresponding to the lowest of the boundaries in figure 6b for Gaussian perturbations. Parameter values are Lx=1.75π and Lz=1.2π. Computations used N=32 Fourier modes in z and were terminated at Tmax=1000.

Comparing figures 3 and 11 there are important qualitative and quantitative differences. Firstly, the boundary between relaminarization and spatio-temporal complexity is much more obvious in figure 3 where we employ the Gaussian initial condition. The boundary in figure 11 shows much more ‘structure’; it is much less obvious that a boundary in the sense, for example, of a countable collection of (piecewise) continuous curves in the (Re,E0) plane, can even be defined. Many deep valleys of short lifetimes persist up to Re=200 and beyond. Secondly, the solid (red) line in figure 11 shows the lowest boundary computed for Gaussian perturbations, corresponding to σ=0.5, see table 2. It appears that the laminar state is substantially more sensitive to Gaussian perturbations than those of the same energy but in the form of low-wavenumber sinusoids. Tentatively, based only on the data in figure 11, we suggest that the lower boundary of the appearance of spatio-temporally complex dynamics arising from the sinusoidal perturbation scales as E0Re−2, i.e. perturbation amplitude scaling as Re−1. This is a typical exponent produced by very many reduced models of a very low order, as summarized and discussed by Baggett & Trefethen (1997).

Within the region of spatio-temporally complex dynamics, the lifetimes in figure 11 and in additional enlargements (not shown here) show a degree of correlation between lifetimes at neighbouring points in the (Re,E0) plane which is not nearly so clearly shown to exist in figure 3 or in the enlargements shown in figure 4.

In summary, in these respects figure 11 is reminiscent of lifetime plots for the very low-dimensional models of Moehlis et al. (2004), see their figures 6 and 8, and figures 5 and 6 in the paper by Eckhardt & Mersmann (1999). We conclude that allowing for the accurate representation of spatially localized initial conditions by extending the spanwise resolution of the model generates results that differ substantially, both qualitatively and quantitatively, from those of previous reduced models.

## 4. Discussion and conclusions

In this paper we have presented an extended version of the Galerkin-truncated model due to Waleffe (1997) for the transition to turbulence (or, at least, spatio-temporally complex dynamics) in sinusoidal shear flow. This model is appealing since it provides an intermediate step between previous analytical work and direct numerical simulation (DNS) of the full Navier–Stokes equations. Preserving full resolution in the spanwise (z) direction and removing the assumption of periodicity allows both the use of spatially localized initial conditions, and the (transient) formation of localized structures in the flow which (although unstable) are known to exist and play a role in mediating the onset of turbulence. The use of a small number of Fourier modes in the wall-normal (y) and streamwise (x) directions provides the simplification of the underlying Navier–Stokes equations, which in turn allows us to perform very detailed investigations of the dynamics of this reduced model.

We compare our results with those of previous authors, in order to see which properties are common to these different approaches, and which are not. For example, we find, in agreement with the results of Moehlis et al. (2004), that the lifetimes of turbulent transients are well described by an exponential distribution. However, our results show that the transition boundary, while exhibiting some of the structured shape observed by many authors (including, in the case of pipe flow, Schneider et al. 2007), appears at lower perturbation energies, and much more abruptly, than for the ODE models investigated by Eckhardt & Mersmann (1999) and Moehlis et al. (2004). The PDE model that we present here is able to represent both spatially localized and spatially extended initial conditions and therefore we are able to make direct comparisons of this kind.

Our key finding is that spatially localized initial conditions are able to provoke complicated behaviour at substantially lower energies than the sinusoidal, spatially extended perturbations used in previous studies. Moreover, the perturbation energy at the lower boundary of the chaotic saddle appears to scale as Rep with the exponent p≈−4.3, rather than Re−2 as in Eckhardt & Mersmann's 19-mode truncated ODE model (note that their figure 5 showing an Re−1 power-law plots Re against mode amplitude which is proportional to ). In the present work, the exponent in this power-law scaling was found to depend only weakly on the width of the Gaussian perturbation used.

In addition, our results are robust to the numerical resolution used in the spanwise direction, and, for a relatively small domain of width Lz=1.2π, point to the necessity of keeping around five Fourier modes in z in order accurately to capture the dynamics of the fully resolved PDE model. One possible explanation of these results is that admitting higher-wavenumber modes generates many more invariant sets within the boundary of the basin of attraction of the laminar state. Then, even small amounts of initial energy in these modes forces the system to spend much longer in the vicinity of these sets before being able to relaminarize. In this sense, ‘holes’ in the basin boundary are filled in. The existence of these new invariant sets, and the lack of ‘holes’, leads to a robustness in the lengths of transients, and therefore to a more clearly defined boundary between monotonic relaminarization and long-lived transients.

It would clearly be of interest in future work to look at the relationship between localized states which have been observed and studied in some detail in DNS for shear flow problems (Schneider et al. 2010a,b) and the dynamics of the reduced model presented here. We anticipate that the reduced model contains such states, and the homoclinic snaking bifurcation diagrams that typically organize them in driven dissipative systems such as shear flows, just as model ODE truncations, for example that discussed in Moehlis et al. (2005), contain equilibria and time-periodic solutions very similar to those located in DNS (Nagata 1990; Gibson et al. 2009). It should be possible systematically to further reduce the model equations presented here in order to make direct connections between theoretical work on localized states (Burke & Knobloch 2006; Chapman & Kozyreff 2009; Dawes 2010) and the DNS results referred to above. In turn, the identification and analysis of additional unstable invariant sets within the boundary of the basin of attraction of the laminar state (as discussed by Lebovitz 2009), and their parameter dependence, will greatly help our understanding of the process of relaminarization.

## Acknowledgements

J.H.P.D. would like to thank Rich Kerswell and Tobias Schneider for useful conversations, and Matthew Chantry for a minor correction. Both authors are grateful to the anonymous referees for very useful comments, and they gratefully acknowledge financial support from the Royal Society; J.H.P.D. currently holds a Royal Society University Research Fellowship.