Purely elastic turbulence in pressure-driven channel flows

Significance Flows of complex fluids (such as polymers, colloids, emulsions, pastes, etc.) are abundant in everyday life. One can think of pouring shampoo or squeezing toothpaste from a tube, but also of fiber-spinning and plastic extrusion. Such flows often exhibit unexpected complexity with an archetypal example being elastic turbulence—chaotic motion exhibited by dilute polymer solutions. While its appearance is understood theoretically in flows with curved streamlines, little is known about elastic turbulence in straight pipes and channels. Here, we report simulations of elastic turbulence in a parallel flow geometry. Our results show that elastic turbulence is very different from its Newtonian counterpart and pave the way to understanding its mechanism.

M ost of the materials around us do not flow like New- tonian fluids do.Polymer rod-climbing (1), colloidal shear-thickening (2), and shear-banding in worm-like micelles (3) are just a few examples of such non-Newtonian behaviour that stems from the microstructure of these materials changing under external deformations.Despite their relevance to a wide spectrum of phenomena, ranging from microorganism propulsion (4) and mechanics of biological tissues (5) to polymer processing (6), time-dependent flows of non-Newtonian fluids are often beyond the reach of the current theory.An archetypal example of such unsteady flows is observed in dilute solutions of long, flexible polymer molecules (7).When subjected to a steady external forcing, polymer solutions exhibit instabilities to unsteady, chaotic flows often referred to as elastic turbulence (8).Unlike its Newtonian counterpart, which is driven by fluid inertia, elastic turbulence originates in strongly anisotropic elastic stresses developing under external forcing, and can occur even in the absence of inertia (9,10).The transition to elastic turbulence is then controlled by a dimensionless parameter, the Weissenberg number Wi -a product of the applied velocity gradient scale and the polymer relaxation time, that determines the strength of the elastic stresses and plays the same role that the Reynolds number, Re, does in Newtonian turbulence.
The route to elastic turbulence is best documented in model geometries with curved streamlines, where anisotropic elastic 'hoop' stresses lead to linear instabilities (11) followed by a transition to chaotic flows at higher Wi (12,13).Much less is known about flows of dilute polymer solutions in straight channels and pipes.Parallel shear flows of model polymer fluids are linearly stable under the experimentally relevant conditions (14), however, recent experiments present strong evidence of the presence of elastic turbulence in parallel shear geometries (15)(16)(17).Indeed, it was proposed that the transition to elastic turbulence in such geometries is a 'bifurcation from infinity' (18), i.e. it requires a finite-amplitude flow perturbation (19).
Recently, the first steps were made in confirming this instability scenario theoretically.The development was spurred by the discovery of a linear instability in parallel shear geometries (20) that only exists in a remote part of the parameter space, with Wi ∼ O( 103 ) and at extremely low polymer concentrations.Although this linear instability might not be accessible experimentally and is not directly responsible for triggering elastic turbulence, it was shown to give rise to two-dimensional nonlinear travelling-wave solutions (21,22).When continued outside this part of the parameter space, the travelling-wave solutions persist even in the absence of the underlying linear instability and can be found for a broad range of experimentally relevant parameters (23).Below, we refer to these as 'narwhal' solutions, reflecting the characteristic spatial profile of the associated polymer stretch.Although they resemble some features observed in experiments (15)(16)(17), the 'narwhal' states cannot be identified with elastic turbulence as they are not chaotic.
Despite these advances, it is presently not known how elastic turbulence is sustained.The main obstacle to further progress is the lack of understanding of the three-dimensional spatial distribution of the polymeric stresses and the associated fluid velocity of elastic turbulence, and how they evolve in time.
In the absence of experimental measurements, it is natural

Significance Statement
Flows of complex fluids (such as polymers, colloids, emulsions, pastes, etc.) are abundant in everyday life.One can think of pouring shampoo or squeezing toothpaste from a tube, but also of fibre-spinning and plastic extrusion.Such flows often exhibit unexpected complexity with an archetypal example being elastic turbulence -chaotic motion exhibited by dilute polymer solutions.While its appearance is understood theoretically in flows with curved streamlines, little is known about elastic turbulence in straight pipes and channels.Here, we report the first simulations of elastic turbulence in a parallel flow geometry.Our results show that elastic turbulence is very different from its Newtonian counterpart and pave the way to understanding its mechanism.

