A low-dimensional model for turbulent shear flows

We analyse a low-dimensional model for turbulent shear flows. The model is based on Fourier modes and describes sinusoidal shear flow, in which the fluid between two free-slip walls experiences a sinusoidal body force. The model contains a total of nine modes, including modes describing the basic mean velocity profile and its modification, downstream vortices, streaks, and instabilities of streaks, with other modes being a consequence of the non-linear interactions. The transition to turbulence for the model is subcritical and intermittent, and the distributions of turbulent lifetimes are exponential, in agreement with observations in many shear flows.


Introduction
The understanding of the transition to turbulence in many flows has benefited from the analysis of simplified models. Beginning with Lorenz's numerical simulations of a three-mode caricature of thermal convection [1], various models for convection and centrifugally unstable systems have been derived and analysed. Such models highlight different kinds of bifurcations and a wide variety of dynamical behaviour in systems with a linear instability [2]- [5]. However, describing the transition to turbulence in shear flows such as plane Couette flow or pipe flow has proven much more difficult: there is no linear instability around which centre manifold and other reduction techniques can be applied; moreover, numerical simulations and experiments indicate that the dynamics are spatially and temporally complicated right above the onset of turbulence [6]- [10].
Although there is considerable complexity in the dynamics of shear flows, there is also some simplicity. Indeed, visual inspection of direct numerical simulations [11] and experiments in the transition region [10,12], proper orthogonal decomposition of the flow field [5,13], and results from linear stability analysis [7,14,15] all suggest that downstream vortices should play a major role. These modes drive streaks, and, as discussed perhaps most strongly by Waleffe [14,16], the unstable modes of these streaks are recycled as downstream vortices. This suggests that it should be possible to close the dynamics of shear flow turbulence with a few modes. What makes the development of such low-dimensional models complicated is the fact that the peculiar form of the non-linearity of the Navier-Stokes equations constrains the possible interactions among modes, so that a poor choice of modes in a model can lead to trivial dynamics.
There have been several previous attempts to derive low-dimensional models for shear flow turbulence. Because of their geometric simplicity, much of this effort has been in the modelling of plane Couette flow, in which fluid is sheared between two infinite parallel no-slip walls moving at the same speed but in opposite directions, and for sinusoidal shear flow, in which fluid between two free-slip walls experiences a sinusoidal body force which gives rise to a sinusoidal laminar profile. For example, building on his analysis of linear instabilities, Waleffe proposed an eight-mode model and a further reduction to a four-mode model for sinusoidal shear flow [16]. The dynamics of the eight-mode model were not studied in detail. For the four-mode model, the 'turbulent state' is represented as a stationary point or a periodic solution [14,16,17]. On the other hand, analysis of a straightforward projection of the sinusoidal shear flow equations onto a set of 19 Fourier modes gives a fluctuating turbulent state [18]. Eliminating a residual symmetry reduces the number of degrees of freedom to nine without a change in qualitative behaviour [19]. Different kinds of dynamical behaviour, including heteroclinic cycles, have been observed in a model for plane Couette flow derived using modes extracted from a proper orthogonal decomposition [13]. We also note that for pipe flow, Brosa and Grossmann [20] studied a model with 42 Galerkin modes selected from the couplings between the linear eigenmodes of the Stokes problem. Finally, phenomena not unlike those observed in shear flows have also been seen in models with couplings that maintain certain characteristics of the flow, but are not derived from the Navier-Stokes equations [21]- [23].
The model that we present below has nine degrees of freedom. All but one of the modes are already contained in the eight-mode model proposed by Waleffe [16], although there are small but significant modifications to some modes. The essential elements are modes for the basic mean velocity profile, the downstream vortices, the streaks, the instabilities of the streaks, and several modes that result from the advection of the previous modes by the mean profile.
The key improvement on Waleffe's eight-mode model [16] is that we include a mode for the modification of the structure of the mean velocity profile: as is well known for turbulent shear flows, the fluctuations change the mean profile, and this is described to lowest order by the new mode. In addition, unlike the modes in [16], the variation in the shear direction in some modes is not a simple sine or cosine function, thus introducing additional couplings between the modes. As a result, we are able to capture desired aspects of shear flow turbulence with a much lower dimensional model than the 17-dimensional extension suggested in [16].
In the present paper we focus on the dynamics of the model. The dynamical system behaviour of the model, and the analysis of the stationary states and periodic solutions and their bifurcations, will be presented in a companion paper [24]. The outline is as follows: after general remarks in section 2 we introduce the modes, their couplings, and the equations of motion in section 3. We analyse the dynamics in the minimal flow unit, then for a wider and longer domain in sections 4.1 and 4.2, respectively. The characteristics of the transition to turbulence are studied in section 5. We conclude with a few remarks in section 6.

