A unified view of elastic and elasto-inertial turbulence in channel flows

Abstract

The chaotic flow of elastic fluids at low Reynolds number ($Re$) is typically distinguished into elasto-inertial and elastic turbulence (EIT/ET). However, the clear separation among these two turbulent regimes in parallel flows with a gradual $Re$ decrease remains elusive. Here, spanning $Re$ over four orders of magnitude, we investigate the statistical and structural nature of the flows computed by fully resolved 3D numerical simulations, revealing that as inertia becomes sub-dominant all flows share similar dynamical properties. Supported by the results, we thus propose a unified view of fully developed ET and EIT in parallel flows.

Turbulence is an open problem in fluid transport applications, as it significantly enhances the drag and energy expenditures. Nevertheless, it is beneficial in mixing problems since it introduces an additional diffusion mechanism on top of the classical viscous one. Indeed, viscous diffusion is often too slow for practical applications, but it is the only mechanism operating in low Reynolds number flows, which are typically laminar. The Reynolds number, $Re$, is the non dimensional parameter conventionally adopted to describe fluid flow, representing the relative importance of inertia compared to viscosity. Micro-scale mixing in low Reynolds number flows is a widely investigated research topic Stroock et al. (2002); Ober et al. (2015); Han et al. (2022) for its high biological relevance, e.g., in drug delivery applications Vargason et al. (2021) and in the analysis of specimens for disease diagnosis Yager et al. (2006). Among the various mixing techniques, one is to exploit elastic turbulence Groisman and Steinberg (2000), a peculiar chaotic state unique of viscoelastic fluids which are characterised by the coexistence of viscous and elastic effects.

The description of viscoelastic fluid flows necessitates of an additional parameter accounting for the ratio of the elastic and flow time scales: the Deborah number, $De$. Comparing the Deborah to the Reynolds number allows to identify distinct regimes of motion: when inertial effects are dominant turbulent drag reduction Serafini et al. (2022) (TDR) occurs, while elastic turbulence Groisman and Steinberg (2000) (ET) is observed for dominant elasticity. In between these two extrema, an intermediate state dubbed elasto-inertial turbulence Dubief et al. (2023) (EIT) has been usually identified, where both elastic and inertial effects are relevant. There is no universal consensus on the definitions of EIT and ET; yet, for the purposes of our investigation, we follow Groisman and Steinberg (2000) and denote with ET the chaotic flow observed in polymer solutions at small values of $Re$ and large values of $De$ ($Re\lesssim De$). Furthermore, according to Samanta et al. (2013), we dub EIT the chaotic motion that sets in at lower Reynolds number than Newtonian turbulence due to the action of fluid elasticity, in a parameter range where inertia cannot be neglected: when $Re$ and $De$ are comparable ($Re\approx De$).

Figure 1: Side views of the channel, with the mean flow going from left to right. We show the spanwise velocity fluctuations normalized with the mean turbulent kinetic energy, $w^{\prime}/\sqrt{K}$, ranging from -0.3 (blue) to 0.3 (green); superimposed are the local isosurfaces of the wall-normal fluctuations, $v^{\prime}/\sqrt{K}$, at $-0.05$ (orange) and $+0.05$ (violet).

