This study provides a numerical characterization of the basin of attraction of the laminar Hagen–Poiseuille flow by measuring the minimal amplitude of a perturbation required to trigger transition. For pressure-driven pipe flow, the analysis presented here covers autonomous and impulsive scenarios where either the flow is perturbed with an initial disturbance with a well-defined norm or perturbed by means of local impulsive forcing that mimics injections through the pipe wall. In both the cases, the exploration is carried out for a wide range of Reynolds numbers by means of a computational method that numerically resolves the transitional dynamics. For , the present work provides critical amplitudes that decay as Re−3/2 and Re−1 for the autonomous and impulsive scenarios, respectively. For Re=2875, accurate threshold amplitudes are found for constant mass-flux pipe by means of a shooting method that provides critical trajectories that never relaminarize or trigger transition. These transient states are used as initial guesses in a damped Newton–Krylov method formulated to find periodic travelling wave solutions that either travel downstream or exhibit a helicoidal advection.
Since the original study of Reynolds (1883), subcritical transition in pipe flow has been an object of analysis by mathematicians, physicists and engineers. Flows in experimental pipes become naturally turbulent for moderately high flow rates, despite the linear stability of the basic flow (Priymak & Miyazaki 1998; Meseguer & Trefethen 2003). As usual the stability of the flow is governed by the Reynolds number, defined as , where is the cross-sectional mean axial speed, a is the pipe radius and ν is the kinematic viscosity of the fluid. There is a long list of experimental, theoretical and numerical studies regarding the stability properties of this flow. We refer the reader to recently published reviews (Kerswell 2005; Eckhardt et al. 2007) and references therein.
Since the Hagen–Poiseuille flow is linearly stable, transition in pipes must be explained in terms of other mechanisms. The initial stage of transition has been mainly ascribed to the linear non-modal transient growth that small perturbations experience in shear flows (Schmid & Henningson 1994a,b; Mellibovsky & Meseguer 2006). Such a transition mechanism is fast and explosive, as reported in many experimental results (Hof et al. 2003; Peixinho & Mullin 2006, 2007). However, linear mechanisms by themselves cannot explain how turbulence sets in. The final stage of transition has been proved to be dominated by the so-called self-sustained process (SSP), consisting of streamwise vortices generating streaks that are eventually destabilized by streamwise-dependent waves capable of regenerating the original vortices via nonlinear selection rules (Waleffe 1997). The use of suitable homotopy transformations in the Navier–Stokes equations emulating the SSP mechanism has been crucial in determining new finite amplitude solutions in shear flows (Waleffe 2003). Several families of unstable travelling waves have been recently identified numerically either using the underlying SSP mechanism (Wedin & Kerswell 2004) or mimicking the key elements of the Nagata–Busse–Clever states (Nagata 1990) in plane Couette flow (Faisst & Eckhardt 2003). Whereas these travelling waves have been shown to play a relevant role in organizing turbulent dynamics (Hof et al. 2004), their involvement in the transition process has not been clarified.
During the last decade, one of the main goals in pipe flow stability research has been to provide a characterization of the basin of attraction of the basic laminar flow. We must think of it as a subset in an infinite dimensional space that contains the basic flow, driving towards this solution any initial perturbation contained in this subset. A still unsolved question is the dependence of the critical amplitude Ac or size of the boundary of that basin of attraction with the Reynolds number. Provided that the minimum amplitude required to trigger transition decreases when Re is increased, it is plausible to assume that its asymptotic behaviour scales with Re according to(1.1)with γ necessarily negative. Asymptotic exponents for plane channel flows have been recently provided within the framework of some particular transition scenarios (Chapman 2002). For pipe flow, recent renormalizations (Trefethen et al. 2000) have been suggested in order to cast different experimental results in terms of a single definition of the amplitude appearing in (1.1), providing lower and upper bounds for the value of this critical exponent that presumably lies within the interval γ∈[−9/5,−6/5]. Recent experiments that have been carried out in Hof et al. (2003), henceforth referred as Hof, Juel and Mullin (HJM), clearly provided a critical exponent γ=−1, whereas the first numerical estimation of the threshold exponent problem (1.1) in pipe flow advanced a distinct value of γ=−3/2 (Meseguer 2003).
The apparent simplicity of expression (1.1) hides many aspects that require an accurate description. First, a mathematical definition of the amplitude A appearing in (1.1) must be provided. Second, depending on the way the fluid is advected downstream (either by a pressure-driven mechanism or at a constant mass-flux rate using a piston) the perturbations may evolve quite differently and particularly for moderate flow rates. Third, the perturbation may develop from an initial disturbance added to the basic flow, the fluid system evolving in an autonomous fashion, or it may develop from an impulsive time-dependent source such as fluid injections from the pipe wall. Fourth, when studying the time evolution of a perturbation in an open flow, advection is crucial, since potential turbulent transients are advected downstream and leave the domain, making it impossible to classify the dynamics for long times. Fifth and last, expression (1.1) is meaningful only for high values of Re.
Even recent experimental results seem to point to a transitory nature of turbulence (Hof et al. 2006), the notion of a threshold separating initial conditions that trigger transition, however short-lived, from others that uneventfully decay still applies (Nusse & Yorke 1989), leading to the question of what a solution wandering about criticality would look like. These so-called critical trajectories have recently been shown to approach a chaotic saddle: the edge state (Schneider & Eckhardt 2006; Eckhardt et al. 2007). The edge state resides in the non-empty boundary separating the laminar and turbulent basin of attraction, i.e. the critical threshold, and therefore seems to govern transition.
We consider the motion of an incompressible viscous fluid of kinematic viscosity ν and density ρ driven through a circular pipe of radius a and length L by means of a uniform axial pressure gradient π0. As usual, the problem is rendered dimensionless using the pipe radius and the maximum Hagen–Poiseuille flow speed at the centre line, , as length and velocity units, respectively. In cylindrical non-dimensional coordinates (r, θ, z), the Hagen–Poiseuille flow reads . Accordingly, the computational domain considered is , with Λ=L/a, i.e. the aspect ratio of the pipe, in radii units. Under the presence of disturbances, the flow in the pipe v is usually expressed as ub+u, where u(r, θ, z, t) is a time-dependent solenoidal disturbance satisfying radial homogeneous boundary conditions, u(1, θ, z, t)=0 and axial–azimuthal periodicity(2.1)Formal substitution of the flow in the Navier–Stokes equations leads to a nonlinear partial differential equation for u and the pressure perturbation q(2.2)with(2.3)and . A highly accurate solenoidal spectral method has been used to discretize (2.2), where the spectral approximation identically satisfies (2.3). As usual in numerical simulation of open flows, periodic boundary conditions are imposed at the ends of the pipe so the Navier–Stokes solver is based on Fourier expansions in (θ, z) and Chebyshev polynomials in r, combined with a fourth-order linearly implicit time-marching scheme (Meseguer & Mellibovsky 2007).
The normalized energy of a perturbation u is measured by means of the volume integral or hermitic product(2.4)so that the amplitude of the perturbation is defined as the square root of its normalized energy(2.5)To better understand how the energy is distributed, it is very convenient to express the flow in the pipe v as a sum of the basic flow ub plus the Fourier components of the perturbation field u(2.6)where u00 contains the azimuthal–axial-averaged perturbation velocity profile, u2D represents the non-axisymmetric streamwise component of the velocity field and u3D represents the remaining streamwise-dependent components. Using (2.4) on the decomposed velocity field, the energies corresponding to the bulk flow, to the streamwise component and to the three-dimensional perturbation can be computed independently as , and , respectively.
3. Critical threshold and exponents for high Re
In this section, critical exponents resulting from two different numerical explorations are provided. In the first case, the disturbances come from an initial array of streamwise vortices that trigger inflectional profiles that are unstable with respect to streamwise-dependent waves. In the second case, the flow is perturbed by means of a localized impulsive forcing that mimics recent experimental studies (Hof et al. 2003).
(a) Autonomous scenario
Former numerical explorations confirm that an efficient way to trigger instability in pipe flow is by adding streamwise vortical perturbations to the basic flow (Zikanov 1996; Meseguer 2003). In particular, streamwise vortices with azimuthal wavenumber nv=1 are very efficient in generating optimal transient growth (Schmid & Henningson 1994b; Meseguer & Trefethen 2003). This energy growth leads to the generation of streaks (inflectional profiles with strong variations of streamwise flow speed) eventually destabilized by streamwise-dependent modes of selected axial periodicity (Zikanov 1996). Two and three pairs of vortices experience lower transient growth, but may provide streaks with higher capability of destabilizing three-dimensional waves. The initial disturbance used in the exploration of the critical threshold consists of a suitable superposition of nv=1, 2, 3-pairs of two-dimensional streamwise vortices (henceforth, N1, N2 and N3 disturbances, respectively), , with random three-dimensional noise added on top(3.1)where c.c. stands for complex conjugated terms. The radial structure of takes the simplest solenoidal form capable of generating streaks(3.2)where σ=1 (2) for nv odd (even), whereas contains a random perturbation (urand) velocity field. The complex constants C2D and C3D in (3.1), which modulate the initial amplitude of the two components of the perturbation, are chosen so that the initial energy of the streamwise vortices, , and of the three-dimensional modes, , take the desired values and , respectively, satisfying .
In order to illustrate the effects of the perturbations described above, figure 1 shows three snapshots of the streak breakdown process for the vortical disturbance N3. In figure 1a, z-averaged cross-sectional contours of the axial speed component of the flow are represented within the range in order to visualize the streaks formation and destabilization. Figure 1b,c correspond to energy density contours of the velocity component appearing in (2.6). More specifically, figure 1b contains z-averaged cross-sectional contours, while figure 1c shows θ-averaged contours on a transversal section . The aforementioned energy density averages are given by(3.3)and(3.4)and their contours are drawn in arbitrary units. In addition, the axial coordinate of the longitudinal sections has been conveniently scaled to aid representation. These series of contours reveal the modal structure of the three-dimensional waves as well as their location and destabilization effects over the streaks. The three-dimensional perturbation organizes itself and grows exponentially in the vicinity of the saddle lines of the streaks-modulated axial velocity profiles, as can be seen in the snapshots at t=100. Once the three-dimensional perturbation has reached a sufficient energy level, nonlinear interaction with the streaks starts (200≤t≤220), destabilizing the laminar profile and leading to turbulence.
In order to establish criteria to decide whether the perturbations (3.1) lead to turbulence or not, it is crucial to run up to a time horizon at which the streaks have fully developed and the three-dimensional perturbations have had enough time to grow. This time has been found to be at least T=1000 advective time units for the lowest ϵ2D at the highest Re explored. After this period, either the streak breakdown or the irreversible onset of viscous decay had taken place. In the present study, a simulation run is considered turbulent if(3.5)otherwise laminar. Condition (3.5) is based on the fact that three-dimensionality is a clear signature of turbulent dynamics and therefore it is required for the streamwise-dependent modes to be still active, and much stronger than initially, at the end of the run.
The critical amplitude threshold exploration covers a wide range of Reynolds numbers, within the interval Re∈[2500,12 600], with a highly resolved spectral grid of and with Λ=20. Results are shown in figure 2. As expected, the critical amplitude Ac is a decreasing function of the Reynolds number. In fact, Ac exhibits a vertical threshold evidenced by the behaviour of the slope that becomes very pronounced at low Re (allegedly converging to a vertical asymptote at Recr≲2000). As soon as Re is increased, the numerical results shown in figure 2 clearly reveal a critical amplitude that decreases according to for N1, for N2 and for N3 (dashed straight lines in figure 2), at least within the studied range. Finally, an enhanced exploration with N1 vortical perturbations but only exciting optimally destabilized three-dimensional modes (Zikanov 1996) produced a further optimized critical threshold that scaled as . No remarkable differences were found when increasing the spectral resolution or the aspect ratio Λ (Mellibovsky & Meseguer 2006).
(b) Impulsive scenario
Very recent experiments that have been carried out in HJM (Hof et al. 2003) explored transition phenomena of pipe flow subjected to finite amplitude impulsive perturbations for a wide range of axial speeds of the flow. The experiments reported in HJM were carried out in a long aspect ratio pipe, with a piston that kept the mass-flux constant during every run and where the disturbances were generated by impulsively injecting fluid into the main flow through six slits azimuthally equispaced on a perimeter around the pipe located at a fixed axial position far downstream from the pipe inlet, so that the base flow was sufficiently developed. Figure 3a is a schematic plot of the injection device used in HJM, where the injected fluid jet penetrates into the basic flow with an angle φ=π/3 with respect to the radial coordinate, in a plane normal to the pipe axis. The injection is activated following a step-like time-dependent function, active for a prescribed injection duration (Hof et al. 2003). The experiments clearly concluded that the minimum amplitude of a perturbation required to trigger transition scaled as Re−1. The experimental procedure of perturbing the basic flow would correspond to the category of an impulsive perturbative scenario. The main goal of this section is to gain some insight on the internal mechanisms responsible for transition by reproducing numerically the experiments of HJM.
We study the perturbative effects generated by introducing in (2.2) an impulsive volume forcing term, f, that acts locally in time and space as an accelerator of the fluid, reproducing the effects of an injection. The equation of the perturbation now reads(3.6)where the forcing field has the following structure:(3.7)with fa a constant amplitude factor and , i.e. a Heaviside double-step function that activates the injection within the lapse . Finally, fs provides the spatial structure of the six-jet injection whose explicit mathematical expression can be found in Mellibovsky & Meseguer (2007). Figure 3b shows the acceleration field inflicted by the modelled forcing upon the stationary basic flow.
The injection amplitude defined in HJM is given by the ratio between the total massflow injected through the Ninj slits and the pipe mass flow upon injection(3.8)Assuming that the injected jets can be described by the axisymmetric jet theory (Schlichting 1968), we can obtain an equivalent expression of the amplitude defined in terms of the asymptotic maximum speed of the jet u measured at a distance g from the slit and the cross-sectional area of the slit Sinj(3.9)with , and (Mellibovsky & Meseguer 2007). To the authors' knowledge, equation (3.9) provides the first quantitative approximation to a law relating an injection property that is mensurable in under-resolved computations and the experimental amplitude.
In HJM, it was observed that the critical amplitude was not altered for injection durations of advective time units, where Δt0 will be considered later as a reference time interval for amplitude renormalization purposes. This phenomenon is reproduced by our numerical simulations. Figure 4a shows the flow speed ug(t) measured at the point of maximum forcing norm for different critical runs carried out for Δtinj=2, 3, 4, 5, 8 and 16 time units. On each of these curves, the forcing was stopped at the indicated instants of time (grey circles), eventually leading to transition for longer times (not shown). Figure 4a also shows the behaviour of the flow speed for the same runs, but with permanent forcing (dashed curves). It can be observed that there is a clear stagnation of ug to a constant value u* in some cases, it being easy to obtain a sharp asymptotic value u* to be used afterwards in relation (3.9) to identify its corresponding amplitude . A consistent comparison between numerics and experiments is carried out by defining the normalized amplitude of a Δtinj-lapse perturbation as(3.10)Therefore, the experimental and numerical threshold amplitudes are normalized independently with respect to their corresponding reference saturation values for Δt0=24, thus rendering the critical amplitudes independent of the injection geometrical features appearing in (3.9). Figure 4b shows the normalized threshold amplitudes Aexp and Anum for different injection lapses.
The critical amplitude threshold has been systematically tracked for Reynolds numbers in the range Re∈[2512,14 125] and the injection duration held fixed to Δtinj=20 in all the explorations. The computations were carried out using a spectral grid within a pipe of Λ=40 radius.
In what follows, geometrical discrepancies between experiments and the actually modelled injection are removed by normalizing the amplitudes according to(3.11)where Re is the Reynolds number of the basic flow before being perturbed and Re0=14 000. The critical amplitude threshold results are shown in figure 5. The amplitude A, normalized according to (3.11), has been plotted as a function of Re for both experiments (grey squares) and computations (white circles), along with a dashed line indicating a slope of γ=−1. It is remarkable how experiments and computations exhibit very similar behaviour at high Re, which seems to show that the numerical model properly captures the transition mechanisms observed in the laboratory. Doubled axial-resolution check runs at high Re have also been plotted (black dots), with no significant change in the critical threshold slope, which reassures us that the resolution chosen for the bulk of the computations is sufficient. At the low Re-range, however, while experiments exhibit the same characteristic asymptotic behaviour from Re values as low as 2000, the computations fail to do so. In fact, the numerical simulations seem to find a vertical stability threshold for Re∼2500, at least for the type of perturbations used. This discrepancy can be mainly ascribed to the constant mass flow and pressure-driven differing natures of the pipes. In our pressure-driven computational pipe, the actual mass flow has a tendency to drop dramatically as soon as the perturbation grows and reorganizes the flow, particularly because of the short length of the domain where intermittency phenomena may fill a high portion of it. Conversely, in the constant mass flow experiments of HJM, some energy may be restituted into the flow through the action of the constant speed piston as soon as any perturbation slightly alters the mass flow. The discrepancy, however, is mitigated by an increase in the Reynolds number, for which subcritical perturbations do not alter the mass flow considerably (Mellibovsky & Meseguer 2007).
4. Critical threshold for moderate Re
For moderate flow speeds, the asymptotic characterization of the basin of attraction of the basic flow is no longer valid, particularly near the range of values of the double threshold that presumably lies at Re∼1800. Triggering transition in that region is possible, although the critical amplitude required may have a very complex dependency with Re. It is now generally accepted that for Re>1800 the Hagen–Poiseuille flow coexists with a complex turbulent state. As in previous sections, transition may be triggered by scaling up the amplitude of the initial perturbation. However, accurate refinements over this initial amplitude may lead to long transient orbits that neither trigger turbulence nor relaminarize. This refinement procedure allows the orbit to transiently approach and land on the stable manifold of a complex invariant set or edge state constituted by a surging amount of bifurcating complex solutions as Re is increased, leading to metastable chaotic dynamics (Nusse & Yorke 1989; Schneider & Eckhardt 2006; Skufca et al. 2006; Eckhardt et al. 2007). These so-called edge states are unstable and reside in the non-empty boundary separating the laminar and turbulent basins of attraction, and thus constitute a sort of chaotic saddle embedded within the critical transition threshold.
The solenoidal spectral approximation of the perturbation field u appearing in (3.6) is(4.1)with . Henceforth, the forcing term appearing in (3.6) is , so that the fluid is now driven at a constant mass flow. As usual, the weak formulation of equation (3.6) leads to a dynamical system for the amplitudes alnm(t) of the form(4.2)where the operators and are the resulting spectral projection of the linear time and space operators appearing in (3.6), is the projection corresponding to the nonlinear advective term and is the projection of the forcing (Meseguer & Mellibovsky 2007).
Unstable finite amplitude solutions of (4.2) cannot be approached by time integration and require suitable techniques to be computed. In this work, we search for exact solutions travelling in the axial direction with unknown streamwise speed c. In this case, the coefficients of the spectral expansion have the form(4.3)Formal substitution of the expansion (4.3) in (4.2), the solenoidal spectral scheme, leads to a time-independent system of nonlinear equations for each axial–azimuthal (l, n) mode,(4.4)with(4.5)The operators , and are the closure of , and over the (l, n) axial–azimuthal Fourier subspace, respectively. The system (4.4) is solved by means of a Newton–Krylov method (Meseguer et al. 2007).
(a) Approaching the stable manifold of the edge
The first stage of the approach to the edge of chaos is accomplished by perturbing the basic flow with disturbances of the form N1,2,3 described in previous sections. Starting from different initial perturbations, and through accurate refinements on their initial amplitudes, trajectories on the edge between turbulence and laminarity have been analysed at Re=2875 on a Λ=10 domain with . For instance, the time evolution of the driving axial pressure gradient (∇p)z, which enforces a constant mass-flux through the pipe, has been represented in figure 6a. It is clear from the plot that the three edge trajectories N1, N2 and N3 follow very similar erratic behaviours with occasional smooth approaching low-friction states followed by sudden outbursts towards more disordered states. In order to detect sporadic approaches of the edge trajectory to travelling wave solutions, an estimation of the phase speed corresponding to the l=1 mode of the state vector alnm over the edge is computed and the norm of in the nonlinear system (4.4) is evaluated as a function of time. This norm should exhibit a minimum whenever the edge trajectory is in the neighbourhood of a travelling wave solution. This phenomenon is clearly seen in figure 6b, where this norm has been plotted for the three runs. For instance, for t=1500, the N3 trajectory shows a clear minimum of ‖F‖.
(b) Underlying solutions
Some of the aforementioned residual drops correspond to approaches to a coherent state dominated by two high-speed streaks in the near-wall region, enclosing a low-speed streak located slightly off-centre, as shown in figure 7a. Together with the streaks, there coexists a strong pair of vortices thrusting low-speed fluid from the near-wall region towards the centre of the pipe, and at the same time feeding the low-speed streak with higher velocity fluid from the surrounding high-speed streaks. These low residual states were systematically used as initial conditions in the Newton–Krylov solver, consistently providing convergence to the same travelling wave solution shown in figure 7b, save for its axial–azimuthal degeneracy. This solution has a speed , posseses a shift and reflect (S&R) symmetry, and has recently been found by volume forcing homotopy (Pringle & Kerswell 2007). A Reynolds continuation of this travelling wave reveals that it originates at a symmetry-breaking pitchfork bifurcation from a reflection symmetric (RS) branch of waves. The continuation curve is shown in figure 8a, together with a snapshot of cross-sectional axial velocity contours of the symmetric travelling wave at Re=2875.
The Newton solver has been adapted for rotating travelling waves, as in Duguet et al. (2008). Feeding the same initial guesses from the critical trajectories to the Newton solver produced, in most of the cases, the same S&R travelling wave already described and shown in figure 7b. In at least one case (N2 critical trajectory, t=890), however, convergence was obtained onto a rotating travelling wave that bore very close resemblance to the non-rotating one in all respects (topology, energy contents and axial phase speed ), but for the S&R symmetry breaking, with a very slight rotation rate of rad. z-averaged axial velocity contours of the non-rotating and the rotating travelling waves are shown in figure 9. The non-rotating travelling wave appears symmetric upon z-averaging due to its S&R symmetry, while the rotating travelling wave is clearly non-symmetric. Their topological resemblance becomes evident when inspecting the velocity fields. Overall, the results presented in this section are in very good agreement with simultaneous work by Duguet et al. (2008), where travelling and also rotating waves were computed with similar methods.
In this work, we have provided a numerical exploration of the boundary of the basin of attraction of the Hagen–Poiseuille flow for high and moderate Reynolds numbers. For high Reynolds numbers, two different scenarios of transition have been addressed. For the autonomous scenario, the critical exponents were found to lie within the interval γ∈[−1.5,−1.0], strongly depending on the topological structure of the initial vortical perturbation. For the impulsive forcing, very good agreement was found with the experiments, providing the same exponent γ=−1. For moderate Reynolds numbers, the asymptotic scaling is not applicable and the exploration has been focused on approaching the critical threshold with a shooting method that allows us to visit the stable manifold of the edge state. Monitoring the edge trajectories, low drag transient coherent states have been identified and eventually converged to travelling and rotating–travelling waves making use of a Newton–Krylov method. Because these travelling waves are embedded within the chaotic saddle lying on the critical threshold, their relevance to transition seems beyond all doubts.
This work was supported by the Spanish Ministry of Science and Education (grant FIS2007-61585) and by the Catalonian Government (grant SGR-2004). We thank Prof. Bruno Eckhardt and Dr Tobias Schneider for enlightening discussions and for their kind hospitality during our short visit to Marburg. We also thank Dr Yohann Duguet for a very inspiring exchange of ideas during his visit to Barcelona.
One contribution of 10 to a Theme Issue ‘Turbulence transition in pipe flow: 125th anniversary of the publication of Reynolds' paper’.
- © 2008 The Royal Society