Sinusoidal shear flows
We use a coordinate system with x pointing downstream, y in the direction of the shear, and z in the spanwise direction. The characteristic velocity U 0 is taken to be the laminar velocity arising due to the forcing at a distance d/4 from the top wall, where d is the distance between the walls; see (5) and (6) with Reynolds number defined to be where ν is the kinematic viscosity. The fluid is incompressible, and there are free-slip boundary conditions at the walls at y = ±1, i.e., Finally, the flow is assumed to be periodic in the streamwise and spanwise directions, with lengths L x and L z , respectively. Following Waleffe [16], we take the non-dimensionalized volume force to be giving the laminar profile The laminar profile is inflectional but remains stable for all Reynolds numbers Re [25]. In the following, we let α = 2π/L x , β = π/2, and γ = 2π/L z , and denote the domain 0 x L x , −1 y 1, 0 z L z by .

A low-dimensional model
Our model is a nine-mode generalization of the eight-mode model considered in Waleffe [16]. The additional ninth mode describes a modification to the basic mean velocity profile by the turbulence. It has wave vector k = (0, 3π/2, 0), to be compared with the wave vector k = (0, π/2, 0) of the basic profile (6). Since, in general, if two modes with wave vectors k 1 and k 2 enter into (u · ∇)u, a mode with wave vector k 1 + k 2 results, generation of such a modification of the basic profile requires modes with y-component of their wave vector equal to π. Allowing such modes is the second major difference from the model of Waleffe.
The modes and wave vectors that should be included appear rather naturally, as follows. A prominent feature of shear flows is the presence of streamwise vortices. Vortices that span the entire gap have a wave vector k = (0, π/2, γ), as for u 3 (the complete modes are given below, equations (7)-(16)). A streak, u 2 , capturing spanwise variation of the streamwise velocity, results from the advection of the basic profile u 1 by the vortex, and hence has wave vectors (0, 0, γ) and (0, π, γ). In fact, these Fourier modes are tied together through a cos 2 (πy/2) = (1 + cos(πy))/2 dependence. Then, advecting the streak by the vortex gives a mode including the wave vector (0, 3π/2, 0), as in the new mode u 9 . In [16], none of the modes, including the streak, have wave vectors with y-component equal to π, and hence that model cannot produce a modification of the shape of the basic profile.
The other modes in [16] are associated with instabilities of the streak mode, with their specific forms determined by physical and symmetry considerations. For example, the streak can develop instabilities with a component representing vortices pointing in the wall-normal direction, described by wave vectors (±α, 0, ±γ). Because our streak mode has an additional cos 2 (πy/2) dependence with respect to Waleffe's streak mode, we modify the appropriate streak instability modes u 4 and u 6 to have the same dependence. The resulting modes then couple to other modes to generate the modification to the basic profile with wave vector (0, 3π/2, 0).
Finally, as in [16], we 'pin'the spanwise and streamwise locations of the modes, which gives only one real degree of freedom for each mode. Because of the translation symmetries in the downstream and spanwise directions, there is actually no preferred location for these structures; proper accounting for these symmetries would require two or more real (or one or more complex) degrees of freedom for each mode in our model, except for u 1 and u 9 , thereby greatly increasing the dimensionality of our model, cf [13,26]. Consideration of such an extension of the present model is deferred to future work.
With the above considerations in mind, the modes for the model are: the basic profile the streak the downstream vortex and the spanwise flows and Furthermore, we have the normal vortex modes and generated from the advection of u 4 and u 5 by both the streak u 2 and the vortex u 3 . These interactions also give rise to a fully three-dimensional mode, πα sin(αx) sin(πy/2) sin(γz) 2(α 2 + γ 2 ) cos(αx) cos(πy/2) sin(γz) −πγ cos(αx) sin(πy/2) cos(γz) Matrix of couplings between the modes: (u i · ∇)u j with indices i from the row and j from the column contributes to mode k, with k the entry in the matrix. Dashes indicate couplings that vanish exactly. Boldface entries indicate couplings that are absent in Waleffe's model. All pairs except for the entry (3,1) and, of course, the dashed ones, couple to higher-order Fourier modes outside the set contained in our model. The lack of coupling within this set is indicated by * .
Finally, we have the new mode mentioned earlier, the modification of the basic profile The modes are orthogonal, and normalized so that Each mode individually satisfies incompressibility and free-slip boundary conditions at the walls. Note that the modes with a cos 2 πy dependence are not pure Fourier modes, but rather superpositions of two Fourier modes. Introducing independent amplitudes for the pure modes would increase the number of equations, but also raise the risk of causing a collapse in the dynamics if the non-linear couplings between the modes are not strong enough. This issue has not been explored further. Table 1 shows how these modes couple through the advection term to generate other modes. In the ith row and jth column, a dash indicates that (u i · ∇)u j vanishes exactly, an entry k indicates a contribution to mode k since and * indicates that (u i · ∇)u j has only a non-trivial projection to modes outside the set contained in our model. Bold entries are couplings that are not contained in Waleffe's eightmode model [16]. The asymmetry in the couplings indicates that sometimes one mode advects another but not vice versa, i.e., (u i · ∇)u j can be different from (u j · ∇)u i . Ordinary differential equations are obtained by Galerkin projection [5,16], i.e. by plugging the ansatz with real amplitudes a m into (1), multiplying by u n (x), integrating over , and using (17). This gives The pressure term makes no contribution here because the modes are divergence-free, their y-component vanishes at the walls, and they satisfy periodic boundary conditions in x and z, cf [5]. For our choice of modes, we obtain the amplitude equations where These equations have a strong similarity to Waleffe's eight-mode model [16], but because some modes differ there are different O(1) factors multiplying some terms, and several additional terms, including all terms that contain a 9 . While the numerical factors can be expected to play a minor role for the qualitative behaviour, all couplings seem to be significant. For example, for Waleffe's eight-mode model we have been unable to locate any unstable, stationary finite amplitude coherent structures consisting of streamwise vortices and streaks, as have been found numerically for shear flows such as plane Couette flow [27,28]. On the other hand, for our nine-mode model, with its additional couplings, we do find such solutions: see figure 1.
The laminar state corresponds to the fixed point with a 1 = 1, and a 2 = · · · = a 9 = 0, which is linearly stable for all Re. In the next section, we will see that turbulence exists for the model despite the stability of the laminar state: the model thus captures the subcritical nature of the transition to turbulence.