EIT occurs at significantly lower Reynolds numbers than TDR, where the laminar flow of a Newtonian fluid would be linearly stable, and it is characterised by different dynamical properties Samanta et al. (2013). When inertial effects vanish, the flow can still become turbulent (ET) due to the onset of different elastic instabilities. It is in fact well established that different instability mechanisms are responsible for the transition to a chaotic state at moderate Garg et al. (2018) and vanishing Berti and Boffetta (2010) values of $Re$, thus motivating the distinction between EIT and ET in the stability framework Dubief et al. (2023). In curvilinear wall bounded flows, ET is triggered and sustained by a force directed as the flow curvature, called hoop stress Bird (1987); Groisman and Steinberg (2001). Recent experiments Jha and Steinberg (2021); Li and Steinberg (2023) and numerical simulations Song et al. (2022); Lellep et al. (2024) have highlighted that ET can also be attained in parallel shear flows, but a comprehensive characterisation of its mechanisms is still elusive. After the discovery of a linear elastic instability in the parallel Kolmogorov flow Boffetta et al. (2005), the weakly non-linear regime following the instability has been studied Bistagnino et al. (2007) and observed to culminate in ET Berti et al. (2008). Arrowhead-shaped traveling elastic waves were first observed in this scenario Berti and Boffetta (2010) and later found also in the bulk of channel flows in the EIT Sid et al. (2018); Garg et al. (2018); Dubief et al. (2022); Buza et al. (2022a) and ET regimes Morozov (2022), suggesting a possible continuous path between the two Dubief et al. (2023) and hinting towards the existence of a unique sustainment mechanism common to EIT and ET Khalid et al. (2021a). Channel flows at moderate to vanishing inertia have been recently observed to share also a different instability Beneitez et al. (2023), localised close to the walls. Nevertheless, due to a lack of measurements, it is still unclear if the fully developed three-dimensional EIT and ET chaotic states are the same, connected or completely decoupled.

In this Letter we show that fully developed EIT and ET in planar channel flows share a similar structure, notwithstanding the instabilities leading to them. Addressing the fundamental question posed by Dubief et al. (2023), we thus clarify the relation between EIT and ET within this setup. Our fully-resolved three-dimensional numerical simulations, spanning Reynolds numbers over four orders of magnitude (as in FIG. 1), allow us to bridge the different regimes (TDR, EIT and ET) and clarify their relation. For every Reynolds number considered in this study we observe the system to converge towards a self-sustained chaotic state even when the flow of a Newtonian fluid would instead be laminar, and we identify only two different regimes, one at large Reynolds number (TDR) and one at smaller ones (ET).

To show this, we simulate the flow of a viscoelastic incompressible fluid between two indefinite planar walls at a distance of $2h$ from each other, described in a Cartesian coordinate frame with the $x$, $y$ and $z$ axes aligned with the streamwise, wall-normal and spanwise directions, respectively. The system is governed by the conservation of mass and momentum,

$$\mathbf{\nabla}\cdot\mathbf{u}=0,$$ (1)

$$\rho\left(\partial_{t}\mathbf{u}+\mathbf{\nabla}\cdot\mathbf{u}\mathbf{u}\right)=-\mathbf{\nabla}p+\mu_{s}\mathbf{\nabla}\cdot(\mathbf{\nabla u}+(\mathbf{\nabla u})^{T})+\rho\mathbf{f}+\mathbf{\nabla}\cdot\mathbf{T},$$ (2)

where $\mathbf{u}$ is the velocity, $\rho$ the density, $p$ the pressure and $\mu_{s}$ is the dynamic viscosity of the Newtonian solvent in which the polymer molecules are dissolved. $\mathbf{f}$ represents the forcing term needed to drive the flow at a constant flow rate, while $\mathbf{T}$ is the tensor accounting for the non-Newtonian behaviour. In our study, $\mathbf{T}$ derives from a finitely extensible nonlinear elastic model with Peterlin’s closure (FENE-P Peterlin (1966), equivalent to a bead-spring-type model Vincenzi et al. (2007)) and relates to the polymer conformation tensor, $\mathbf{C}$, as $\mathbf{T}(\mathbf{C})=(\mu_{p}/\lambda)(f(tr(\mathbf{C}))\mathbf{C}-(L_{max}^{2}\mathbf{I})/(L_{max}^{2}-3))$, where $L_{max}$ denotes the maximum extensibility of the polymer chains, $\mathbf{I}$ is the tensorial identity, $\mu_{p}$ is the dynamic viscosity of the polymer, and $\lambda$ is its relaxation time. In particular, $f(s)=L_{max}^{2}/(L_{max}^{2}-s)$. Finally, $\mathbf{C}$ evolves in time according to its own transport equation