D R A F T
to rely on direct numerical simulations to provide structural understanding of the flow, and this route has proved to be very fruitful in studying Newtonian turbulence (24).However, simulations of polymeric flows suffer from numerical instabilities and are notoriously difficult to perform (25).Despite a significant effort, there are no reports of successful simulations of three-dimensional elastic turbulence in confined geometries in the absence of a linear instability.
To fill this gap, here we report direct numerical simulations of elastic turbulence in pressure-driven flow of a model dilute polymer solution.We study a model polymer fluid in a threedimensional straight channel with x, y, and z being Cartesian coordinates along the streamwise, wall-normal, and spanwise directions, respectively.The fluid fills the gap between two parallel plates and is driven by a constant pressure gradient applied externally (Fig. 1a).We employ periodic boundary conditions in the x-and z-directions, with Lx and Lz denoting the simulation box lengths along the respective axes.A dilute polymer solution is modelled by the simplified Phan-Thien-Tanner (sPTT) constitutive relation (26) given by where c is the polymer conformation tensor, v is the fluid velocity, p is the pressure, and x denotes a unit vector in the streamwise direction.Equations are rendered dimensionless by using d, U, d/U, and ηpU/d, and (ηs + ηp)U/d as the units of length, velocity, time, stress, and pressure, respectively.Here, d is the half distance between the channel plates, ηs and ηp are the solvent and polymeric contributions to the viscosity, and U is the maximal value of the laminar fluid velocity of a Newtonian fluid with the viscosity ηs + ηp at the same value of the applied pressure gradient.The flow is characterised by the following dimensionless quantities: the Reynolds number, Re = ρ Ud/(ηs + ηp), the Weissenberg number, Wi = λ U/d, the viscosity ratio that acts as a proxy for the polymer concentration, β = ηs/(ηs +ηp), the parameter controlling shear-thinning in the sPTT model, ϵ, and the stress diffusivity, κ.Here, ρ is the density of the fluid, and λ is its For all Wi studied here, the laminar flow is linearly stable (14).To study its non-linear stability, we introduce a random perturbation of magnitude δ (Materials and Methods), and follow its time evolution.For sufficiently low Weissenberg numbers, Wi < 80, simulations return to the laminar state independent of the perturbation magnitude δ.For Wi ≥ 80, there exists a finite value of δ, which is required to destabilise the flow (Figs.1b and SI Appendix Fig. S3a).Following the non-linear evolution of the flow in this regime, we observe irregular temporal oscillations of the kinetic energy (Fig. 1c) commensurate with unsteady, chaotic dynamics observed in elastic turbulence (8).These oscillations persist for many polymer relaxation times, measured in terms of t/Wi in our units (Figs.1c), signifying the presence of non-linear, chaotic flow states for Wi ≥ 80.The transition to purely elastic turbulence thus follows a bifurcation-from-infinity scenario (18,19), with finite-amplitude jumps in the observables separating the two states.The amplitude of the kinetic energy jump, however, is small (Fig. 1d), in line with a weakly sub-critical bifurcation that has recently been reported in experiments (16,28).We note that it can easily be mistaken for a linear instability.The amplitude of the jump of the midplane polymer stretch, on the other hand, is significantly higher (Fig. 1e), consistent with the polymer stress being the key dynamical variable of purely elastic turbulence.
Notably, close to the onset of elastic turbulence, the long time evolution of the non-linear flow states is followed by a sudden return to the laminar state (see Supplementary Movie 1), and we observed such events for Wi = 80 and Wi = 90, but not for Wi = 85 (Fig. 1c).Sudden relaminarisation events are a key characteristics of the transition scenario in linearly stable Newtonian parallel shear flows (29), indicative of a fractal laminar-turbulent boundary (30), finite turbulent lifetimes (31), and localised flow structures (32).Our observations, supported by further evidence below and by recent experimental reports of sudden splitting of localised structures in elastic pressure-driven channel flows (33), suggest a similar transition scenario for linearly stable purely elastic flows.
Our simulations provide direct access to the spatial distribution of the flow velocity and the polymeric stress and allow us to understand the structural features of elastic turbulence in straight channels.For all Wi studied, we observe that the largest deviations of the polymer stress from its laminar profile were mainly localised in a thin sheet around the channel centreline (Fig. 2a-c), while the flow is almost laminar close to the walls.This is in a stark contrast with Newtonian turbulence that exhibits the strongest fluctuations close to confining walls (34).For low values of the Weissenberg number, close to the saddle-node bifurcation at Wi ≈ 80, the polymer extension presents a spatially localised profile, similar to turbulent puffs and spots observed in Newtonian pipe and channel flows (35), respectively.For larger values of Wi, these localised structures proliferate throughout the domain in what resembles a percolation transition (29), exhibiting chaotic splitting and merging (see Supplementary Movies 1-3).Such events are reflected in the strongly intermittent temporal behaviour of the kinetic energy (Fig. 1c) and the midplane polymer stretch (SI Appendix Fig. S3a).Throughout these dynamics, the stress profile in the streamwise-wall-normal plane (Fig. 2d-f) preserves the overall features of the two-dimensional 'narwhal' states.This suggests that although they are linearly unstable (36), the subcritical 'narwhal' states organise the three-dimensional chaotic dynamics close to the onset of elastic turbulence.
We observe that the most prominent velocity features of elastic turbulence are associated with the streamwise velocity component.Since the chaotic stress dynamics are largely confined to a region around the centreline, the deviations of the mean streamwise velocity U (y) from its laminar profile are only noticeable in the middle of the channel (Fig. 1f) and constitute a few percent of the laminar centreline value.These observations are consistent with previous experimental measurements ( 16) that identified the centreline velocity fluctuations as an observable sensitive to the onset of elastic turbulence in channel flows.The spatial distribution of the streamwise velocity on the midplane (Fig. 2g-i) exhibits lowand high-velocity streaks similar to the structures reported in recent experiments (17).In contrast to Newtonian turbulence, the presence of streaks is not associated with the streamwiseoriented vortices, characteristic of the near-wall cycle (37) (Fig. 2j-l).Instead, we observe that the midplane streaky features and the chevron-like stress profile are reminiscent of the flow patters in viscoelastic Kolmogorov flow (38), and we hypothesise that a shear-layer instability plays a role in sustaining three-dimensional elastic turbulence in channel flows.