Dynamics of the model
In [11], plane Couette flow is considered for a minimal flow unit, the smallest domain which is numerically found to sustain turbulence. It elucidates a 'self-sustaining process' in which streamwise vortices cause streak formation, then the streaks break down to give x-dependent flow, then the streamwise vortices regenerate and the process repeats. It has been suggested that such a self-sustaining process is a universal characteristic of shear flow turbulence [11,16,29,30]. In section 4.1, we analyse the model in the highly constrained geometry of a minimal flow unit.
On the other hand, the analysis of Nagata [27] and Busse and Clever [28] for the optimal domain size for the formation of stationary coherent structures suggests considering a longer and wider domain. The dynamics for that domain are given in section 4.2.

Domain with
We first consider a flow domain with L x = 1.75π and L z = 1.2π. For plane Couette flow, this corresponds to the minimal flow unit, the smallest domain which is found numerically to sustain turbulence [11]. Previous models for shear flows for this domain size include [26] for plane Couette flow. See also [31] for a discussion of models for minimal flow unit channel flow.
The dynamics for the model at Re = 400 are shown in figure 2. The movie clearly shows the following series of events: if we start with a vortex (indicated in the upper right panel by the arrows for the velocity field in the plane), the formation of a streak can be observed next (indicated by the colours for the downstream component). The streak then goes unstable by the sinusoidal mode that is the basic instability kept for the model [16]: it shows up prominently in the midplane cross section of the velocity field. The sinusoidal modulation grows, breaks up into patches and rearranges with the regions of vertical motions shifted by half a wave length: the consequences for the vortices are that they disappear and then reform with a reversed sense of rotation. the eventual decay to the laminar state shows it does not belong to a turbulent attractor but rather to a transient state on a chaotic saddle. Notice that the amplitude of a 1 becomes large near t ≈ 1600, indicating an attempt to leave the turbulent state and to relax to the laminar one. However, the conditions for relaxation to the laminar state are not only that a 1 has to become large, but also that the other modes have to be small, which is not satisfied until t ≈ 2600 when the trajectory finally escapes the turbulent region.

