A universal transition to turbulence in channel flow

Abstract

Transition from laminar to turbulent flow drastically changes the mixing, transport, and drag properties of fluids, yet when and how turbulence emerges is elusive even for simple flow within pipes and rectangular channels^1,2. Unlike the onset of temporal disorder, which is identified as the universal route to chaos in confined flows^3,4, characterization of the onset of spatiotemporal disorder has been an outstanding challenge because turbulent domains irregularly decay or spread as they propagate downstream. Here, through extensive experimental investigation of channel flow, we identify a distinctive transition with critical behaviour. Turbulent domains continuously injected from an inlet ultimately decayed, or in contrast, spread depending on flow rates. Near a transition point, critical behaviour was observed. We investigate both spatial and temporal dynamics of turbulent clusters, measuring four critical exponents, a universal scaling function and a scaling relation, all in agreement with the (2 + 1)-dimensional directed percolation universality class.

Main

Transition to turbulence in open shear flows such as pipe flow and channel flow has been a difficult puzzle for over 130 years¹. In such flows, laminar flow becomes turbulent despite its linear stability^5,6,7. Also, turbulent structures tend to be localized; laminar states do not break up into turbulent states unless they are invaded by turbulent neighbours. If the tendency for invasion by a turbulent state increases, the turbulent state will eventually spread over the entire space. It is this behaviour that led Pomeau to conjecture that the spatiotemporal intermittency observed at the transition from laminar flow to turbulence belongs to the directed percolation (DP) universality class^8,9. DP is a stochastic spreading process of an active (turbulent) state with a single absorbing state¹⁰, to which diverse phenomena such as spreading of epidemics, fires, synchronization¹¹, and granular flows potentially belong¹⁰. Thus, if the transition is continuous and the interaction is short ranged, then universal critical exponents are expected^10,12. The linear stability of the laminar flow and recent experimental findings of two competing processes (namely decaying and splitting of a turbulent puff) in pipe flow¹³ qualitatively support this analogy including other shear flows such as plane Couette and Taylor–Couette flows^{14,15,16,17,18,19}. However, direct characterization of the transition has been lacking. This situation is presumably due to the extremely long timescale of pipe flow, thereby requiring experiments with extraordinarily long pipes to observe the critical phenomena. To overcome this difficulty, we chose a quasi-two-dimensional channel flow and forced the inlet boundary condition to be an active (turbulent) state. This enabled us to study the transition to turbulence as a surface critical phenomena. As a result, a clear transition between decay and penetration of the injected turbulent flow was observed. Quantification of the order parameter and the correlation length revealed critical behaviour of the transition in the experiment on shear flows; three independent critical exponents support the notion that the transition to turbulence in channel flow belongs to the DP universality class.

In channel flow, the Reynolds number (Re) is defined as Re = Uh/ν_K, where U is the centreline velocity of the parabolic profile, h is the half-height of the channel, and ν_K is the kinematic viscosity of the fluid (Fig. 1). Laminar channel flow (plane Poiseuille flow) is linearly stable up to a Reynolds number of Re_L = 5,772 (ref. 20). However a turbulent spot excited by a finite perturbation can grow and split to spread into extended spatial regions because of a global nonlinear instability even if Re is much smaller than Re_L(refs 21,22). To study this transition, an experimental set-up was configured. The flow channel has a length of 5,880 mm in the streamwise (x) direction, a cross-section of 5 mm in depth (the y direction), and a width of 900 mm in the spanwise (z) direction. Thus the aspect ratio of the channel is 2,352h × 2h × 360h. The flow dynamics in the (x, z) plane was visualized and recorded using a visualization technique and three charge-coupled device (CCD) cameras (see Methods). Instead of triggering turbulent spots by a local perturbation for each measurement, as in the previous experiments^13,21, turbulent flow is continuously excited in the buffering box through the use of a grid and injected from the inlet (x = 0), otherwise the flow remained laminar up to much higher Reynolds numbers (see Methods). Figure 1 shows the visualization of turbulent spots observed near the middle of the channel (x/h = 1,200) at Re = 810. Note that most of the turbulent flow injected at the inlet decayed quickly and became a laminar flow. Hence, any surviving turbulent flow tends to be visible as localized turbulent spots characterized by finer-scale disordered eddies surrounded by several streaks and clear laminar flows²¹. The typical size of the turbulent spot is about 40–80h.

Figure 1: Apparatus and snapshot of turbulent spots.