$$\partial_{t}\mathbf{C}+(\mathbf{u}\cdot\mathbf{\nabla})\mathbf{C}-\mathbf{C}\cdot\mathbf{\nabla u}-(\mathbf{\nabla u})^{T}\cdot\mathbf{C}=-\displaystyle\frac{1}{\lambda}(f(tr(\mathbf{C}))\mathbf{C}-\mathbf{I}).$$ (3)

The regime of motion of the viscoelastic fluid is identified by the specification of the Reynolds and the Deborah numbers, along with the solvent-to-total viscosity ratio, $\beta$, and $L_{max}$. In this work we adopt the definition of “bulk” Reynolds number, $Re=hU/\nu$, based on the mean velocity within the channel, $U$, and the total kinematic viscosity of the polymeric solution, $\nu=(\mu_{s}+\mu_{p})/\rho$; we therefore decrease $Re$ in $\{2800,1500,500,50,5,0.5\}$. Furthermore, we complement our results with those available in the literature by Lellep et al. (2024) at even smaller Reynolds number, $Re=0.01$. The Deborah number is defined as $De=\lambda U/h$ and imposed equal to $50$ throughout the whole study, while $L_{max}$ and $\beta=\mu_{s}/(\mu_{s}+\mu_{p})$ are set to $500$ and $0.9$, respectively. We start each simulation from a fully developed field of the closest case at a higher Reynolds number, and advance it in time according to the methods detailed in the Supplemental Material . Our parameters are chosen consistently with previous stability analysis Buza et al. (2022a) and yield values of the elasticity number $El=De/Re$ in $\{0.03,0.1,1,10,100\}$, denoting how the investigated range equally encompasses cases dominated by inertia ($El\ll 1$) and elasticity ($El\gg 1$), as well as the intermediate range. The results presented in the following have been replicated also with the Oldroyd-B liquid, thus proving their robustness against the polymer model adopted; furthermore, we have tested the effect of different viscosity ratios $\beta$, as detailed in the Supplemental Material (Note, 2).

Figure 2: Flow statistics. In $(\mathbf{a})$ we report the viscous (dashed), turbulent (dashed-dotted) and polymer (dotted) shears, summing up to the total shear stress $\tau$ (solid line). At $Re\leq 1500$ viscosity and elasticity dominate throughout the whole half-channel height $h$. All the contributions integrated along the wall-normal direction, in ($\mathbf{b}$), are thus dominated by the viscous and polymer shears, but for the case at $Re=2800$, dominated by the turbulent shear. Consistently, the polymer stretching is enhanced, as shown in $(\mathbf{c})$ by the PDFs of the trace of $\mathbf{C}$ at the channel centerline: their barycenters are far from $tr(\mathbf{C})=3$. $(\mathbf{d})$ shows the J-PDFs of the normalized streamwise and wall-normal velocity fluctuations at $y=h/2$, at $\{0.1,0.3,0.5,0.7\}$ of their maximum. The typical turbulent shape is recovered at $Re=2800$, while in all the other cases the fluctuations decorrelate. The TKE spectra $(\mathbf{e})$ show a good collapse for all the low-Reynolds number cases. We show additional data with Oldroyd-B at $Re=5$ (black circles) and from Lellep et al. (2024) at $Re=0.01$ and $De=80$/$150$ as dashed/dashed-dotted lines, respectively. The $-4$ slope contrasts sharply with the typical $-5/3$ decay, recovered in TDR, pointing towards a unified view of EIT and ET.

