1 Introduction
The transition to turbulence in wall-bounded shear flows has been studied for well over a century, and yet only recently have experiments, numerical simulations and theory advanced to the point of providing a comprehensive understanding of the route to turbulence in such flows. Of late, research has focused on how turbulence first appears and becomes sustained. The issue is that typically wall-bounded shear flows undergo subcritical transition, meaning that as the Reynolds number is increased, turbulence does not arise through a linear instability of laminar flow, but instead appears directly as a highly nonlinear state. Moreover, the flow does not simply become everywhere turbulent beyond a certain Reynolds number. Rather, turbulence initially appears as localised patches interspersed within laminar flow. The resulting flow takes on a complex spatiotemporal form with competing turbulent and laminar domains. This, in turn, greatly complicates the quantitative analysis of turbulent transition in subcritical shear flows. See Barkley (Reference Barkley2016) and Manneville (Reference Manneville2016) for recent reviews.
In the 1980s the connection was developed between spatially extended dynamical systems and subcritical turbulent flows. This provided a broad and useful context in which to view turbulent–laminar intermittency. Kaneko (Reference Kaneko1985) constructed minimal models that demonstrated how dynamical systems with chaotic (‘turbulent’) and steady (‘laminar’) phases would naturally generate complex spatiotemporal patterns. Simple models were further studied by Chaté & Manneville (Reference Chaté and Manneville1988) amongst others. At the same time, Pomeau (Reference Pomeau1986) observed that subcritical fluid flows have the characteristics of non-equilibrium systems, exhibiting what is known as an absorbing state transition. Based on this, he postulated that these flows might fall into the universality class of directed percolation. This would imply that the turbulence fraction varies continuously with Reynolds number, going from zero to non-zero at a critical Reynolds number, with certain very specific power laws holding at the onset of turbulence. (These concepts will be explained further in § 2.) Since then considerable effort has been devoted to investigating these issues. The first experimental observation of directed percolation was reported by Takeuchi et al. (Reference Takeuchi, Kuroda, Chaté and Sano2007, Reference Takeuchi, Kuroda, Chaté and Sano2009) for electroconvection in nematic liquid crystals.
The status of our understanding for prototypical subcritical shear flows is as follows. For pipe flow, there are extensive measurements of the localised turbulent patches (puffs) that drive the transition to turbulence and we have a good estimate of the critical point for the onset of sustained turbulence (Avila et al. Reference Avila, Moxey, de Lozar, Avila, Barkley and Hof2011). However, currently there is no experimental or computational measurement of the scalings from which to determine whether the flow is, or is not, in the universality class of directed percolation, although model systems suggest that the transition is in this class (Barkley Reference Barkley2011, Reference Barkley2016; Shih, Hsieh & Goldenfeld Reference Shih, Hsieh and Goldenfeld2016). The scaling exponents depend on the spatial dimension of the system. Lemoult et al. (Reference Lemoult, Shi, Avila, Jalikop, Avila and Hof2016) recently carried out a study of Couette flow highly confined in two directions so that large-scale turbulent–laminar intermittency could manifest itself only along one spatial dimension. In both experiments and numerical simulations, they measured turbulence fraction as a function of Reynolds number and analysed the spatial and temporal correlations close to the critical Reynolds number. The results support a continuous variation of the turbulence fraction, from zero to non-zero at the onset of turbulence, with scaling laws consistent with the expectations for directed percolation in one spatial dimension.
In systems in which the flow is free to evolve in two large spatial directions, such as Couette and channel flow, the problem is much more difficult and the situation is less clear. Past work has suggested that the turbulence fraction varies discontinuously in plane Couette flow, and hence that transition in the flow is not of directed-percolation type (Bottin & Chaté Reference Bottin and Chaté1998; Bottin et al.
Reference Bottin, Daviaud, Manneville and Dauchot1998; Duguet, Schlatter & Henningson Reference Duguet, Schlatter and Henningson2010). More recently, Avila (Reference Avila2013) conducted experiments in a counter-rotating circular Couette geometry (radius ratio
$\unicode[STIX]{x1D702}=0.98$
) of large aspect ratio, and observed a variation of turbulence fraction with Reynolds number, suggesting a continuous transition to turbulence. Further investigation would be needed to determine whether the transition is in the universality class of directed percolation. Sano & Tamai (Reference Sano and Tamai2016) performed experiments on plane channel flow and concluded that this flow exhibits a continuous transition to turbulence in the universality class of directed percolation. However, they report a critical Reynolds number (based on the centreline velocity of the equivalent laminar flow) of
$Re=830$
, whereas other researchers (Xiong et al.
Reference Xiong, Tao, Chen and Brandt2015; Paranjape et al.
Reference Paranjape, Vasudevan, Duguet and Hof2017; Kanazawa, Shimizu & Kawahara Reference Kanazawa, Shimizu and Kawahara2017; T. Tsukahara & T. Ishida 2017, Private communication) observe sustained turbulent patches below 700. These later authors do not address the question of whether the transition is continuous or discontinuous, and further study is needed.
The goal of the present paper is threefold. Firstly, using a coupled-map lattice, we present the essential issues surrounding the onset of turbulence in a spatiotemporal setting, with particular emphasis on the case of two space dimensions. Secondly, we present a numerical study of a planar shear flow of unprecedented lateral extent and show that the onset of turbulence in this flow is continuous and is in the universality class of directed percolation. Finally, we discuss the issue of scales in the current and past studies, and we offer guidance to future investigations.
2 Coupled-map lattices and directed percolation revisited
Before discussing the planar shear flow, we revisit some important issues concerning spatiotemporal intermittency and directed percolation. The issues can most easily be illustrated using a coupled-map-lattice (CML) model. Such discrete-space, discrete-time models have been widely used to study the generic behaviour arising in spatially extended chaotic dynamical systems (e.g. Kaneko Reference Kaneko1985; Chaté & Manneville Reference Chaté and Manneville1988; Rolf, Bohr & Jensen Reference Rolf, Bohr and Jensen1998). Most notably they have been used as minimal models for describing the transition to turbulence in plane Couette flow (Bottin & Chaté Reference Bottin and Chaté1998).
The CML model is illustrated in figure 1. A state variable
$u$
is defined on a discrete square lattice (figure 1
a). Time evolution is given by discrete updates on the lattice. Specifically, letting
$u_{ij}$
denote the state variable at point
$(i,j)$
, the update rule for
$u$
is
$$\begin{eqnarray}u_{ij}\leftarrow f(u_{ij})+d\unicode[STIX]{x1D6E5}_{f}u_{ij},\end{eqnarray}$$
where the first term,
$f(u_{ij})$
, is the local dynamics and the second term is a nearest-neighbour diffusive-like coupling. Lattice sites are updated asynchronously by cycling through
$(i,j)$
in a random order for each step, as described by Rolf et al. (Reference Rolf, Bohr and Jensen1998). The control parameter is the coupling strength
$d$
.
The local dynamics are given by the map
$f$
shown in figure 1(b). It has a ‘turbulent’ tent region for
$0\leqslant u\leqslant 1$
and a ‘laminar’ region
$u>1$
surrounding the stable fixed point at
$u^{\ast }$
. In the absence of coupling, a turbulent site evolves chaotically and eventually makes a transition to the laminar state. Once a site becomes laminar, it will remain so indefinitely. The laminar fixed point is referred to as an absorbing state. Hence, the local dynamics are a simple caricature of a subcritical shear flow with coexisting turbulent and laminar flow states. Because the model is only a slight generalisation of those appearing in numerous past studies, we relegate the details to the appendix.
Figure 1. Intermittent transition in a coupled-map lattice. (a) Illustration of the lattice with nodes coloured according to whether the system is locally laminar (white) or turbulent (black). (b) The map defining the local dynamics at each node.
$u^{\ast }$
is a stable fixed point. (c) Typical time evolution, seen in a slice through the lattice at constant
$j$
, initialised with all sites in the turbulent state. (The spatial and temporal laminar gaps
$\ell _{x}$
and
$\ell _{t}$
are discussed in § 4.) (d) Equilibrium turbulence fraction
$F_{t}$
as a function of the coupling strength
$d$
in two cases:
$u^{\ast }=1.25$
(top) and
$u^{\ast }=1.1$
(bottom). In the top case the transition to turbulence is discontinuous while in the bottom case it is continuous. In the continuous case, close to the critical value
$d_{c}$
,
$F_{t}$
increases from zero with the universal power law for directed percolation in two space dimensions:
$F_{t}\sim (d-d_{c})^{\unicode[STIX]{x1D6FD}}$
, where
$\unicode[STIX]{x1D6FD}\simeq 0.583$
. The red curves in the main plot and inset show this power law.
We are primarily interested in the long-time dynamics of the system. We start from an initial condition with randomly selected values within the turbulent region. Figure 1(c) shows the evolution as seen in a one-dimensional slice through the lattice. The main quantity of interest is the turbulence fraction
$F_{t}$
, which is the fraction of sites in the turbulent state. After some time, the system will reach a statistical equilibrium and we can obtain the equilibrium value of
$F_{t}$
. If this is zero, then the system is everywhere in the laminar (absorbing) state. If it is non-zero, then at least some turbulence persists indefinitely.
A basic question is: how does the turbulence fraction at equilibrium depend on the coupling strength
$d$
, and in particular, how does it go from zero to non-zero? Figure 1(d) shows the two distinct cases: one discontinuous and one continuous. In the discontinuous case, there is a gap in the possible values of
$F_{t}$
. Long-lived transients with small turbulence fraction can be observed for values of
$d$
below
$d_{c}$
, the critical value of
$d$
, but the system simply cannot indefinitely maintain a small level of turbulence, no matter how large the system size. On the other hand, in the continuous case,
$F_{t}$
becomes arbitrarily small (in the limit of infinite system size) as
$d$
approaches
$d_{c}$
from above. In this case the system behaves in accordance with the power laws of directed percolation. In particular, as is shown, the turbulence fraction grows as
$F_{t}\sim (d-d_{c})^{\unicode[STIX]{x1D6FD}}$
, where
$\unicode[STIX]{x1D6FD}\simeq 0.583$
; see Lübeck (Reference Lübeck2004). We will discuss the other important power laws later when we analyse the planar fluid flow.
Note that on any finite lattice the minimum possible non-zero turbulence fraction is
$1/K$
(i.e. just one turbulent site), where
$K$
is the total number of lattice points. This means that even if the transition is continuous in principle, some discontinuity in the turbulence fraction from finite-size effects will be present in any numerical study. It is by investigating scaling behaviour, such as the log–log plot in figure 1(d), that one gains confidence in the nature of the transition.
2.1 Connection to turbulent transition and directed percolation
The difference between the continuous and discontinuous cases presented in figure 1(d) is only the location of the laminar fixed point
$u^{\ast }$
in the map
$f$
. Hence, either case could in principle correspond to a shear flow and there is no way to know a priori what type of transition could be expected. This point was well understood by the Saclay group in their early studies on transition in plane Couette flow (e.g. Bergé, Pomeau & Vidal Reference Bergé, Pomeau and Vidal1998; Bottin & Chaté Reference Bottin and Chaté1998; Bottin et al.
Reference Bottin, Daviaud, Manneville and Dauchot1998; Manneville Reference Manneville2016). Those experiments suggested a discontinuous transition to turbulence, based not only on the turbulence fraction, but also on the nature of transients below the critical point. Although we will argue that the physical size of those experiments was too small to produce a continuous transition, the conclusion reached was reasonable at the time.
More generally, directed percolation describes a stochastic process involving active and absorbing states (or equivalently bonds between sites that are randomly open or closed). As Manneville (Reference Manneville2016, § 4.2) notes, deterministic iterations of continuous variables coupled by diffusion will not necessarily behave in the same way as the directed-percolation process. The CML model presented here demonstrates this point. Depending on parameters, the system might, or might not, show a continuous transition in the universality class of directed percolation. Notwithstanding Pomeau’s conjecture, it is even less immediately evident that the full Navier–Stokes equations will behave in the same way, owing to the global nature of the pressure field for example. (For more technical details on absorbing state transitions, we refer the reader to Lübeck (Reference Lübeck2004) and references therein. One can find there details of the Janssen–Grassberger conjecture (Janssen Reference Janssen1981; Grassberger Reference Grassberger1982) concerning the ubiquity of the directed-percolation universality class.)
From hereon we shall use the notation of directed percolation and refer to the case of two space dimensions as (2+1)-D, meaning two spatial and one temporal dimension. In this notation, the spatial dimensions are referred to as perpendicular (
$\bot$
) and the temporal dimension as parallel (
$\Vert$
).
3 Waleffe flow
Pinning down the details of transition requires very large system sizes. For example, in the quasi-one-dimensional experiments of Lemoult et al. (Reference Lemoult, Shi, Avila, Jalikop, Avila and Hof2016), the long direction was more than
$2700$
times the fluid gap. Our goal is to achieve something approaching this size, but in two spatial directions and in a computational framework. To this end, we shall study a cousin of Couette flow, commonly referred to as Waleffe flow. This is the shear flow between parallel stress-free boundaries, driven by a sinusoidal body force. The two related computational advantages of this flow are that it lacks high-shear boundary layers near the walls and that the wall-normal dependence of the flow can be accurately represented by a few trigonometric functions. As shown in Chantry, Tuckerman & Barkley (Reference Chantry, Tuckerman and Barkley2016), a poloidal–toroidal representation with at most four trigonometric modes in the wall-normal direction,
$y$
, is capable of capturing turbulent bands and spots, the building blocks of turbulent–laminar intermittency. A Fourier representation is used for the large streamwise,
$x$
, and spanwise,
$z$
, directions.
In Chantry et al. (Reference Chantry, Tuckerman and Barkley2016), we showed that Waleffe flow corresponds closely to the interior of plane Couette flow, leading to a change in length scales from
$2h$
(the gap between walls in plane Couette flow) to
$1.25h$
(the Couette interior region) for Waleffe flow. Furthermore, this argument regarding the interior region leads to a comparable velocity scale
$U=1.6V$
, with
$V$
the maximum velocity of laminar Waleffe flow. The Reynolds number of the flow is then
$Re=Uh/\unicode[STIX]{x1D708}$
, where
$\unicode[STIX]{x1D708}$
is the kinematic viscosity. The sole change from Chantry et al. (Reference Chantry, Tuckerman and Barkley2016) is the addition of a small horizontal drag force
$-\unicode[STIX]{x1D70E}(u\boldsymbol{e}_{x}+w\boldsymbol{e}_{z})$
to the Navier–Stokes equation. Such a term, usually called Rayleigh or Ekman friction, is used in many hydrodynamic modelling contexts to approximate the effect of friction due to a solid boundary that has been omitted from the model. In geophysics (Marcus & Lee Reference Marcus and Lee1998; Pedlosky Reference Pedlosky2012, chap. 4) the inclusion of this term is the standard method of including the first-order departure from geostrophic flow due to the Ekman boundary layer between a stationary bottom and a rotating bulk. In their study of electromagnetically driven Kolmogorov flow in an electrolyte, Suri et al. (Reference Suri, Tithof, Mitchell, Grigoriev and Schatz2014) include such a term in their depth-averaged model of an assumed Poiseuille-like profile in order to account for the presence in their experiment of a solid boundary at the bottom of the fluid layer.
In our case, we introduce this force in order to damp flows with no curvature in
$y$
and very little curvature in
$x$
and
$z$
, which decay extremely slowly in Waleffe flow and which are not present at all in Couette flow. Our purpose is to use Waleffe flow to mimic the bulk region of Couette flow, which it does very well except for this point. The value
$\unicode[STIX]{x1D70E}=10^{-2}$
reproduces the damping to which these modes would be subjected in the wall regions of the corresponding Couette flow. In very large domains, without this damping the recovery of the laminar flow after spot decay is very slow. Beyond this, the damping has no effects on the phenomenology of Waleffe flow.
Figure 2. Intermittent turbulence typical of that found slightly above the onset of sustained turbulence. Visualised is streamwise velocity in the midplane at
$Re=173.824$
after
$1.2\times 10^{6}$
time units. Laminar flow is seen as white. The streamwise and spanwise size of the computational domain is
$2560h\times 2560h$
. The turbulence fraction is
$F_{t}\approx 0.1$
and the reduced Reynolds number is
$\unicode[STIX]{x1D716}=(Re-Re_{c})/Re_{c}=1.4\times 10^{-4}$
. The supplementary movie shows the time evolution in a domain of the same size for
$\unicode[STIX]{x1D716}=5.2\times 10^{-4}$
and
$\unicode[STIX]{x1D716}=5.7\times 10^{-5}$
, slightly above and below this value. For reference, the Couette domains of Bottin et al. (Reference Bottin, Daviaud, Manneville and Dauchot1998), Duguet et al. (Reference Duguet, Schlatter and Henningson2010), Avila (Reference Avila2013) and Lemoult et al. (Reference Lemoult, Shi, Avila, Jalikop, Avila and Hof2016) are overlaid in red, blue, orange and green respectively. The full streamwise length of the experiment of Lemoult et al. (Reference Lemoult, Shi, Avila, Jalikop, Avila and Hof2016) exceeds the figure size and is not fully shown. The spanwise width of the channel experiment of Sano & Tamai (Reference Sano and Tamai2016) is indicated in purple on the right.
We shall present results for domains of size [
$1280h$
,
$1.25h$
,
$1280h$
], [
$2560h$
,
$1.25h$
,
$2560h$
] and [
$5120h$
,
$1.25h$
,
$1280h$
]. Our largest square domain is plotted in figure 2, where we show a representative turbulent state slightly above the onset of turbulence. See also the supplementary movie available at https://doi.org/10.1017/jfm.2017.405, which shows the time evolution in a domain of this size. For this domain the highest resolved wavenumber in each horizontal direction is 2047, with
$3/2$
dealiasing used (leading to a grid spacing of 0.42). The same turbulence fractions are found in simulations with twice the resolution.
For context, the experiments of Bottin et al. (Reference Bottin, Daviaud, Manneville and Dauchot1998) used a domain of size
$[380h,2h,70h]$
, Prigent et al. (Reference Prigent, Grégoire, Chaté and Dauchot2003) used a domain of size
$[770h,2h,340h]$
and Avila (Reference Avila2013) used a domain of size
$[622h,2h,526h]$
. To date, the largest simulations have been those of Duguet et al. (Reference Duguet, Schlatter and Henningson2010), who considered a domain of size
$[800h,2h,356h]$
. Both Bottin et al. (Reference Bottin, Daviaud, Manneville and Dauchot1998) and Duguet et al. (Reference Duguet, Schlatter and Henningson2010) report evidence of a discontinuous transition, unable to sustain turbulence fractions significantly below 0.4, while Avila (Reference Avila2013) observed evidence of a continuous transition, with sustained turbulence fractions as small as
$0.07$
.
The flow at
$(x,z,t)$
is defined as turbulent if
$E(x,z,t)>E_{T}$
, where
$E(x,z,t)$
is the
$y$
-integrated energy of the velocity deviation from the laminar state and
$E_{T}=0.01$
is a threshold. Varying
$E_{T}$
between 0.001 and 0.05 changes only slightly the size of patches deemed turbulent and has no effect on any of the scaling relationships to follow.
4 Results
Figure 3. (a) Turbulence fraction as a function of time for a range of Reynolds numbers with an initial condition of uniform turbulence. Above criticality, the turbulence fraction saturates at a finite value, and below it falls to zero. At criticality, the turbulence fraction decays in time as a power law
$F_{t}\sim t^{-\unicode[STIX]{x1D6FC}}$
with the (2+1)-D directed-percolation exponent
$\unicode[STIX]{x1D6FC}\simeq 0.4505$
(dashed line). Coloured lines for decreasing turbulence fractions correspond to Reynolds numbers
$[173.952,173.888,173.840,173.824,173.792,173.773,173.696,173.568]$
. The supplementary movie shows the decay of the turbulence fraction at Reynolds numbers 173.888 and 173.808. (b) Data above and below criticality collapse onto two scalings (black dashed curves) when the directed-percolation exponents are used to rescale time and turbulence fraction.
In figure 3(a), we plot the time evolution of the turbulence fraction
$F_{t}$
for a series of Reynolds numbers. Each run was initialised from uniform turbulence and run until a saturated turbulence fraction was reached (quench protocol, see Bottin & Chaté (Reference Bottin and Chaté1998)). Below a critical value
$Re_{c}=173.80$
(to five significant figures), the turbulence fraction eventually falls off to zero, while above
$Re_{c}$
it saturates at a finite value. (
$Re_{c}$
for this system differs from that of plane Couette flow.)
In figure 4(a) we plot the equilibrium turbulence fraction as a function of
$Re$
. We find clear evidence for a continuous transition in which small
$F_{t}$
can be sustained given a sufficiently large domain. The saturated turbulence fraction follows a power law
$F_{t}\sim \unicode[STIX]{x1D716}^{\unicode[STIX]{x1D6FD}}$
where
$\unicode[STIX]{x1D716}\equiv (Re-Re_{c})/Re_{c}$
.
$Re_{c}$
is determined as the value of
$Re$
that minimises the mean squared error of a linear fit of the logarithms of
$\unicode[STIX]{x1D716}$
and
$F_{t}$
(see figure 4
b). This linear fit estimates
$\unicode[STIX]{x1D6FD}=0.58\pm 0.04$
with a
$95\%$
confidence interval. This agrees with the (2+1)-D directed-percolation value
$\unicode[STIX]{x1D6FD}\simeq 0.583$
(dashed line). Figure 4(d) shows the turbulence fraction obtained from our system in a domain whose size is that of the experiments of Bottin & Chaté (Reference Bottin and Chaté1998). (See also figure 2.) As was observed experimentally, below a turbulence fraction of about
$0.5$
, turbulence appears only as a long-lived transient state, and hence the equilibrium turbulence fraction exhibits a discontinuous transition. This strongly suggests that the discontinuous transitions reported for plane Couette flow (Bottin & Chaté Reference Bottin and Chaté1998; Bottin et al.
Reference Bottin, Daviaud, Manneville and Dauchot1998; Duguet et al.
Reference Duguet, Schlatter and Henningson2010) are due to finite-size effects. Interestingly, long-lived transient states in the small system have turbulence fractions close to those for equilibrium states in our large domain; they just are not sustained states. Also note that the critical Reynolds number is not greatly affected by system size.
Figure 4. Bifurcation diagrams for the transition to turbulence. (a) Continuous transition in a large domain:
$[2560h,1.25h,2560h]$
. Equilibrium turbulence fraction
$F_{t}$
is plotted as a function of
$Re$
. Points and error bars denote mean and standard deviation of
$F_{t}$
. The black dashed curved shows the directed-percolation power law. (b) Log–log plot of the same data in terms of
$\unicode[STIX]{x1D716}=(Re-Re_{c})/Re_{c}$
, where
$Re_{c}=173.80$
. Near criticality the data are consistent with
$F_{t}\sim \unicode[STIX]{x1D716}^{\unicode[STIX]{x1D6FD}}$
with
$\unicode[STIX]{x1D6FD}\simeq 0.583$
. (c)
$F^{1/\unicode[STIX]{x1D6FD}}$
against
$Re$
showing linear behaviour. (d) Discontinuous transition in a domain of size
$[380h,1.25h,70h]$
, approximately that of the experiments by Bottin & Chaté (Reference Bottin and Chaté1998). (See figure 2.) Filled points denote sustained turbulence, while open points denote the turbulence fraction of long-lived transient turbulence. The dashed curve is the directed-percolation power law from the large domain.
Figure 5. Exponents for temporal and spatial correlations. (a–c) Distributions of laminar gaps in (a) time and the two spatial directions, (b)
$x$
and (c)
$z$
, for a variety of domain sizes at
$Re=173.824$
(
$\unicode[STIX]{x1D716}=1.4\times 10^{-4}$
).
$N$
is the gap count, normalised by the shortest gap count. The directed-percolation scalings (
$\unicode[STIX]{x1D707}_{\Vert }\simeq 1.5495$
and
$\unicode[STIX]{x1D707}_{\bot }\simeq 1.204$
) are plotted as dashed lines and show excellent agreement with the
$t$
- and
$z$
-gaps. For the
$x$
-gaps a power law closer to
$-1$
is observed. (d–f) Exponential tails of the gap distributions for several values of
$\unicode[STIX]{x1D716}$
just above criticality. Increasing domain sizes are used for points closer to criticality. Insets show scaling of correlation lengths
$\unicode[STIX]{x1D709}$
with
$\unicode[STIX]{x1D716}$
together with the directed-percolation exponents:
$\unicode[STIX]{x1D708}_{\Vert }\simeq 1.295$
and
$\unicode[STIX]{x1D708}_{\bot }\simeq 0.733$
. (g–i) Collapse of the data using directed-percolation power laws
$\unicode[STIX]{x1D707}$
and
$\unicode[STIX]{x1D708}$
(see axis labels/text).
To substantiate whether a system is in the directed-percolation universality class, it is necessary to verify three independent power-law scalings close to criticality; see Takeuchi et al. (Reference Takeuchi, Kuroda, Chaté and Sano2009), whose approach we will follow closely. Having demonstrated the scaling of
$F_{t}$
(exponent
$\unicode[STIX]{x1D6FD}$
), we now turn to scalings associated with temporal and spatial correlations.
One approach to determining the correlations is via the distribution of laminar gaps at
$Re\simeq Re_{c}$
(figure 5
a–c). The flow has a temporal laminar gap of length
$\ell _{t}$
if
$E(x,z,t)>E_{T}$
and
$E(x,z,t+\ell _{t})>E_{T}$
but
$E(x,z,t^{\prime })<E_{T}$
for
$0<t^{\prime }<\ell _{t}$
. Spatial gaps
$\ell _{x}$
and
$\ell _{z}$
are defined similarly. Such gaps are illustrated for the CML model in figure 1(c). From simulations just above
$Re_{c}$
, we generate gap distributions by measuring and binning the laminar gaps within the intermittent flow once the turbulence fraction has saturated. Given the anisotropy between the streamwise and spanwise directions, we measure gaps in these directions separately. At criticality, a system within the directed-percolation universality class displays power-law behaviour,
$N\sim \ell ^{-\unicode[STIX]{x1D707}}$
, where
$N$
is the number of gaps of length
$\ell$
. The temporal gaps (figure 5
a) show excellent scaling with the directed-percolation temporal exponent
$\unicode[STIX]{x1D707}_{\Vert }\simeq 1.5495$
. This is also true of the spanwise gaps (figure 5
c) with the spatial exponent
$\unicode[STIX]{x1D707}_{\bot }\simeq 1.204$
. However, the streamwise laminar gaps (figure 5
b) do not show a clear, extended power law, and to the extent that there is a power law, the exponent is closer to 1 than to
$\unicode[STIX]{x1D707}_{\bot }\simeq 1.204$
. Indeed, Takeuchi et al. (Reference Takeuchi, Kuroda, Chaté and Sano2007) observed that the laminar gap distribution in one direction of a liquid crystal layer had an exponent closer to 1 than to
$\unicode[STIX]{x1D707}_{\bot }\simeq 1.204$
, This is also true for our simulations of a non-isotropic CML very slightly above the critical point, thus indicating that the issue here is that the flow is not exactly at
$Re_{c}$
, as it should be for the scaling to hold. Although the gap distribution in each spatial direction should show power-law behaviour, these may converge at different rates as
$Re\rightarrow Re_{c}$
.
Given the poor agreement for the streamwise gaps, we use a second approach to measure the percolation exponents, which does not rely on simulations at
$Re_{c}$
. Away from criticality, power-law behaviour will be seen only over a finite range of temporal and spatial gap lengths. Beyond these lengths, exponential tails are expected of the form
$N\sim \exp (-\ell /\unicode[STIX]{x1D709})$
, with correlation lengths
$\unicode[STIX]{x1D709}$
diverging as
$\unicode[STIX]{x1D716}$
goes to zero:
$\unicode[STIX]{x1D709}\sim \unicode[STIX]{x1D716}^{-\unicode[STIX]{x1D708}}$
. In figure 5(d–f) we fit exponential tails for several values of
$\unicode[STIX]{x1D716}$
. In the insets we plot
$\unicode[STIX]{x1D709}$
as a function of
$\unicode[STIX]{x1D716}$
and compare with the expected exponents for directed percolation. The exponents
$\unicode[STIX]{x1D707}$
and
$\unicode[STIX]{x1D708}$
are exactly related via
$\unicode[STIX]{x1D707}=2-\unicode[STIX]{x1D6FD}/\unicode[STIX]{x1D708}$
, thus giving
$\unicode[STIX]{x1D708}_{\Vert }\simeq 1.295$
and
$\unicode[STIX]{x1D708}_{\bot }\simeq 0.733$
(Lübeck Reference Lübeck2004). Because the exponents
$\unicode[STIX]{x1D707}$
and
$\unicode[STIX]{x1D708}$
are linked, the power laws in the middle row of figure 5 are not independent of the corresponding power laws in the top row. However, they rely on different data and hence are not limited by the same finite-size and finite-distance-from-critical effects. The power law for the
$x$
direction is now seen to be in clear agreement with the directed-percolation exponent, as are those in
$t$
and
$z$
.
Combining these power laws, the laminar gap distributions can be collapsed using the relationship
$N\ell ^{\unicode[STIX]{x1D707}}=G(\ell \unicode[STIX]{x1D716}^{\unicode[STIX]{x1D708}})$
, where
$G$
is an unknown function (Henkel et al.
Reference Henkel, Hinrichsen, Lübeck and Pleimling2008, pp. 111–112). In figure 5(g–i) we plot our data in collapsing coordinates using the (2+1)-D percolation exponents. This collapse is well illustrated by the temporal gap distribution, a culmination of the excellent fits of
$\unicode[STIX]{x1D707}_{\Vert }$
and
$\unicode[STIX]{x1D708}_{\Vert }$
. For the
$x$
-gaps, only gaps of length 100 or greater are counted, corresponding to the start of the power law in figure 5(b). The collapse of the
$z$
-gaps is hindered by the
$\unicode[STIX]{x1D708}_{\bot }$
scaling seen in figure 5(f), but close to criticality (last three lines) the data begin to show collapse.
We have also run simulations in a quasi-1D, streamwise-oriented domain similar in spirit to the experiments of Lemoult et al. (Reference Lemoult, Shi, Avila, Jalikop, Avila and Hof2016) (see the experimental domain in figure 2). In a domain of size
$[1280h,1.25h,40h]$
, the distribution of streamwise laminar gaps near criticality exhibits a clear power law, in contrast to the poor power-law behaviour found for streamwise laminar gaps in the full planar system (figure 5
b). The exponent is
$\unicode[STIX]{x1D707}_{\bot }\simeq 1.748$
, as predicted for systems with a single spatial dimension ((1+1)-D directed percolation).
We return to the time evolution shown in figure 3(a). Between the evolution at
$Re<Re_{c}$
and
$Re>Re_{c}$
is the power-law decay predicted for a directed-percolation process at criticality:
$F_{t}\sim t^{-\unicode[STIX]{x1D6FC}}$
, where
$\unicode[STIX]{x1D6FC}=2-\unicode[STIX]{x1D707}_{\Vert }\simeq 0.4505$
. Close to criticality, we observe evidence for this power law in the data. This plot highlights a major challenge to simulations near criticality – well over
$10^{5}$
time units are required to reach even the moderately small turbulence fractions simulated here. As was noted by Avila (Reference Avila2013, p. 32), these long time scales proved an issue in the work of Duguet et al. (Reference Duguet, Schlatter and Henningson2010), who in
$10^{4}$
time units of simulation were unable to converge turbulence fractions much below 0.4. As in the present study, Avila (Reference Avila2013, p. 90) let the system evolve for
$O(10^{6})$
advective time units close to transition. Hence both simulation time and domain size can be limiting factors in observing the hallmarks of percolation. Using directed-percolation scalings (Takeuchi et al.
Reference Takeuchi, Kuroda, Chaté and Sano2009), the data above and below criticality collapse onto two curves (figure 3
b), highlighting the universality of directed percolation in the transition to turbulence.
5 Discussion
Over the years several attempts have been made to quantify the transition to turbulence and to determine whether or not subcritical shear flows follow the spatiotemporal scenario of directed percolation. Such attempts have consistently been frustrated by the large system sizes required to address the issue. Here we have performed simulations of a planar flow of sufficient size that we have been able to eliminate significant finite-size, finite-time effects, and thereby to examine in full detail the onset of turbulence in a planar example. We have demonstrated both that the equilibrium turbulence fraction increases continuously from zero above a critical Reynolds number and that statistics of the turbulent structures exhibit the power-law scalings of the directed-percolation universality class. Meeting such demands has necessitated not only turning to the stress-free boundaries of Waleffe flow, but further truncating the simulations to just four wall-normal modes. Performing a comparable computational study directly on plane Couette flow is currently far beyond available resources.
In light of what we now understand about the scales needed to capture sparse turbulent structures near the onset of turbulence, we have re-examined the apparent discontinuous transition to turbulence reported in past studies of plane Couette flow. The conclusions of those studies were reasonable at that time, but our results indicate that prior experimental system sizes were too small, and prior simulation times were too short, to accurately capture sustained turbulence close to onset – both space and time constraints can limit estimates of true equilibrium dynamics. Apart from overall issues of scale, we have observed that the scaling relations in the streamwise and spanwise directions may converge at different rates. Because shear flows are non-isotropic, it is important to monitor these directions separately. These considerations should guide the design of experiments and computations. In this regard, it should be noted that while the Reynolds numbers in our study differ from those of plane Couette flow, the length and time scales are closely comparable to those of plane Couette flow. The efficiency with which our system can be simulated offers potential for use in conjunction with future investigation.
While we cannot rule out the possibility that other subcritical shear flows follow some different route to turbulence, we know that truncated Waleffe flow contains the essential self-sustaining mechanism of wall-bounded turbulence and that it produces the oblique turbulent bands that characterise transitional turbulence in plane Couette and plane channel flow (Chantry et al. Reference Chantry, Tuckerman and Barkley2016). The closeness of these phenomena suggests that all of these flows exhibit the same route to turbulence.
Acknowledgements
We thank Y. Duguet and P. Manneville for useful discussions. We are grateful to A. Lemaître for advice on the CML model and specifically for suggesting the use of
$u^{\ast }$
as a parameter to vary between discontinuous and continuous transitions. M.C. was supported by the grant TRANSFLOW, provided by the Agence Nationale de la Recherche (ANR). This research was supported in part by the National Science Foundation under grant no. NSF PHY11-25915. This work was performed using high performance computing resources provided by the Institut du Developpement et des Ressources en Informatique Scientifique (IDRIS) of the Centre National de la Recherche Scientifique (CNRS), coordinated by GENCI (Grand Équipement National de Calcul Intensif).
Supplementary movie
A supplementary movie is available at https://doi.org/10.1017/jfm.2017.405.
Appendix A. Details of the CML
We provide here details of the CML model and simulations shown in § 2. For the most part, these follow previous works (e.g. Bottin & Chaté Reference Bottin and Chaté1998; Rolf et al.
Reference Rolf, Bohr and Jensen1998). The local dynamics is given by the map
$f$
:
$$\begin{eqnarray}f(u)=\left\{\begin{array}{@{}ll@{}}ru,\quad & u\leqslant 1/2\\ r(1-u),\quad & 1/2<u\leqslant 1\\ k(u-u^{\ast })+u^{\ast },\quad & 1<u,\end{array}\right.\end{eqnarray}$$
where
$r$
,
$k$
and
$u^{\ast }$
are parameters. Here we fix
$r=3.0$
and
$k=0.8$
. The only difference between this map and those used previously is that here
$u^{\ast }$
is a free parameter, rather than being set by the value of
$r$
to
$u^{\ast }=(r+2)/4$
.
The spatial coupling is given by
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6E5}_{f}u_{ij} & = & \displaystyle \frac{1-\unicode[STIX]{x1D6FF}}{4}(f(u_{i-1,j})-2f(u_{i,j})+f(u_{i+1,j}))\nonumber\\ \displaystyle & & \displaystyle +\,\frac{1+\unicode[STIX]{x1D6FF}}{4}(f(u_{i,j-1})-2f(u_{i,j})+f(u_{i,j+1}))\end{eqnarray}$$
subject to periodic boundary conditions. This term differs from the standard coupling only in that the parameter
$\unicode[STIX]{x1D6FF}$
permits different coupling strengths in the
$i$
and
$j$
directions, to mimic the anisotropy of planar shear flows. We use
$\unicode[STIX]{x1D6FF}=0.6$
. This anisotropy has no significance for the results presented in this paper since continuous and discontinuous transitions occur also in the isotropic case
$\unicode[STIX]{x1D6FF}=0$
. We show results at
$\unicode[STIX]{x1D6FF}=0.6$
only for consistency with future publications.