a, Schematic of the apparatus. The aspect ratio of the channel is 2,352h × 2h × 360h, where the depth 2h is 5 mm. b, Turbulent spots are visualized near the middle (x = 3 m) downstream location of the channel at Re = 810. The turbulent flows are injected by using a grid at the inlet (x = 0) of the channel. Visualization was assisted by means of micro-platelets and grazing angle illumination. Scale bar, 100 mm.

Figure 2 shows normalized intensity images of the flow pattern for three different Reynolds numbers. As shown in Fig. 2a for Re = 798, the injected turbulent structure separated into localized turbulent spots which quickly decayed as they propagated with the mean flow, and ultimately disappeared before reaching the channel exit. For Re ≥ 830, splitting and spreading of turbulent spots were clearly observed (see Fig. 2b for Re = 842). These processes contributed to the creation of turbulent clusters whose dynamics exhibited an intermittent stochastic nature in space and time. For sufficiently large Re values (for example, Re > 900), turbulent flow was sustained (see Fig. 2c for Re = 1,005).

Figure 2: Spatial variation of the flow across the transition.

a, Snapshots taken by three CCD cameras from left to right, respectively, at Re = 798. Quick decay of turbulent flow is evident. Colour represents the normalized image intensity. A black colour is assigned to the point where the image intensity is close to the laminar state (see the colour map). b, Snapshots at Re = 842. The intermittent nature of the turbulent spots can be seen. c, Snapshots at Re = 1,005. Saturation of the turbulent fraction is evident.

This set-up enabled a steady-state measurement of the area fraction of the turbulent region (the turbulent fraction ρ) for various values of x. The value of ρ, estimated by measuring the time fraction occupied by turbulent flow averaged over a protracted time period (approximately 40 min; that is, 100 times the length of the flow circulation time), was found to saturate for higher Re and for larger x, as shown in Fig. 3a. Therefore, the turbulent fraction was measured as a function of Re at several distant locations, x, satisfying x/h > 1,280 (see Fig. 3a). The area fraction of the active (turbulent) region is an order parameter in the DP transition which increases continuously from zero to positive values. Thus, the curves are fitted by the function

in the inset of Fig. 3a, where ɛ is the reduced Reynolds number. As a result, β = 0.58(3) and Re^c = 830(4) were obtained as the best fit values. The value of β was very close to the universal exponent of (2 + 1)-dimensional (that is, two-dimensional in space and one-dimensional in time) DP, β^DP = 0.583(3). Furthermore, the result Re^c = 830(4) is consistent with the results of direct numerical simulations for a channel flow, in which the global instability was reported as Re^c < 840 (ref. 23).

Figure 3: Critical behaviour of the turbulent fraction.

a, The turbulent fraction ρ versus Re is plotted at different downstream locations: x/h = 1,292 (orange square), x/h = 1,880 (blue diamond) and x/h = 2,096 (green square). Error bars represent standard deviation. Inset: a log–log plot of ρ as a function of reduced Reynolds number ɛ, where É› ≡ (Re − Re^c)/Re^c, with Re^c = 830(4). The solid blue lines are the best fit, ɛ^β with Î² = 0.58(3), for the data in 10^−3 < ɛ < 10^−1. Here, numbers in the parentheses denote 95% confidence intervals in the sense of the Student’s t distribution. The same applies to the following. Note that data points below Re^c are removed for fitting. A non-vanishing order parameter below Re^c due to a finite size effect exists as usual; however, relatively small systems can show remarkably clear power-law behaviour in numerical models exhibiting a DP transition (see Supplementary Information for simulation). b, The turbulent fraction as a function of distance x from the inlet where turbulence is created by a grid. Measurements were performed for six different x locations where the incidence angles and the reflected angles of the light were identical. The solid lines show the exponential fittings, ρ(x) ∼ exp(−x/L), applied for the data satisfying x/h > 1,040. Error bars represent standard deviation. Inset: log–log plot for L versus É›â€² ≡ (Re^c − Re)/Re^c. Error bars of the fitted values L are 95% confidence limits. The solid line is the best fit, L ∼ |ɛ′|^−ν with Î½ = 1.1(3).