Domain with
The moderately longer and wider domain with L x = 4π and L z = 2π has been considered in many previous studies of shear flow turbulence, including plane Couette flow [13,19,27,28] and sinusoidal shear flow [19]. This corresponds to an optimal domain size for plane Couette flow in the sense that the steady finite amplitude solutions appear at the smallest value of Re [28,32]. The dynamics in the wider domain is similar to that for the smaller one, but perhaps a bit more violent. The reversals of the vortices are still present, but the dynamics of the sinusoidal instability in the midplane representation is less regular. In a full numerical simulation, more modes would be present and the variations in space and time can be expected to be even more irregular.
The timeseries for all the nine amplitudes in the model for the solution shown in figure 4 are displayed in figure 5. As for the minimal flow unit domain, the fluctuations are compatible with a non-periodic, irregular, transient turbulent state. Notice that the amplitude of a 1 becomes large several times, indicating an attempt to leave the turbulent state and to relax to the laminar one, e.g. near the times t = 800, 2700, 4500, 7800 and, most clearly, near t = 6500, where it almost succeeds. Finally, near t = 10 000 the flow laminarizes.

Transition behaviour
The finite lifetime of the trajectories for our model shown in the previous section indicates that the turbulent state is not an attractor. Similar behaviour has been seen in the other shear flows without a linear instability, and it has been suggested that the turbulent state is a chaotic saddle [9,18,33]. The typical characteristic of a chaotic saddle is an exponential lifetime distribution [34]- [36] that has also been seen in experiments and simulations [9,37].
For low values of the Reynolds number Re, the model shows an apparent fractal dependence of lifetimes on Re and on initial energy E of a perturbation ( figure 6). Here initial conditions are chosen so that a 2 = a 3 = a 4 = a 5 = √ E − 1/2, a 1 = a 6 = a 7 = a 8 = a 9 = 0. The strong variations of the lifetime with initial conditions suggests a probabilistic description, e.g. the study of the survival probability P(t) to be still turbulent after a time t. For the determination of P(t) for a set of Reynolds numbers we used ensembles of 1000 randomly selected initial conditions on an energy shell i a 2 i = 0.3 that were integrated for a maximum time of 10 000 units. The results in figure 7 support an exponential lifetime distribution, with a characteristic mean lifetime that increases rapidly with Reynolds number. However, for Re > 320 we have clear deviations from exponential decay. They appear first for long lifetimes, but set in earlier for increasing Reynolds number. The simplest explanation for this change is that a stable attractor appears and absorbs a finite fraction of initial conditions. Analysis of the periodic orbits in the system indeed shows this to be the case [24]. The lifetime distributions for the longer and wider domain with L x = 4π and L z = 2π are shown in figure 8; initial conditions are specified as for figure 6. The cuts at fixed Reynolds number in figure 9 and the successive magnifications in figure 8 indicate how strongly the lifetimes vary under tiny variations in parameters and how sensitively the system responds to changes in initial conditions. Survival probabilities for the larger domain are shown in figure 10. An exponential form is a good representation of the data for all Reynolds numbers. There does not seem to be a co-existing non-trivial stable state, as there was for the minimal flow unit, cf [24].
For both geometries, the median lifetimes increase rapidly with Reynolds number, perhaps exponentially, as indicated in figure 11. A rapid increase of lifetimes has also been observed in experiments [9,10], numerical simulations [37,38] and in other models for spatially extended systems [36,39].

Conclusions
The calculations and simulations show that nine judiciously chosen modes can capture the dynamics of vortices and streaks in turbulent shear flows. Interactions between the modes are suffice to keep the dynamics alive and non-periodic, as expected for a turbulent flow. For a wide range of parameters, the turbulent states do not persist indefinitely, but decay after a while. The survival probability is exponential, as expected for a chaotic saddle. The appearance of an attracting state in the minimal flow unit for large Reynolds numbers shows that a transition to an attractor is possible, but this depends on the domain size and most likely also on the number of modes included.
The model captures the dynamics of vortices and shows characteristics of the self-sustaining process elucidated in [11,16]. As suggested there, the cycle closes and the instability of the streaks drives modes that can be recycled to form vortices, thus completing the self-sustaining regeneration cycle.   Even the Reynolds number at which long-lived turbulence is found to be present is not too different from that found in full numerical simulations. For the free-slip boundary conditions, a critical Reynolds number of about 130 has been found in [40]. As can be expected for the more flexible free-slip boundary conditions, it is about 2/5 of the values for the transition in plane Couette flow [19,37]. That it is so close to the value we find here must be considered fortuitous, since experiments with other low-dimensional models shows wide variations in behaviour under changes in the number and characteristics of modes kept.