At all the Reynolds numbers but the highest one, qualitatively similar coherent fluctuations of the velocity are superimposed to a mean flow which only slightly departs from the laminar one (more details on this in the Supplemental Material Note, 2). As visible in FIG. 1, large structures span the center of the domain, while smaller ones characterize the region near the wall. All of them exhibit a quasi-periodic spatial arrangement superimposed to random perturbations. A more quantitative characterization of the flow fields is attained looking at the balance of the total shear stress, $\tau$, along the half-channel height, $h$. The conservation of linear momentum requires the wall-normal derivative of $\tau$ to balance the driving pressure gradient, $d_{y}\tau=d_{x}P$, imposing a linear trend of $\tau$ with $y$. Three different contributions sum together in the total shear stress: a polymeric one due to the shearing component of the viscoelastic stresses $\mathbf{T}$ ($\langle T_{xy}\rangle$, where angle brackets denote averaging in time and along the homogeneous directions), a purely viscous one coming from the mean velocity profile $\langle u\rangle$ ($\rho\nu\partial_{y}\langle u\rangle$), and a turbulent one caused by the correlated fluctuations of the streamwise and wall-normal velocities ($-\rho\langle u^{\prime}v^{\prime}\rangle$, where the primes indicate fluctuations). These contributions are reported in the plots of FIG. 2$a$.

As expected for the TDR regime at $Re=2800$, the viscoelastic contribution remains small almost everywhere and viscosity is mostly relevant close to the wall, while the turbulent contribution dominates the center of the channel. For $Re\leq 1500$, instead, the shear stress balance is completely different, as the polymer and viscous shears become relevant throughout the channel and the turbulent one is depleted. To better weight the relative importance of the different contributions of the shear stress balance, we integrate the momentum balance along the wall-normal direction ($\int{(\langle T_{xy}\rangle+\rho\nu\partial_{y}\langle u\rangle-\rho\langle u^{\prime}v^{\prime}\rangle)dy}=h^{2}\enspace d_{x}P$), as represented graphically in FIG. 2$b$. The turbulent shear is the most relevant component at $Re=2800$, while it is practically zero in all the other cases; those are dominated by the viscous shear, along with a non-negligible polymeric contribution (individually shown in the Supplemental Material Note, 2). Significantly, the relative amplitude of the different contributions remains unchanged for all the Reynolds numbers (ET) but the highest one (TDR).

The previous results warrant two important clarifications. First, the relative small contribution of the polymers does not mean that the polymers are not stretched. Indeed, the polymer molecules are significantly stretched regardless of the case considered, as shown by the probability density functions (PDF) of the trace of $\mathbf{C}$ reported in FIG. 2$c$. Also here, a clear separation can be appreciated between the case at $Re=2800$ and all the others. Second, the vanishing magnitude of the turbulent shear when $Re\leq 1500$ does not mean that the turbulent fluctuations disappear, but rather suggests that the streamwise and wall-normal velocity fluctuations decorrelate. We therefore compute the joint probability density function (J-PDF) of the streamwise and wall-normal velocity fluctuations at a wall distance of $h/2$, reported in FIG. 2$d$, where velocities are made dimensionless with the mean turbulent kinetic energy in the channel, $K$, for a better comparison of the different Reynolds numbers. In the high Reynolds case the isolines exhibit an oval shape with the major axis inclined downward, denoting a negative correlation among the two velocity components, and ejections ($u^{\prime}<0$, $v^{\prime}>0$) dominate over sweeps ($u^{\prime}>0$, $v^{\prime}<0$), analogously to the Newtonian case Tardu (2014). At this wall distance, the J-PDF peaks in the fourth quadrant. This picture completely changes when decreasing the Reynolds number, as the streamwise velocity decorrelates from the wall-normal one and its fluctuations heavily dominate in intensity; the J-PDF consequently appears flattened along the x-axis, explaining the apparent disappearance of the turbulent shear stress from FIG. 2.

The features of the low Reynolds flow fields described up to this point would appear compatible with an incomplete relaminarisation of the flow; nevertheless, its turbulent nature is confirmed by the full multi-scale Fourier’s spectra $E$ of the turbulent kinetic energy (TKE) at the centerline of the channel. Those are characterised by a power-law decay of $E(k_{x})\sim k_{x}^{-4}$ over roughly a decade of streamwise wave numbers $k_{x}$ (FIG. 2$e$), in good agreement with recent observations in a similar setup at lower Reynolds number Lellep et al. (2024) and in homogeneous isotropic elastic turbulence Singh et al. (2024). Contrarily, in the TDR regime at $Re=2800$, our data exhibit a tendency towards the conventional $E(k_{x})\sim k_{x}^{-5/3}$ trend expected for Newtonian turbulence at larger values of the Reynolds number. Once more, no significant difference is observed among all the low Reynolds number cases; the robustness of this statement to significant variations of the polymer concentration is detailed in the Supplemental Material Note (2). Our results therefore point towards a unified view of EIT and ET in three-dimensional fully developed planar channel flows.