Moreover, spatial variations of ρ(x) over space were investigated. The turbulent fraction ρ(x) showed clear exponential decays for Re values smaller than 803, whereas ρ showed saturations at constant values in space for Re = 904, as shown in Fig. 3b. Hence, transition between decay and penetration is evident. Thus, we fit ρ(x) with an exponential decay; ρ(x) ∼ exp(−x/L) for the data taken at Re < Re^c. The decay length L increased as Re^c was approached, with a power-law relationship, L ∼ |ɛ|^−ν. As a result, ν = 1.1(3) was obtained as the best fit (see the inset of Fig. 3b). This value is close to the critical exponent characterizing the divergence of temporal correlation length, ν_∥^DP = 1.295(6). For the temporal correlation length ξ_∥ and the spatial correlation length Î¾_⊥, the relations ${\mathrm{ξ}}_{âˆ\yen }∼{\mathrm{Î\mu }}^{−v_{âˆ\yen }^{\text{DP}}}$ and ${\mathrm{ξ}}_{âŠ\yen }∼{\mathrm{Î\mu }}^{−v_{âŠ\yen }^{\text{DP}}}$ hold, respectively, in DP. As ɛ approaches 0, the spatial correlation from the active wall becomes irrelevant compared with the temporal correlation as a result of the relation ν_∥^DP > ν_⊥^DP (and thereby ξ_∥ ≫ ξ_⊥, for ɛ ≪ 1), which holds in DP. Thus, the examination of spatial variation of ρ(x) is actually equivalent to the examination of a quenching dynamics of turbulence injected from the inlet that is conveyed downstream by the flow. This is the very reason ν_∥ was observed instead of ν_⊥ with respect to ρ(x). Therefore, the decay length L coincides with the survival length of the active cluster which defines the temporal correlation length ξ_∥ in DP (refs 24,25,26).

There are three independent static exponents that characterize the DP universality class: β, ν_∥ and ν_⊥. Numerical simulation on a simple directed bond percolation model with advection indicates that one can estimate the remaining exponent ν_⊥ by measuring distributions of the durations τ of the laminar state (laminar interval distribution) N(τ) at fixed downstream locations for Re > Re^c (see Supplementary Figs 6 and 7). Therefore, the distributions N(τ) at x = 3,200 mm were accumulated for 40 different z-positions within a half-span width (±225 mm) around the mid-height. For small τ values, a power-law distribution is expected near Re = Re^c reflecting the scale invariance of critical clusters¹⁰. Figure 4a shows the resulting N(τ) dependence for several different Reynolds numbers. We fit this by the power law N(τ) ∼ τ^−μ, with μ = 1.25(5), which is close to the universal exponent in DP, μ_⊥^DP = 1.204(2).

Figure 4: Critical behaviour of correlation length and universal scaling.

a, Laminar interval distribution N(τ) measured at a fixed downstream location, x/h = 1,280, for different Re. Power-law fit for Re = 846–832 approaches N(τ) ∼ τ^−μ, with μ = 1.25(5). The solid blue line corresponding to μ_⊥^DP = 1.204 is a guide for the eye. b, Complementary cumulative laminar interval distribution P(τ). The tails show exponential decays: P(τ) ∼ exp(−τ/ξ). c, Correlation length ξ versus reduced Reynolds number ɛ. The solid line is the fit, $\mathrm{ξ}∼{\mathrm{Î\mu }}^{-v_{âŠ\yen }}$, with ν_⊥ = 0.72(6). The dashed line corresponding to ν_⊥^DP = 0.733 is a guide for the eye. d, The data collapse for several different Re according to the scaling hypothesis. For Re^c, the value estimated in the experiment was used. For ν_⊥ and Î¼_⊥, the theoretical values for (2 + 1)-dimensional DP were used to assess whether the phenomena belong to the DP class or not.

To observe the tail of the distributions, a complementary cumulative probability, P(τ) ≡ ∫ _τ^∞N(t)dt/∫ ₀^∞N(t)dt was calculated, as shown in Fig. 4b. We defined the correlation length, ξ, by fitting the tail of P(τ) with an exponential function, P(τ) ∼ exp(−τ/ξ). As the transition point (Re^c) is approached, ξ increases substantially (Fig. 4c). Thus a best fit was determined in the form ξ ∼ ɛ^−ν, with an exponent ν = 0.72(6) for a small ɛ region (0.005 < ɛ < 0.06) in accordance with ν_⊥^DP = 0.733(3). Although the range of the power law is limited owing to the finite size of the system, the obtained exponents of μ_⊥, β and ν_⊥ consistently satisfy the universal scaling relation μ_⊥ = 2 − β/ν_⊥. As such, these results encourage the further exploration of universal features for the subject phenomena. Thus, a universal scaling hypothesis,

for P was introduced (see Supplementary Information for a detailed discussion and numerical validation of this hypothesis) with a universal scaling function g(x). By plotting the rescaled probability ${\mathrm{Î\mu }}^{-v_{âŠ\yen }({\mathrm{μ}}_{âŠ\yen }-1)}$ P(τ) as a function of the rescaled duration ${\mathrm{Î\mu }}^{v_{âŠ\yen }}\mathrm{τ}$, we find that several curves overlap (see Fig. 4d) when we choose Re^c = 830 in accordance with the previous result shown in Fig. 3a. All these results support that the transition can be understood as the DP process conveyed downstream by the flow.