D R A F T
The final velocity features of elastic turbulence in our simulations are the spanwise-oriented vortices placed symmetrically above and below the midplane (Fig. 2g-i).Such structures are remarkably similar to the vortices reported in elasto-inertial turbulence at high Re (39) and, thus, our simulations provide support to recent claims that elastic and elasto-inertial turbulence have the same physical origins (14,40) and that there is a continuous connection between these phenomena in the Re − Wi space (20).
To analyse the statistical features of the ET dynamics, we study one-dimensional spatial spectra of the streamwise velocity fluctuations (Fig. 3a) at the channel midplane.We observe that the kinetic energy spectra are broadly consistent with the E u ′ u ′ ∝ k −4  x behaviour, with the extent of the scaling region increasing with Wi (SI Appendix Fig. S1).This power law, considered to be a hallmark of polymer-induced chaotic flows (10,38,41,42), indicates that the velocity field of elastic turbulence observed in our simulations is spatially smooth, unlike its Newtonian counterpart (43).As can be further seen from Fig. 3a, the error bars, indicative of one standard deviation, become large at high wavenumbers kx.We observe that there are instances in time where velocity and polymer stress fluctuations are confined to large and intermediate scales, with the corresponding spectra dropping off steeply at sufficiently large kx (SI Appendix Fig. S2).At other instances in time, the small scales fluctuate intensely.This indicates an intermittent breakdown of small-scale dynamics.We hypothesise that the dynamics at small scales are important to sustaining elastic turbulence; this, in turn, would explain why our simulations require large spatial resolution to successfully capture polymer-induced chaotic flows.Temporal statistical features are captured by the spectra of the centerline velocity (Fig. 3c), that scale with the frequency f as SU max (f ) ∝ f −2 .Spectral exponents in line with this value have been measured experimentally around the transition to turbulence in viscoelastic channel flow (16,28,44).

D R A F T
At small Re, the main dynamical variable of elastic turbulence is the polymer conformation tensor, while the velocity field is linearly enslaved to its time-evolution by virtue of Stokes' equation.To study its statistical behaviour, we plot one-dimensional spectra of Tr c (Fig. 3b) and observe that x .A spectral exponent of −2 is readily obtained by dimensional analysis from the velocity spectrum reported above and Stokes' equation.However, a small but finite deviation of the spectral exponent of the stress fluctuation from predictions based on dimensional analysis indicates the presence of additional effects determining the statistics of the stress fluctuations.This motivates further analyses, for instance concerning statistical self-similarity as a function of Wi.
Our simulations also provide indirect evidence that the chaotic dynamics of ET close to its onset is organised by unstable coherent structures, such as relative periodic orbits.In Fig. 3d, we plot the phase-space projections of high-dimensional chaotic trajectories corresponding to the time-series of the kinetic energy shown in Fig. 1c onto a two-dimensional subspace spanned by the deviation of the kinetic energy from the laminar profile, 1 − E(t)/E lam , and the normalised mean of the trace of the polymer stress, A(t) (Materials and Methods).We observe that for each Wi the phase-space trajectories evolve around points and closed curves in the two-dimensional subspace for long times.Similar observations in Newtonian flows were shown to indicate the presence of coherent states in the vicinity of the trajectories (31,32,45,46).As noted above, the spatial stress configuration resembling the two-dimensional 'narwhal' states (23) is clearly visible in Fig. 2d-f key part of the ET coherent structures.Finally, we note that our simulations exhibit a high degree of spatio-temporal intermittency, with the number of spatially localised states fluctuating in time (Supplementary Movies 2 & 3).These dynamics are a hallmark of the dynamical process based on splitting and decaying of spatially localised structures.Future work is required to demonstrate whether this process belongs to the directed percolation universality class (47), as is the case in Newtonian parallel shear flows (29,48,49).