Figure 3: Modal decomposition, showing the temporal frequencies $(\mathbf{a})$, $\omega$, of the two most dynamically relevant modes $(b,c)$ for five different values of the Reynolds number. The error bars represent the uncertainty of the decomposition algorithm calibration process Note (2). The colored panels $\mathbf{b,c}$ report the typical shapes of the streamwise, $u$, and spanwise, $w$, velocity components of the modes (sampled from the case at $Re=50$), where positive isosurfaces are colored in purple and negative isosurfaces in orange. The identification of two dominant modes common to all the low Reynolds number cases (consistent with previous stability analysis results) corroborates the unified view of EIT and ET.

To further confirm this scenario, we apply a purely data-driven algorithm to identify and extract the coherent structures in the flow at different Reynolds numbers. In particular, we adopt the higher order dynamic mode decomposition Vega and Le Clainche (2020), a dimensionality reduction algorithm which provides a finite set of modes, each oscillating in time with a specific frequency $\omega$. Past stability analysis of viscoelastic two-dimensional channel flows suggests that travelling waves originated at a relatively high-$Re$ bifurcation Garg et al. (2018) persist in EIT Page et al. (2020); Khalid et al. (2021b); Dubief et al. (2022); Buza et al. (2022a) and ET Buza et al. (2022b); Morozov (2022). These structures are mainly located away from the walls, in the central region of the channel, and are sometimes termed “arrowheads” Page et al. (2020). In ET, these coexist with a recently discovered wall-localised linear instability caused by the presence of small but non-zero diffusivity of the polymer stresses Beneitez et al. (2023); Couchman et al. (2024). The diffusive instability perturbs the flow close to the walls, in the region where the travelling waves convert most of the kinetic energy into elastic energy Beneitez et al. (2024). Our analysis suggests that all the turbulent flows at $Re\leq 1500$ can be described by the superposition of two main sets of modes for all the Reynolds numbers but the highest one, once more corroborating the unified view of EIT and ET suggested by the analysis above. The frequencies of these modes are shown in FIG. 3$a$ as a function of the Reynolds number, and their typical shapes are reported in FIG. 3$b$ and $c$. We observe that the fast modes are mainly located in the center of the channel, while the slow modes are localised close to the walls; their temporal frequencies do not undergo a significant variation with the Reynolds number. The shapes and locations of the modes are reminiscent of the results from simulations Lellep et al. (2024) and stability analysis Buza et al. (2022b); Beneitez et al. (2023).

In this work we have simulated the turbulent flow of a viscoelastic fluid in a planar channel setup across a wide range of Reynolds numbers from TDR to ET, reporting instabilities and flow structures consistent with the available literature Page et al. (2020); Beneitez et al. (2023); Dubief et al. (2023). Our results identify two distinct turbulent regimes. First, one at large Reynolds number dominated by the turbulent shear stresses with the classical energy spectra decay predicted by Kolmogorov’s theory Kolmogorov (1941) $E\sim k_{x}^{-5/3}$, which corresponds to the extensively investigated TDR regime Lumley (1969). Next, one at lower Reynolds number with highly elongated polymer chains producing strong polymer stresses with the energy spectra showing a fast fall-off $E\sim k_{x}^{-4}$, consistently with recent results in ET even at Reynolds numbers lower than those here Singh et al. (2024); Lellep et al. (2024). We do not find any other fully-developed state in between these two, and conclude that the peculiar chaotic state attained in viscoelastic planar channel flows remains unchanged by a variation in the Reynolds number once inertial effects are sub-dominant. Thus, within this set-up, fully-developed EIT Dubief et al. (2023) and ET Steinberg (2021) share a similar structure.