In conclusion, the present result strongly supports the notion that the transitions to turbulence in shear flows belong to the (2 + 1)D DP universality class (the critical exponents obtained are summarized in Table 1). In fact, a similar conclusion has been reported recently by Shi et al. for shallow height Taylor–Couette flow, where (1 + 1)D DP universality was observed²⁷. Unveiling the ‘dynamical origin’^28,29,30 of the critical behaviour and quantifying dynamics by a spreading experiment are future challenges towards a deeper insight into the onset of turbulence.

Table 1 Summary of the critical exponents measured in this experiment.

Methods

Construction of the flow channel.

The channel walls were made of 25-mm-thick polymethyl methacrylate (PMMA) plates of optical surface quality. The entire 6 m (5,880 mm) channel comprised three pieces with 1,960 mm × 1,000 mm slots (see Supplementary Fig. 1). Both ends of each slot were reinforced by welding 50-mm-thick flanges to ensure the precision of the joint between two slots using an O-ring. The side walls were made of PMMA strips with dimensions 50 mm × 5 mm × 1,000 mm. When constructed in this way, the precision of the depth was ±0.1 mm. To avoid further deflection due to static pressure load in the channel, cross-braces were placed at 425-mm intervals along the channel. The working fluid is water. The channel inlet was connected to a buffering box by means of a smoothly curved contracting joint whose area contraction ratio was 1:20. To set a turbulent boundary condition, we placed a grid near the inlet. (When the grid is covered with seven layers of mesh screens, the flow remained laminar in a whole channel at least up to Re = 1,400. As the covering by mesh screen was not sufficient at the edge, turbulent flows did not decay near either end of the buffering box near z = 0 mm and z = 900 mm at Re = 1,400. Those turbulent flows injected from the inlet gradually grew and spread. Even in that case, there was no spontaneous nucleation of turbulent spots from the laminar state in the middle of the channel.) Velocity control was attained by electronically controlling the speed of the pump and the opening of the valve. We monitored the pressure gradient across the channel. The pressure gradient was almost constant during each measurement. The flow rate was measured by a flow meter (FD-UH40G, Keyence). The temperature of the water was controlled at 25 °C within an accuracy of ±0.1 °C.

Visualization.

As the measurement of the spatiotemporal dynamics of turbulent spots in a large space is problematic, we used a simple visualization using tracer particles. Metal-coated mica platelets (10–20 μm in diameter and 3 μm in thickness, Iriodin, Merck) were added to water for visualization. The concentration of the tracer was reduced to 0.04% in weight to keep the change of viscosity negligible (<0.1% according to Einstein’s law²³). Thin platelets tend to align perpendicular to shear stress, which implies parallel to the x–z surface in laminar flow states, whereas they rotate in turbulent spots. The grazing angle illumination gave moderate light reflections from the laminar regions towards the front, whereas scattering from the turbulent spots is omni-directional and its intensity deviated significantly from that of the laminar regions (see Supplementary Fig. 3). Six projectors (PJ4114NW, 3000 lumen, ultra-short focal length, Ricoh) were attached 250 mm apart and 300 mm above the (x, z) surfaces to illuminate the channel surface at a grazing angle to attain a reasonably uniform intensity of illumination. Three CCD cameras (1,608 pixels × 1,208 pixels, 10 frames s^−1) facing the centre of the x–z plane of each slot synchronously captured movies of the spatiotemporal dynamics of the flow of each slot. For the evaluation of the turbulent fraction, p(x) was measured at six positions (x = 0.65 m. 1.27 m, 2.68 m, 3.23 m, 4.70 m, 5.24 m) where the incidence angles from each of the six light sources to each measuring position in the channel were almost equal, and simultaneously the reflection angles from the measuring position to each of the three CCD cameras were almost equivalent. This choice was made to avoid unwanted inhomogeneity in the turbulent fraction ρ(x) due to the anisotropic nature of the light scattered from the platelets.

Image analysis.

Original images from movies captured by the three CCD cameras were normalized by the background images. Background images were obtained for the laminar state before and after each daily measurement. Next, histograms of normalized images were calculated at each pixel by accumulating time series from the movies for each Reynolds number. If the fluctuations of the normalized intensities at each pixel exceeded three times the standard deviation of the fluctuations of the laminar state, we determined that the pixel point lies in a turbulent cluster. For the resulting cluster structures in space, we used a criterion that the minimum size of a cluster is larger than h². A more detailed explanation of the image analysis is given in the Supplementary Information.