Materials and Methods
Definition of the observables.The velocity field and the conformation tensor are decomposed into the mean profile and fluctuations, v = ⟨v(y)⟩x,z,t + v ′ (x, y, z, t), [4] c = ⟨c(y)⟩x,z,t + c ′ (x, y, z, t), [5] where ⟨. . .⟩x,z,t denotes the spatial average along the streamwise and spanwise directions, and the temporal average in the statistically stationary state.In particular, we study the mean streamwise velocity profile, U (y), and the instantaneous deviation of the streamwise velocity from its mean profile, u ′ .The maximal value of the mean streamwise velocity is denoted by Umax and always corresponds to U (y = 0).The laminar values of the streamwise velocity and the conformation tensor are denoted by U lam (y) and c lam (y), respectively.The instantaneous ratio of the kinetic energy to its laminar value is defined through This quantity is used in Fig. 1b (without the normalisation), Fig. 1c, and Fig. 3d, while Fig. 1d shows its value averaged over time in the statistically stationary state.Similarly, the instantaneous ratio of the volume-averaged polymer stretch to its laminar value is defined through To aid the visualisation of the flow structures, in Fig. 2 we introduce the relative polymer stretch, defined as , [8] that allows for comparison between simulations at different Wi.One-dimensional spatial spectra of the streamwise velocity and the trace of the conformation tensor fluctuations at the channel centre-plane are defined as where f is the frequency, and F denotes the temporal Fourier transform.

Q-criterion.
The Q-criterion ( 50) is a Galilean-invariant vortex identification measure that is based on the second invariant of the velocity gradient tensor where are the strain rate and the vorticity tensor, respectively, and indices denote Cartesian components of the vectors.A vortex is defined as a region in the flow where Q > 0.
Direct numerical simulations.Equations Eq. ( 1)-Eq.( 3) are solved numerically with an in-house MPI-parallel fully-dealiased pseudospectral code developed within the Dedalus framework (51).The velocity, conformation tensor, and pressure are represented by a Fourier-Chebyshev-Fourier spectral decomposition along the streamwise, wall-normal, and spanwise directions, respectively, with the spectral resolution being given by the number of modes in each direction (Nx, Ny, Nz).Production runs are carried out using Nx = Nz = 256 and Ny = 1024, comprising approximately 7•10 8 degrees of freedom.Temporal discretisation employs the semi-implicit backward differentiation scheme of order four ( 52) with time step of 0.005.All simulations have been evolved for at least 130 polymer relaxation times (10 4 -10 5 time units depending on Wi) in a statistically stationary state.Data are sampled at intervals of 10 time units for full-state snapshots and 0.1 time units for the kinetic energy.A summary of all simulations is provided in SI Appendix Table S1.A typical production run was carried out on 16384 cores for approximately 260 hours.We employed three types of initial conditions in our simulations.In Fig. 1b, we perturbed the laminar state by adding δ|ξ| to the cxx component of the conformation tensor.Here, ξ is a positiondependent random Gaussian noise with zero mean and unit variance, and δ is the amplitude of the noise.The absolute value ensured that the perturbation preserved the positive-definiteness of the conformation tensor.This strategy resulted in very long runs to reach the statistically stationary state, we have only used it in Fig. 1b to demonstrate the existence of a finite-amplitude threshold, and in SI Appendix Fig. S3.For the second strategy, we embedded the two-dimensional 'narwhal' states (23) into our three-dimensional domain with a small amount of random noise added to break the translational symmetry along the spanwise direction.This strategy was employed for Wi = 100, and the quick destabilisation of the translationally invariant 'narwhal' state is visible at early times in Supplementary Movie 2. All other simulations presented here were started from the chaotic states obtained from the previous runs at a new value of Wi.
To demonstrate that different initial conditions produce the same statistical steady states, we continued the simulation at Wi = 150 with the highest noise amplitude (δ = 625 in Fig. 1b) until it reached a statistically steady-state.Additionally, to demonstrate convergence with respect to spatial discretisation, we repeated the same run at a reduced spatial resolution Nx = Nz = 128 and Ny = 512.In SI Appendix Fig. S3a, these runs are compared to the Wi = 150 simulation from Fig. 1c.Although the dynamics are strongly intermittent in all cases, visual inspection confirms convergence to a statistically steady state.This observation is further supported by SI Appendix Fig. S3b where we find no statistically significant differences between the midplane polymer stretch in these simulations.

Fig. 1 .
Fig.1.Characterising the transition to elastic turbulence.a, Flow geometry.b, Temporal evolution of the time derivative of the kinetic energy in simulations at Wi = 150 perturbed by random noise of amplitude δ (Materials and Methods).For small δ, the derivative reaches zero up to machine precision, indicating relaminarisation, while at larger values of δ it grows in time.The subsequent time evolution of the run with the largest noise amplitude (SI Appendix Fig.S3a) shows that it develops into a fully chaotic state.c, Temporal evolution of the kinetic energy normalised by its laminar value.Time is measured in terms of polymer relaxation time λ (t/Wi in our dimensionless units).d and e, Bifurcation-from-infinity: the finite-amplitude jumps in the kinetic energy (d) and the midplane polymer stretch (e) from their laminar values.f, Mean streamwise velocity profiles.(Inset) Deviation of the mean streamwise velocity from its laminar profile.

Fig. 2 .
Fig. 2. Instantaneous spatial structure of elastic turbulence at various Wi.The mean flow direction is indicated by the white arrows.a-c, Midplane profile of the relative polymer stretch, Ψ(c) (Materials and Methods).d-f, Vertical slices of Ψ(c) at several positions in the channel demonstrating variations in the y-direction.g-i, The midplane profile of u ′ /Umax and the isosurfaces of the Q-criterion (Materials and Methods) visualised for values −0.002 and 0.002 in blue and red, respectively.The simulation domain is repeated twice in the streamwise direction to improve readability of the figure.j-l, Streamwise slices of u ′ /Umax with the in-plane streamlines at x = 5.

Fig. 3 .
Fig. 3. Spatio-temporal properties of elastic turbulence at various Wi.One-dimensional spatial spectra of the streamwise velocity component (a) and the trace of the conformation tensor (b).c, Temporal spectrum of the centerline velocity.In (a-c), the dashed lines indicate tentative power-law scaling laws, the shaded regions indicate one standard deviation, and the data have been shifted vertically by the same offset to improve the readability of the figure.d, Phase-space projections of the chaotic trajectories.

]
Sc(kx) =dkz Tr ĉ′ (kx, y = 0, kz)2 t , [10]where ⟨. . .⟩t denotes the temporal average calculated in the statistically stationary state, and • denotes the two-dimensional Fourier transform along the x-and z-directions with kx and kz being theD R A F Tcorresponding components of the wavevector.The temporal spectra of the centerline velocity Umax are calculated as

Fig. S3 .
Fig.S3.Studying the influence of the initial condition and spatial resolution for Wi = 150.a, Time evolution of the midplane polymer stretch.Run 1 employed a configuration taken from the Wi = 100 simulation as its initial condition.The data is re-plotted from Fig.1cof the main text.Runs 2 and 3 were started from the laminar state perturbed by random noise with the amplitude δ = 625 (Materials and Methods).The small-t data of Run 2 are plotted in Fig.1b.Runs 1 and 2 employed the production spatial resolution, while Run 3 employed a reduced (halfed) resolution (Materials and Methods).b, Statistical analysis of the data from (a) evaluated in the statistical steady state.The circles denote the averages, while the lines indicate one standard deviation.In the absence of statistically-significant or systematic discrepancies between the runs, we conclude that all three simulations have converged to the same chaotic steady state.

Movie S1 .
Elastic turbulence at Wi = 80.The simulation rapidly converges to a localised turbulence structure that persists for a long time before suddenly relaminarising.Movie S2.Elastic turbulence at Wi = 100.The simulation is started from a two-dimensional 'narwhal' state, translationally-invariant along the spanwise direction, perturbed by a small amount of noise.Early time evolution shows that this state is unstable and the simulation quickly reaches a chaotic steady state.Movie S3.Elastic turbulence at Wi = 150.A strongly intermittent simulation exhibiting splitting and merging of localised coherent structures.

Table S1 . Parameters and key observables used in the statistical analyses. E and ⟨Tr c⟩ are the volume-and time-averaged kinetic energy and the trace of the conformation tensor
, respectively, while E lam and Tr c lam are the corresponding laminar values.tstart (