1. Introduction
The reduction of turbulent drag using polymer additives has been an active area of research due to its application to save energy during liquid transport (Xi Reference Xi2019) ever since it was discovered that the addition of a tiny amount (
$\sim O(10)$
ppm) of high molecular weight polymer suppresses inertial turbulence and dramatically reduces frictional drag (Toms Reference Toms1949). Although polymer additives can significantly reduce frictional drag, the flow does not become laminar; hence the drag is always higher than the laminar value. This phenomenon is attributed to the presence of a self-sustaining weakly fluctuating flow state known as the maximum drag reduction (MDR) asymptote (Virk et al. Reference Virk, Mickley and Smith1970). At MDR, the frictional drag is nearly independent of the type and concentration of polymers; understanding MDR is an important area of research because it sets a limit for the attainable amount of drag reduction. Earlier studies of MDR mainly focused on the ensemble statistics of the velocity fluctuations; however, recent advancements in the field attempt to understand the dynamic evolution of the MDR state (Xi Reference Xi2019). The discovery of elastoinertial turbulence (EIT), a self-sustaining chaotic state resulting from the interplay between inertia and elasticity, provides some understanding of the dynamics of the MDR state (Samanta et al. Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013; Datta et al. Reference Datta2022; Dubief et al. Reference Dubief, Terrapon and Hof2023). In this paper, we investigate the dynamics of EIT in pipe flow, and discover dominant structures underlying its dynamics.
The EIT can sustain in the parameter regime characterised by the Reynolds number
$\textit{Re}$
(a dimensionless number representing the ratio of inertial to viscous forces) and the Weissenberg number
$\textit{Wi}$
(the product of polymer relaxation time and characteristic strain rate), where the flow is linearly stable, suggesting a nonlinear route for transition to EIT. Viscoelastic two-dimensional (2-D) channel flow becomes linearly unstable due to a wall mode instability at sufficiently large
$\textit{Re}$
(Zhang et al. Reference Zhang, Lashgari, Zaki and Brandt2013) and a centre mode instability at large
$\textit{Wi}$
(Page et al. Reference Page, Dubief and Kerswell2020; Khalid et al. Reference Khalid, Chaudhary, Garg, Shankar and Subramanian2021), with parameter values significantly higher than the EIT parameter regime of interest here. The wall mode instability leads to the emergence of a travelling wave known as the Tollmien–Schlichting (TS) wave (Drazin & Reid Reference Drazin and Reid1981), and the centre mode instability leads to the emergence of another distinct travelling wave having an ‘arrowhead’ structure (Page et al. Reference Page, Dubief and Kerswell2020). The origins of these travelling waves are highly subcritical, and they can persist subcritically to the parameter regimes relevant to EIT. Inspired by these facts, it was speculated that the nonlinear excitation of either the wall mode (Shekar et al. Reference Shekar, McMullen, Wang, McKeon and Graham2019, Reference Shekar, McMullen, McKeon and Graham2020, Reference Shekar, McMullen, McKeon and Graham2021) or the centre mode (Page et al. Reference Page, Dubief and Kerswell2020; Dubief et al. Reference Dubief, Page, Kerswell, Terrapon and Steinberg2022) could lead to EIT in two dimensions. Recent advancements indicate that in the elastoinertial regime, wall mode structures dominate the dynamics (Beneitez et al. Reference Beneitez, Page, Dubief and Kerswell2024a
; Kumar & Graham Reference Kumar and Graham2024; Kumar et al. Reference Kumar, Constante-Amores and Graham2025), while at vanishing
$\textit{Re}$
, the dynamics tends to be organised around the centre mode structures (Morozov Reference Morozov2022; Beneitez et al. Reference Beneitez, Page, Dubief and Kerswell2024b
; Lewy & Kerswell Reference Lewy and Kerswell2025).
Recently, using spectral proper orthogonal decomposition (SPOD) (Towne et al. Reference Towne, Schmidt and Colonius2018), we have revealed the dominant coherent structures underlying the dynamics of EIT in 2-D channel flow (Kumar & Graham Reference Kumar and Graham2024). The SPOD is a temporal variant of proper orthogonal decomposition (POD) (Holmes et al. Reference Holmes, Lumley, Berkooz and Rowley2012), and decomposes time-dependent state variables into oscillating orthogonal basis functions such that the basis functions are coherent in both space and time. We demonstrated that the dynamics of EIT in channel flow is predominantly composed of a family of non-harmonic travelling waves, where the most dominant travelling wave of the family resembles the structure of the TS wave (Kumar & Graham Reference Kumar and Graham2024). The wall-normal velocity of the most dominant travelling wave has large-scale regular structures spanning the channel height, the streamwise velocity component has regular structures in the vicinity of the channel walls, and the polymeric stress field contains thin sheets of high polymeric stress fluctuations localised at the critical layers of the travelling wave. The critical layer is the location where the wave speed of the travelling wave matches the mean streamwise velocity. The remaining travelling waves of the family exhibit similar structures, but nested in the region bounded by the polymeric sheets, hence critical layers of the immediately slower travelling wave. The sheets of high polymeric stress act as barriers for flow and hence prevent velocity fluctuations from penetrating them (Kumar et al. Reference Kumar, Guasto and Ardekani2023). Thus the dynamics of EIT in channel flow originates from a nonlinear excitation of a wall mode, and is dominated by a family of self-similar nested travelling waves where the polymeric sheets of each travelling wave confine the velocity fluctuations of their immediate faster wave. Strengthening the notion of the connection of EIT with wall mode, it has also been found that the arrowhead structure resulting from the centre mode instability is disconnected from EIT, and it does not play a role in the self-sustenance of the dynamics of EIT (Beneitez et al. Reference Beneitez, Page, Dubief and Kerswell2024a ).
The numerical investigation of the dynamics of EIT in pipe flow has received relatively little attention (Lopez et al. Reference Lopez, Choueiri and Hof2019) and is the topic of the present study. Similar to channel flow, the centre mode becomes linearly unstable even in viscoelastic pipe flow at
$\textit{Re}=O(100)$
and
$\textit{Wi}=O(10)$
(Garg et al. Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018; Chaudhary et al. Reference Chaudhary, Garg, Subramanian and Shankar2021). Weakly nonlinear analysis around the linearly unstable centre mode demonstrates that nonlinearity stabilises the flow at high polymer concentration while it destabilises the flow at low polymer concentration, indicating supercritical and subcritical bifurcations at high and low polymer concentrations, respectively (Wan et al. Reference Wan, Sun and Zhang2021). In contrast, the wall mode in the pipe flow remains linearly stable at all
$\textit{Re}$
. However, using resolvent analysis, it has been found that the structure of the most strongly amplified mode by the linearised dynamics in viscoelastic pipe flow contains strong stress fluctuations localised at the critical layer (Zhang Reference Zhang2021), similar to the viscoelastic channel flow (Shekar et al. Reference Shekar, McMullen, Wang, McKeon and Graham2019), indicating a critical layer mechanism of the dynamics of EIT even in pipe flow. Further, there are common characteristics of the structures of EIT in both channel and pipe flows. For example, EIT in both channel and pipe flows exhibits near-wall vortical structures oriented in the direction perpendicular to the mean flow direction, and near-wall elongated streaks of streamwise velocity fluctuations aligned in the flow direction with a slight slope (Samanta et al. Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013; Lopez et al. Reference Lopez, Choueiri and Hof2019; Choueiri et al. Reference Choueiri, Lopez, Varshney, Sankar and Hof2021).
The basic characteristics of EIT in channel and pipe flows are fundamentally 2-D (Sid et al. Reference Sid, Terrapon and Dubief2018; Lopez et al. Reference Lopez, Choueiri and Hof2019). Therefore, in the present study, we investigate the dynamics of EIT in axisymmetric pipe flow. We perform direct numerical simulations of EIT using a spectral method. To analyse the dynamics, we introduce a viscoelastic variant of SPOD (Towne et al. Reference Towne, Schmidt and Colonius2018), which we refer to as ‘VESPOD’, that enables us to investigate the well-defined structures underlying EIT that evolve coherently in time. As detailed below, VESPOD is based on the total mechanical energy (kinetic plus elastic) of the flow, enabling an energy-based modal decomposition that contains both velocity and polymer stress fields. The eigenvalue spectrum of VESPOD exhibits distinct peaks, and by analysing those peaks, we find that the dynamics of EIT in pipe flow, in the parameter regime studied, is dominated by three distinct families of travelling waves. We also explore the effects of the domain length,
$\textit{Re}$
and
$\textit{Wi}$
on the dynamics of EIT and the travelling waves underlying its dynamics, and comment on the relation between EIT in pipe and channel flow. The remaining part of the paper has been organised as follows: § 2 discusses the methodology, § 3 focuses on results and discussion, and § 4 concludes with the findings of the study.
2. Formulation and governing equations
2.1. Direct numerical simulations of EIT
The flow is governed by the conservation of momentum and mass, described by the equations
\begin{equation} \frac {\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u} =-\boldsymbol{\nabla }\!p+\frac {\beta }{Re} {\nabla} ^2 \boldsymbol{u} + \frac {1-\beta }{Re} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau }_{\!p}+F(t)\,\boldsymbol{e}_z,\quad \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=0, \end{equation}
where
$\boldsymbol{u}$
and
$p$
represent the non-dimensional velocity field and pressure field, respectively. The Reynolds number is defined as
$\textit{Re}=\rho U_c R/\eta$
, where
$\rho$
and
$\eta$
are the density and zero shear rate viscosity of the fluid. The pipe radius (
$R$
) and the Newtonian laminar centreline velocity (
$U_c$
) are the characteristic length and velocity scales, respectively. The viscosity ratio
$\beta$
represents the contribution of solvent viscosity (
$\eta _s$
) to the zero shear solution viscosity, and is defined as
$\beta =\eta _s/\eta$
. The contribution of polymer chains to the stress tensor in the momentum equation is denoted as
$\boldsymbol{\tau }_{\!p}$
, and we use the FENE-P constitutive equation having polymer molecular diffusion to model it (Bird et al. Reference Bird, Curtiss, Armstrong and Hassager1987):
\begin{equation} \frac {\partial \boldsymbol{\alpha }}{\partial t}+\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\alpha } - \boldsymbol{\alpha }\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}-(\boldsymbol{\alpha }\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u})^{\rm T}=-\boldsymbol{\tau }_{\!p}+\frac {1}{ Re{Sc}} {\nabla} ^2 \boldsymbol{\alpha }, \end{equation}
\begin{equation} \boldsymbol{\tau }_{\!p}= \frac {1}{\textit{Wi}} \left ( \frac {\boldsymbol{\alpha }}{1-\mathrm{tr}(\boldsymbol{\alpha })/b}- \unicode{x1D644}\right )\!, \end{equation}
where
$\boldsymbol{\alpha }$
,
$ \unicode{x1D644}$
and
$b$
are the conformation tensor, the identity tensor and the maximum extensibility of the polymer chain, respectively. The Weissenberg number is defined as
$\textit{Wi}=\lambda U_c/R$
, where
$\lambda$
is the polymer relaxation time. We keep the molecular diffusion in the evolution equation of the conformation tensor (2.2) to ensure numerical stability. The strength of molecular diffusion is controlled by the Schmidt number
$\textit{Sc}=\eta /\rho D$
, where
$D$
is the diffusion coefficient.
We perform numerical simulations for the axisymmetric pipe flow EIT. We consider no-slip boundary conditions on the wall of the pipe (
$r=1$
), and axisymmetric boundary conditions at the centre of the pipe (
$r=0$
). The finite molecular diffusion of the polymer chains requires boundary conditions for the governing equation of the conformation tensor (2.2). At the wall (
$r=1$
), we solve (2.2) considering
$1/Sc=0$
, and use the respective value of
$\boldsymbol{\alpha }$
as the boundary condition. In the streamwise direction, we consider periodic boundary conditions. The volumetric flow through the pipe is kept constant at the Newtonian laminar value by tuning the body force (
$F(t)$
) at each time step.
We used Dedalus (Burns et al. Reference Burns, Vasil, Oishi, Lecoanet and Brown2020) to perform direct numerical simulations using the spectral method. The code validation can be found in Appendix A. We consider viscosity ratio
$\beta =0.97$
representing the extreme dilute limit, and
$b=6400$
corresponding to high-molecular-weight polymer chains (
$500$
kDa polyacrylamide). The length of the computational domain is
$L=5 \times n$
, where
$n$
is an integer. We used
$N_z \times n = 256 \times n$
Fourier basis functions in the streamwise direction, and
$N_r = 512$
Chebyshev basis functions in the radial direction. For time stepping, we use a second-order two-stage implicit–explicit Runge–Kutta scheme with fixed time step size
$\Delta t =0.001$
(unless stated otherwise). We consider molecular diffusion (
$\textit{Sc}=250$
) similar to previous studies (Sid et al. Reference Sid, Terrapon and Dubief2018; Kumar & Graham Reference Kumar and Graham2024). The simulations have been initiated with sufficiently large random perturbations in the state variables to trigger EIT. The simulation data corresponding to the initial 200 time units are dropped to make sure that the dynamics is fully developed. We also performed a sensitivity analysis of the EIT dynamics on
$\textit{Sc}$
and mesh resolution (Appendix B), and continue with the parameters mentioned above as they capture the essential physics with affordable computational cost.
2.2. The VESPOD
To analyse the EIT dataset for pipe flow, we use a viscoelastic variant of SPOD (VESPOD), which allows the modal decomposition of the velocity and polymeric stress fields together. We note that the flow displays translation symmetry in the streamwise (
$z$
) direction, in which case the Fourier basis functions are optimal in terms of energy (Holmes et al. Reference Holmes, Lumley, Berkooz and Rowley2012). To embed this symmetry in our VESPOD framework, first we Fourier-transform state variables in the streamwise direction, and then perform VESPOD decomposition wavenumber-by-wavenumber. The VESPOD represents the frequency-by-frequency VEPOD decomposition of a Fourier-transformed (in the temporal direction) time-dependent dataset (Holmes et al. Reference Holmes, Lumley, Berkooz and Rowley2012; Wang et al. Reference Wang, Graham, Hahn and Xi2014), where VEPOD is a viscoelastic variant of POD (Kumar et al. Reference Kumar, Constante-Amores and Graham2025). Hence VESPOD seeks to find the basis function
$\boldsymbol{\psi }(\kappa ,r,f)$
that maximises the objective function
\begin{equation} E \big\{ | \langle \boldsymbol{q}(\kappa ,r,f), \boldsymbol{\psi }(\kappa ,r,f) \rangle |^2 \big\}, \end{equation}
given the constraint
$\langle \boldsymbol{\psi }(\kappa ,r,f), \boldsymbol{\psi }(\kappa ,r,f) \rangle =1$
, where
$E\{{\cdot }\}$
and
$|{\cdot }|$
represent expectation and modulus operations on the data ensemble, respectively. The quantity
$\boldsymbol{q}(\kappa ,r,f)$
represents a vector containing state variables having scaled streamwise wavenumber
$\kappa$
(defined as the number of wavelengths per domain length) and temporal oscillation frequency
$f$
. In VESPOD,
$\boldsymbol{q}$
is defined as
\begin{equation} \boldsymbol{q}=\left [\boldsymbol{u}', \unicode{x1D64F}'\right ]\!, \end{equation}
where the stretch tensor is
$ \unicode{x1D64F}=\sqrt {( {1-\beta })/({\textit{Re}\textit{Wi}})}\,\boldsymbol{\theta }$
, and
$\boldsymbol{\theta } \boldsymbol{\cdot }\boldsymbol{\theta }=\boldsymbol{\alpha }/(1-\mathrm{tr}(\boldsymbol{\alpha })/b)$
(Wang et al. Reference Wang, Graham, Hahn and Xi2014). The superscript prime (
$'$
) represents the perturbation from the mean (averaged over
$t$
and
$z$
). The inner product
$\langle {\cdot }, {\cdot }\rangle$
in the objective function (2.4) has been defined as
\begin{equation} \langle \boldsymbol{q},\boldsymbol{\psi }\rangle = 2\pi \int _{0}^{1} \boldsymbol{q}(r) \boldsymbol{\cdot }\boldsymbol{\psi }(r)\, r {\rm d}r. \end{equation}
The specific definition of variable
$\boldsymbol{q}$
leads to the following inner product of
$\boldsymbol{q}$
with itself:
\begin{equation} \langle \boldsymbol{q},\boldsymbol{q}\rangle = 2\pi \int _{0}^1 \left \{\boldsymbol{u}' \boldsymbol{\cdot }\boldsymbol{u}' + \frac {1-\beta }{ Re \textit{Wi}}\boldsymbol{\theta }' : \boldsymbol{\theta }'\right \} r{\rm d}r. \end{equation}
This inner product and choice of
$\boldsymbol{q}$
are chosen because in the limit
$b\rightarrow \infty$
, (2.7) becomes the total mechanical energy of the flow and is thus a natural and interpretable inner product for analysis of viscoelastic flows. In this limit, VESPOD yields an optimal decomposition in terms of total mechanical energy. For
$b$
finite, the identification with total mechanical energy is approximate but is nevertheless useful. We also plot the probability distribution function (PDF) of
$\mathrm{tr}(\boldsymbol{\alpha })/b$
to find the distribution of stretching of polymeric chains in EIT (Appendix C). The probability of a large value of
$\mathrm{tr}(\boldsymbol{\alpha })/b$
is small, indicating that the inner product used in the present study provides a good representative of the total mechanical energy.
The maximisation of the objective function (2.4) leads to the self-adjoint eigenvalue problem
\begin{equation} 2\pi \int _{0}^1 E\{\boldsymbol{q}(\kappa ,r,f)\, \boldsymbol{q}^*(\kappa ,\rho ,f) \}\, \boldsymbol{\psi }(\kappa ,\rho ,f)\, \rho {\rm d} \rho =\sigma \boldsymbol{\psi }(\kappa ,r,f), \end{equation}
which gives an infinite set of eigenmodes
$\{\sigma _{\!j}, \boldsymbol{\psi }_{\!j}(\kappa ,r,f)\}$
, which are arranged in decreasing value of (approximate) mechanical energy. The state variables can be reconstructed using VESPOD modes as
\begin{equation} \tilde {\boldsymbol{q}}(\kappa , r, f)=\sum _{j=1}^{\infty }a_{\!j}(\kappa ,f)\,\boldsymbol{\psi }_{\!j}(\kappa , r, f), \end{equation}
where
$a_{\!j}(\kappa ,f)=\langle \boldsymbol{q}(\kappa , r, f), \boldsymbol{\psi }_{\!j}(\kappa , r, f) \rangle$
.
To obtain the VESPOD of a discrete time series, the data matrix containing the Fourier coefficients of the structures having wavenumber
$\kappa$
of the discrete time series of snapshots
$N_t$
is constructed as
\begin{equation} \unicode{x1D64C}_{\kappa }=[\boldsymbol{q}_{1,\kappa }, \boldsymbol{q}_{2,\kappa }, \ldots ., \boldsymbol{q}_{N_t, \kappa }], \end{equation}
where
$\boldsymbol{q}_{i,\kappa }=\boldsymbol{q}(\kappa ,r,t_i)$
. The computation of VESPOD using (2.8) requires an ensemble of independent realisations. In practice, when only a single time-resolved data point is available, but the flow is chaotic (adequate temporally decorrelated), the ensemble can be constructed by using data from different time windows of a very long time series (Welch Reference Welch1967). Therefore, multiple realisations of the flow time series have been generated by dividing the data matrix into overlapping blocks as
\begin{equation} \unicode{x1D64C}^{{{\kern1pt}}n}_{\kappa }=[\boldsymbol{q}_{1,\kappa }^n, \boldsymbol{q}_{2,\kappa }^n, \ldots,\boldsymbol{q}_{m,\kappa }^n,\ldots, \boldsymbol{q}_{N_f,\kappa }^n], \quad n=1,2,\ldots , N_b, \end{equation}
where
$N_b$
and
$N_f$
are the number of blocks and the number of snapshots in each block. For
$N_o$
overlapping snapshots in each block, the entries of the data matrices
$ \unicode{x1D64C}^{{{\kern1pt}}n}_{\kappa }$
and
$ \unicode{x1D64C}_{\kappa }$
are related as
$\boldsymbol{q}_{m,\kappa }^n=\boldsymbol{q}_{m+(n-1)(N_f-N_o), \kappa }$
. Next, we compute the discrete Fourier transform (DFT) of each block in the temporal direction. To minimise the spectral leakage resulting from the non-periodicity of the data in each block during the estimation of the DFT, we compute the DFT of the windowed data:
\begin{equation} \unicode{x1D64C}^{{{\kern1pt}}n,w}_{\kappa }=[w_1\boldsymbol{q}_{1,\kappa }^n, w_2\boldsymbol{q}_{2,\kappa }^n,\ldots,w_m\boldsymbol{q}_{m,\kappa }^n,\ldots, w_{N_f}\boldsymbol{q}_{N_f,\kappa }^n], \end{equation}
where
$w_m$
represents the nodal value of the symmetric Hamming window function and is given as
\begin{equation} w_m= 0.54-0.46\cos \left (\frac {2\pi (m-1)}{N_f-1} \right )\!. \end{equation}
The Fourier coefficients of the temporal DFT of
$ \unicode{x1D64C}^{{{\kern1pt}}n,w}_{\kappa }$
are given as
\begin{equation} \tilde { \kern-3pt\unicode{x1D64C}}^{{{\kern1pt}}n}_{\kappa }=[\tilde {\boldsymbol{q}}_{1,\kappa }^n, \tilde {\boldsymbol{q}}_{2,\kappa }^n, \ldots,\tilde {\boldsymbol{q}}_{m,\kappa }^n,\ldots, \tilde {\boldsymbol{q}}_{N_f,\kappa }^n], \end{equation}
where
$\tilde {\boldsymbol{q}}_{m,\kappa }^n$
represents the Fourier coefficient of wavenumber
$\kappa$
at frequency
$f_m$
in the
$n$
th block. Now the Fourier coefficients at frequency
$f_m$
are collected from each block, and a new data matrix
$\kern2pt\tilde { \kern-3pt\unicode{x1D64C}}_{m,\kappa }$
is constructed as
\begin{equation} \tilde { \kern-3pt\unicode{x1D64C}}_{m,\kappa }=[\tilde {\boldsymbol{q}}_{m,\kappa }^1, \tilde {\boldsymbol{q}}_{m,\kappa }^2,\ldots, \tilde {\boldsymbol{q}}_{m,\kappa }^n,\ldots, \tilde {\boldsymbol{q}}_{m,\kappa }^{N_b}]. \end{equation}
The VESPOD modes at frequency
$f_m$
can be obtained by calculating the eigenmodes of the discretised cross-spectral density matrix
$ \unicode{x1D64E}_{m\kappa }=\,\tilde { \kern-3pt\unicode{x1D64C}}_{m,\kappa }\,\tilde { \kern-3pt\unicode{x1D64C}}_{m,\kappa }^*$
by solving the eigenvalue problem
\begin{equation} \unicode{x1D64E}_{m,\kappa }\kern-2pt \unicode{x1D652} \boldsymbol{\psi }_{m,\kappa } = \boldsymbol{\psi }_{m,\kappa } \boldsymbol{\sigma }_{m,\kappa }, \end{equation}
where
$\boldsymbol{\psi }_{m,\kappa }=\boldsymbol{\psi }(\kappa ,r,f_m)$
represents the VESPOD mode structures at wavenumber
$\kappa$
and frequency
$f_m$
, and diagonal matrix
$\boldsymbol{\sigma }_{m,\kappa }$
represents mechanical energy associated with respective mode structures. The matrix
$ \unicode{x1D652}$
is a positive-definite weighting matrix, which properly accounts for the integration on a non-uniform discrete grid.
The eigenmode computation of (2.16) is expensive. To overcome this challenge, the following eigenvalue problem is solved:
\begin{equation} \tilde { \kern-3pt\unicode{x1D64C}}_{m,\kappa }^* \unicode{x1D652}\, \tilde { \kern-3pt\unicode{x1D64C}}_{m,\kappa } \boldsymbol{\varTheta }_{m,\kappa } = \boldsymbol{\varTheta }_{m,\kappa } \boldsymbol{\sigma }_{m,\kappa }. \end{equation}
This is faster to compute because the number of flow realisations (
$N_b$
) is much smaller than the number of grid points. This analogous problem gives the same non-zero eigenvalues as (2.16), and their eigenvectors are related as
\begin{equation} \boldsymbol{\psi }_{m,\kappa } = \,\tilde { \kern-3pt\unicode{x1D64C}}_{m,\kappa }\boldsymbol{\varTheta }_{m,\kappa } \boldsymbol{\sigma }_{m,\kappa }^{-1/2}. \end{equation}
The VESPOD mode structures having wavenumber
$\kappa$
and frequency
$f_m$
can be projected in the physical space as
\begin{equation} \boldsymbol{\psi }_m(z,r) = \boldsymbol{\psi }_{m,\kappa }\,{\rm e}^{{\rm i}2\pi \kappa z/L}. \end{equation}
We modify the MATLAB tool originally developed by Schmidt (Reference Schmidt2022) to perform VESPOD analysis of EIT in pipe flow. We used data corresponding to
$600$
time units sampled at time interval
$\Delta t_s=0.025$
time units, leading to a total of
$24\,000$
snapshots. The multiple realisations of the flow for most of the cases have been constructed by dividing the data into blocks having
$N_f=2000$
snapshots in each block with
$50\,\%$
overlap. This leads to a total of
$23$
blocks, which is sufficient for a converged computation of VESPOD. The current choices of
$\Delta t_s$
and
$N_f$
lead to frequency resolution
$\Delta f=1/(\Delta t_s\, N_f)=0.02$
. However, for
$\kappa =1$
structures, we have used
$\Delta f=0.01$
by considering
$N_f=4000$
to properly resolve distinct peaks in the eigenvalue spectra. A coarse frequency resolution may fail to resolve distinct spectral peaks, whereas an excessively fine frequency resolution can reduce statistical convergence and can generate numerous closely spaced peaks by resolving weak temporal variability within broadband structures, thereby making it difficult to isolate the dominant energetic structures. Therefore, an optimal value of
$\Delta f$
is chosen to balance frequency resolution and statistical convergence (Appendix D).
3. Results and discussion
We use VESPOD analysis to uncover the structures underlying the dynamics of EIT in axisymmetric pipe flow at several parameter regimes and domain lengths. First, we analyse dynamics at
$\textit{Re}=3000$
,
$\textit{Wi}=35$
and
$L=5$
(§ 3.1). Then we discuss dynamics of EIT in a longer domain (
$L=10$
) (§ 3.2) and the effect of
$\textit{Re}$
and
$\textit{Wi}$
on EIT (§ 3.3). Finally, we discuss the relationship between EIT in the pipe and channel flows, and their link to linear dynamics (§ 3.4).
Instantaneous (
$a$
) radial velocity, (
$b$
) streamwise velocity, and (
$c$
)
$zz$
-component of the stretch tensor of EIT in axisymmetric pipe flow at
$\textit{Re}=3000$
,
$\textit{Wi}=35$
and
$L=5$
. Means (averaged over
$z$
and
$t$
) of (
$d$
) velocity components and (e) stretch components. Dotted lines in (d,e) represent the laminar profiles at the same parameter.
3.1. The EIT at
$\textit{Re}=3000$
,
$\textit{Wi}=35$
and
$L=5$
Figure 1 depicts instantaneous snapshots of state variables and their temporal means for the dynamics of EIT at
$\textit{Re}=3000$
,
$\textit{Wi}=35$
and
$L=5$
. The dynamics of EIT has well-defined means (figures 1
d,e). Therefore, we report the perturbations of different state variables from their temporal means. The dynamics consists of weak velocity fluctuations, consistent with the nature of EIT (figures 1
a,b). The radial velocity field (
$u_r'$
) has distinct large-scale dynamic patterns dominating the region between the pipe centreline and the wall (figure 1(a), supplementary movie 1 are available at https://doi.org/10.1017/jfm.2026.11909). The streamwise velocity fluctuation (
$u_z'$
) has thin elongated streaks close to the wall; however, it contains large-scale structures in the vicinity of the pipe centreline (figure 1(b), supplementary movie 2). The polymeric stress field contains the formation of thin sheets of high polymeric stress in most of the pipe, with slight clustering near the pipe wall (Appendix B, supplementary movie 3). In the present paper, we use the scaled stretch field (
$ \unicode{x1D64F}$
) to represent the polymeric stress field (figure 1
$c$
), and due to the appropriate scaling of the stretch field, which is based on its contribution to the mechanical energy, the velocity and stretch fields have similar magnitude (figures 1
$a{-}c$
). This observation provides quantitative evidence that for EIT, both inertia and elasticity are equally important. The
$zz$
-component of the polymeric stress (and stretch) tensor dominates the stress field. Therefore, we report only
$T_{zz}'$
in the present study. The stretch field (
$T_{zz}'$
) contains thin inclined sheets of high stretch fluctuations situated well away from the pipe centre
$r=0$
(figure 1
$c$
).
The VESPOD energy spectrum of the structures having wavenumber
$\kappa =1$
has been depicted in figure 2(
$a$
). The leading eigenvalue (integrated over all frequencies) contributes approximately
$78 \,\%$
to the total variance (approximate mechanical energy), and there are three distinct peaks in the leading eigenvalue, suggesting that the mode structures corresponding to these peaks have a large contribution to the overall dynamics. In figure 3, we visualise the mode structures corresponding to different peaks in the energy spectrum. The radial velocity component has large-scale structures located away from the pipe centreline, which move from the location close to the pipe wall for the peak at low frequency to the vicinity of the pipe centreline for the peak at high frequency (figures 3
$a{-}c$
). The streamwise velocity component has large-scale structures centred at the pipe centreline, which shrink in the radial direction as the frequency of the peak increases; however, it has elongated patches of velocity fluctuations close to the wall (figures 3
$d{-}f$
). The stretch field contains a layer having thin inclined sheets of high stretch fluctuations, and the location of this layer moves to the centreline of the pipe as the frequency increases (figures 3
$g{-}i$
). These three structures are distinct travelling waves; their wave speeds
$c_w=fL/\kappa$
can be given as
$c_{w,1}$
,
$c_{w,2}$
and
$c_{w,3}$
in ascending order of their wave speed, and corresponding peaks in the energy spectrum as peak index 1, 2 and 3, respectively.
(
$a$
) Leading eigenvalue of the VESPOD energy spectrum of the structures having wavenumber
$\kappa =1$
underlying EIT at
$\textit{Re}=3000$
,
$\textit{Wi}=35$
and
$L=5$
; inset shows spectrum in log–log scale. (
$b$
) Leading eigenvalues of the VESPOD energy spectra for structures having different wavenumbers. Red symbols represent the peaks in the leading eigenvalues (circle, star and triangle are first, second and third peaks, respectively).
The VESPOD mode structures of (
$a{-}c$
) radial velocity (
$u_r'$
), (
$d{-}f$
) streamwise velocity (
$u_z'$
), and (
$g{-}i$
) the
$zz$
-component of the stretch tensor (
$T_{zz}'$
) corresponding to the peaks in the leading eigenvalue of the VESPOD energy spectrum of
$\kappa =1$
. Locations of critical layers have been denoted with white dotted lines in (
$g{-}i$
). Other parameters are
$\textit{Re}=3000$
,
$\textit{Wi}=35$
and
$L=5$
.
The concept of critical layer has proven useful in the understanding of transitional flows, as for travelling fluctuations it provides the most favourable place to exchange energy between the base flow and fluctuations (Drazin & Reid Reference Drazin and Reid1981). Therefore, we also explore the relationship between the critical layer of the travelling wave and the mode structures. As the mean velocity profile is close to the laminar parabolic profile (figure 1
$d$
), the location of the critical layer can be given as
$r_c=\sqrt {1-c_w}$
. The sheets of high polymeric stretch fluctuations in the VESPOD modes are localised around the location of the critical layer (figures 3
$g{-}i$
). This observation is consistent with the structures of travelling waves underlying EIT in channel flow (Shekar et al. Reference Shekar, McMullen, Wang, McKeon and Graham2019; Kumar & Graham Reference Kumar and Graham2024).
To explore the mechanism of the emergence of sheets in the polymeric stress field, in figure 4 we plot the streamlines of velocity fields that are constructed from the mean flow plus a small contribution from the peak VESPOD mode at the peaks in figure 2(
$a$
), in the reference frame moving with their respective wave speeds, along with
$T_{zz}'$
. There is a hyperbolic stagnation point in the velocity field (in the reference frame moving with the wave speed) at the critical layer of the travelling wave, where the flow is extension dominated. Polymeric chains get highly stretched as they approach and depart such a stagnation point, leading to the formation of thin polymeric sheets at the critical layer. This critical layer mechanism is consistent with the critical layer mechanism of the origin of polymeric sheets in channel flow EIT (Shekar et al. Reference Shekar, McMullen, Wang, McKeon and Graham2019).
Streamlines of a velocity field consisting of the mean flow and a small component of the peak VESPOD mode at (
$a$
)
$f=0.08\ (c_{w,1}=0.40)$
, (
$b$
)
$f=0.13\ (c_{w,2}=0.65)$
and (
$c$
)
$f=0.18\ (c_{w,3}=0.90)$
from figure 2(
$a$
) are shown in reference frames moving with their respective wave speeds, and are superimposed on the
$zz$
-component of the stretch tensor (
$T_{zz}'$
). The red circles denote hyperbolic stagnation points, and the streamlines passing through them are shown in white.
The VESPOD energy spectra of higher wavenumber structures have characteristics similar to those at
$\kappa =1$
. The leading eigenvalues contain most of the energy across the spectra of different wavenumber structures, and have three distinct peaks. In figure 2(
$b$
), we visualise the leading eigenvalues of VESPOD corresponding to different
$\kappa$
as a function of wave speed. For a given peak index, the wave speeds of the mode structures are nearly independent of wavenumber, implying that the higher wavenumber structures (
$\kappa \geqslant 2$
) are simple harmonics of the fundamental wavenumber structure (
$\kappa =1$
). This is also evident from figure 5. The travelling waves corresponding to each peak index belong to the respective unique family of travelling waves having unique wave speed (figure 5
$a$
), and the wave speed of the family of travelling waves increases almost linearly with the peak index (figure 5
$b$
).
Wave speeds of the VESPOD mode structures corresponding to the peaks in the energy spectra: (
$a$
) wave speed versus wavenumber, and (
$b$
) wave speed versus peak index. Parameter values are
$\textit{Re}=3000$
,
$\textit{Wi}=35$
and
$L=5$
.
The VESPOD mode structures of
$u_r'$
corresponding to different peaks in the leading eigenvalue of energy spectrum of wavenumbers (
$a{-}c$
)
$\kappa =2$
, (
$d{-}f$
)
$\kappa =3$
, (
$g{-}i$
)
$\kappa =4$
, and (
$j{-}l$
)
$\kappa =5$
. Other parameters are
$\textit{Re}=3000$
,
$\textit{Wi}=35$
and
$L=5$
.
In addition, we depict
$u_r'$
and
$T_{zz}'$
of mode structures corresponding to the three distinct peaks in the VESPOD energy spectra of higher wavenumber structures (
$\kappa \geqslant 2$
) in figures 6 and 7, respectively. The radial velocity component has large-scale structures, and these structures have a partition line oriented in the streamwise direction that divides the structures into two parts (figure 6). For each wavenumber, the location of the partition line moves to the centreline of the pipe as the wave speed increases. However, for a given wave speed (i.e. peak index), the location of the partition line is insensitive to the wavenumber. The stretch field is characterised by the formation of a layer, again localised at the critical layer position, having thin inclined sheets of high stretch fluctuations (figure 7). For a given wave speed, the locations of this layer for different wavenumber structures coincide with each other. However, they move to the pipe centreline as the wave speed increases for each wavenumber structure. Further, the location of the partition line of the structure in radial velocity coincides with the critical layer, and hence the polymeric sheets.
The discussion so far demonstrates that the dynamics of EIT in pipe flow at
$\textit{Re}=3000$
,
$\textit{Wi}=35$
and
$L=5$
is dominated by three distinct families of travelling waves, where higher wavenumber structures (
$\kappa \geqslant 2$
) are the simple harmonics of the respective fundamental waves (
$\kappa =1$
). There also exist some relationships among different families of travelling waves. For example, the sheets of high polymeric stretch fluctuations of
$c_{w,1}$
are located in the vicinity of the wall of the pipe, and the polymeric stretch sheets of
$c_{w,2}$
are confined in the region bounded by the pipe centreline and polymeric stretch sheets of
$c_{w,1}$
. Similarly, the sheets of
$c_{w,3}$
are bounded by the pipe centreline and the sheets of
$c_{w,2}$
. Thus the polymeric sheets of different families of travelling waves exhibit a nested arrangement, where the polymeric sheets of faster waves are confined by the polymeric sheets of immediately slower waves. A similar nesting relationship among the dominant travelling waves underlying EIT in channel flow has been observed (Kumar & Graham Reference Kumar and Graham2024).
The VESPOD mode structures of
$T_{zz}'$
corresponding to different peaks in the leading eigenvalue of energy spectrum of wavenumbers (
$a{-}c$
)
$\kappa =2$
, (
$d{-}f$
)
$\kappa =3$
, (
$g{-}i$
)
$\kappa =4$
, and (
$j{-}l$
)
$\kappa =5$
. Other parameters are
$\textit{Re}=3000$
,
$\textit{Wi}=35$
and
$L=5$
. These structures correspond to the velocity structures shown in figure 6.
3.2. The EIT at
$\textit{Re}=3000$
,
$\textit{Wi}=35$
and
$L=10$
To explore the effects of pipe length on EIT dynamics, we performed a simulation in a domain twice the length of the pipe discussed in § 3.1, but at the same
$\textit{Re}$
and
$\textit{Wi}$
. The velocity and stretch fields corresponding to EIT in the longer pipe are shown in figure 8. We analyse VESPOD spectra to identify the dominant coherent structures underlying them, and how pipe length affects them. The leading eigenvalue of the VESPOD spectrum of each wavenumber structure contains a dominant contribution to the total mechanical energy and has three distinct peaks, similar to the shorter pipe. We visualise the leading eigenvalue of different wavenumber structures as a function of wave speed, and find that the wave speeds corresponding to peaks in the leading eigenvalues are nearly independent of the wavenumber of structures (figure 9
$a$
). This shows that travelling waves with high wavenumber (
$\kappa \geqslant 2$
) are harmonics of the fundamental wave (
$\kappa =1$
). This is also evident from figure 9(
$b$
), where we see three distinct families of travelling waves originating from the simple harmonic excitation of respective fundamental waves. We have also visualised the structures of different travelling waves in figure 10. We emphasise that the travelling waves having an even wavenumber (
$\kappa =2,4, \ldots$
) in the pipe with
$L=10$
are the same travelling waves discovered in the pipe with
$L=5$
. However, the travelling waves with an odd wavenumber (
$\kappa =1,3,5, \ldots$
) in the pipe with
$L=10$
are new and do not exist in the pipe with
$L=5$
because they do not satisfy the periodicity in the shorter domain. We also note that the energy associated with
$\kappa =2$
structures is higher than that of
$\kappa = 1$
structures (at least for the first two families of travelling waves), suggesting that the structures having wavelength
$5$
units still dominate even in the pipe with
$L=10$
(figure 9
$a$
). This indicates that the pipe with
$L=5$
is long enough to capture the dominant structures.
Snapshots of (
$a$
) radial velocity, (
$b$
) streamwise velocity, and (
$c$
) the
$zz$
-component of the stretch tensor of EIT in axisymmetric pipe flow at
$\textit{Re}=3000$
,
$\textit{Wi}=35$
and
$L=10$
.
(
$a$
) Leading eigenvalues of the VESPOD energy spectra of different wavenumber structures at
$\textit{Re}=3000$
,
$\textit{Wi}=35$
and
$L=10$
. (
$b$
) Wave speeds of the travelling waves corresponding to the peaks in the leading eigenvalues of VESPOD energy spectra.
The VESPOD mode structures of
$u_r^{\prime }$
corresponding to different peaks in the leading eigenvalue of the energy spectrum of different wavenumber structures in the pipe of length
$L=10$
(
$\textit{Re}=3000$
and
$\textit{Wi}=35$
).
3.3. Effects of
$Re$
and
$\textit{Wi}$
To explore the effects of
$\textit{Wi}$
and
$\textit{Re}$
on the dynamics of EIT, we analyse EIT at
$\textit{Re}=3000,\ \textit{Wi}=70,\ L=5$
(figures 11
$a{-}c$
) and
$\textit{Re}=6000,\ \textit{Wi}=35,\ L=5$
(figures 11
$d{-}f$
). The VESPOD energy spectrum at
$\textit{Re}=3000,\ \textit{Wi}=70$
(figure 12
a) exhibits characteristics similar to those of the spectrum at
$\textit{Re}=3000,\ \textit{Wi}=35$
(figure 2
b). The three distinct harmonic families of travelling waves persist even at
$\textit{Wi}=70$
, which is evident from the three distinct peaks in the leading eigenvalues of the energy spectra of different wavenumber structures. For EIT at
$\textit{Wi}=70$
, the speeds of these three families of travelling waves are the same as for the waves at
$\textit{Wi}=35$
, suggesting that
$\textit{Wi}$
does not affect the speeds of different travelling waves at a fixed
$\textit{Re}$
. At higher value of
$\textit{Wi}$
, along with three harmonic peaks in the leading eigenvalues of VESPOD energy spectra, we also observe additional distinct peaks for some wavenumber structures (figure 12(
$a$
), purple diamond symbol). These are non-harmonic modes, having energy comparable to the harmonic modes. The mode structures corresponding to the peaks at
$\textit{Wi}=70$
are similar to those for
$\textit{Wi}=35$
. Therefore, here we only include the structures corresponding to the first three peaks in the leading eigenvalue of
$\kappa =1$
at
$\textit{Re}=3000$
and
$\textit{Wi}=70$
(figure 13). The mode structures shown in figures 13(
$a,d$
) and 13(
$c,f$
) correspond to the fundamental wave of the first (
$c_{w,1}$
) and second (
$c_{w,2}$
) families of harmonic travelling waves. However, the structures in figures 13(
$b,e$
) correspond to the non-harmonic peak in the leading eigenvalue of the
$\kappa =1$
structure at
$\textit{Re}=3000$
and
$\textit{Wi}=70$
(figure 12
$a$
). We do not see any fundamental difference between the mode structures of the harmonic and non-harmonic peaks.
Snapshots of (
$a,d$
) radial velocity, (
$b,e$
) streamwise velocity, and (
$c,f$
) the
$zz$
-component of the stretch tensor of EIT in axisymmetric pipe flow at (
$a{-}c$
)
$\textit{Re}=3000,\ \textit{Wi}=70$
and (
$d{-}f$
)
$\textit{Re}=6000, \textit{Wi}=35$
.
Leading eigenvalues of VESPOD energy spectra for different wavenumber structures at (
$a$
)
$\textit{Re}=3000,\ \textit{Wi}=70$
and (
$b$
)
$\textit{Re}=6000,\ \textit{Wi}=35$
. Red symbols represent harmonic peaks, while a purple diamond symbol represents a non-harmonic peak.
Further, we analyse dynamics at higher
$\textit{Re}$
(
$\textit{Re}=6000$
), while keeping
$\textit{Wi}=35$
, the same as in § 3.1 (figure 12
$b$
). Again, there exist three dominant families of harmonic travelling waves. However, the speeds of the travelling waves at
$\textit{Re}=6000$
are slightly different from those at
$\textit{Re}=3000$
. For example, the speed of the second family of travelling waves (
$c_{w,2}$
) at
$\textit{Re}=6000$
is higher than that at
$\textit{Re}=3000$
. At larger
$\textit{Re}$
, non-harmonic modes become increasingly important. As
$\textit{Wi}$
or
$\textit{Re}$
increases, the EIT dynamics becomes increasingly broadband and temporally complex, resulting in less distinct spectral peaks and weaker coherent travelling wave signatures. Accordingly, additional peaks emerge alongside the harmonic peaks at higher
$\textit{Re}$
and
$\textit{Wi}$
, indicating increasingly irregular nonlinear dynamics.
The VESPOD mode structures of (
$a{-}c$
)
$u_r^{\prime }$
and (
$d{-}f$
)
$T_{zz}^{\prime }$
corresponding to the first three peaks in the leading eigenvalue of
$\kappa =1$
at
$\textit{Re}=3000$
and
$\textit{Wi}=70$
(figure 12
$a$
). Here, (b,e) correspond to the structures of the non-harmonic peak (purple diamond in figure 12
$a$
); however, (a,d) and (c, f) represent harmonic peaks.
3.4. Relation to channel flow EIT and linear dynamics
The dynamics of EIT in channel and pipe flows share several similarities, but also exhibit important differences. Even though both dynamics have similar first glance features such as the large-scale structures in wall-normal velocity and the formation of thin inclined sheets of high stress in the polymeric stress field, the dominant coherent structures underlying them are distinct. The EIT in channel flow is dominated by a family of non-harmonic travelling waves, whereas in pipe flow it is dominated by the three distinct families of harmonic travelling waves. The wall-normal velocity component of these travelling waves exhibits large-scale structures centred at the centreline in channel flow, whereas they are located away from the centreline in pipe flow. The streamwise velocity component of these waves exhibits regular structures in the vicinity of the centreline in pipe flow, whereas in channel flow, such distinct structures are located away from the centreline. However, the polymeric stress field in both channel and pipe flows contains sheets of high stress fluctuations at the critical layers of the travelling waves. Further, the polymeric sheets of faster travelling waves in both channel and pipe flows are confined by the polymeric sheets of their immediately slower travelling waves, leading to a nested arrangement of polymeric sheets in both dynamics.
The resolvent analysis, which is based on the response of linearised governing equations to an external forcing, provides the modes that optimally describe the linear growth. It has been demonstrated that the SPOD modes are the same as the resolvent modes for uncorrelated resolvent-mode expansion coefficients, which is the case for the white-noise forcing of linear dynamics (Towne et al. Reference Towne, Schmidt and Colonius2018). The resolvent analysis for viscoelastic pipe flow considering the white noise forcing has found that the most strongly amplified resolvent mode exhibits thin tilted sheets like fluctuations in the polymeric stress field at the critical layer (Zhang Reference Zhang2021), which closely resembles the structure of the most dominant travelling wave discovered here using the VESPOD analysis, suggesting a possible connection of the dominant coherent structures underlying EIT in pipe flow with its linear dynamics. For channel flow, it has already been demonstrated that the TS mode resulting from the linear stability analysis closely resembles the most strongly amplified resolvent mode (Shekar et al. Reference Shekar, McMullen, Wang, McKeon and Graham2019) and the most dominant SPOD mode underlying EIT (Kumar & Graham Reference Kumar and Graham2024).
The differences between the velocity fields of the dominant travelling waves underlying EIT in channel and pipe flows result at least in part from the differences in their underlying Newtonian linear stability problems and symmetries imposed for 2-D dynamics. The linear instability of a wall mode leads to the emergence of the TS wave in Newtonian channel flow, which sets the structure and symmetry (shift-reflect) of the most dominant travelling wave underlying EIT in channel flow. In contrast, no analogous wall mode linear instability arises in pipe flow. Additionally, the constraint that the wall-normal velocity has to vanish at the axis for an axisymmetric flow precludes the shift-reflect symmetry of the TS wave for travelling waves underlying EIT in pipe flow. However, the common critical layer mechanism responsible for the formation of thin sheets in the polymeric stress field can explain the similarities observed in the structures of the polymeric stress field of the dominant travelling waves underlying EIT in both channel and pipe flows.
4. Conclusions
In the present study, we reveal the dominant well-defined structures underlying the dynamics of elastoinertial turbulence (EIT) in axisymmetric pipe flow using a viscoelastic variant of spectral proper orthogonal decomposition (VESPOD). The VESPOD enables us to perform modal decomposition of velocity and polymeric stress fields together, and identify coherently evolving mode structures that contribute the most to the total mechanical energy of the dynamics of EIT. The leading eigenvalue of the VESPOD spectrum contains most of the mechanical energy and hence dominates its dynamics. The leading eigenvalue of the VESPOD energy spectrum of each wavenumber (
$\kappa$
) structure exhibits three distinct peaks. Further analysis of the mode structures corresponding to these peaks reveals three distinct families of well-defined travelling waves, which are characterised by their wave speed, and dominate the dynamics of EIT in axisymmetric pipe flow. The higher wavenumber structures (
$\kappa \geqslant 2$
) of each family of travelling waves are harmonics of their respective fundamental waves (
$\kappa =1$
). The radial velocity component of these travelling waves contains large-scale structures spanning the radial direction; however, the polymeric stress field contains thin sheets of high stress fluctuations localised at the critical layers of the travelling waves. The critical layers contain hyperbolic stagnation points (in the reference frame moving with the wave speeds), which lead to the emergence of thin polymeric sheets. The substantial portion of fluctuations in EIT travel downstream with the three wave speeds as the form of highly coherent travelling wave structures, and the radial positions corresponding to their critical layers are the preferred locations for the formation of thin sheets in the polymeric stress field. Further, the polymeric sheets of a given travelling wave confine the polymeric sheets of immediately faster waves, leading to a nested arrangement of the polymeric sheets. At larger values of the Reynolds numbers and Weissenberg numbers, the EIT dynamics becomes more chaotic, leading to the emergence of non-harmonic travelling waves, along with the three families of harmonic travelling waves, having energy content comparable to the harmonic waves.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2026.11909.
Funding
This research was supported under Office of Naval Research grant N00014-18-1-2865 (Vannevar Bush Faculty Fellowship) and National Science Foundation grant CBET-2437151.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Validation of numerical tool
We use the method of the manufactured solution for the validation of the numerical tool. We assume a solution having velocity and conformation tensor fields
\begin{equation} u_z=(1-r^2)\sin(z), \quad u_r=0, \end{equation}
and
\begin{equation} \boldsymbol{\alpha }=(1+\sin(z))\unicode{x1D644}. \end{equation}
Next, we analytically obtain forcing functions to the governing equations (2.1) and (2.3) that satisfy the assumed solution (A1)–(A2). Using these forcing terms, we perform a numerical simulation, and find good agreement between the computed solution and the analytical solution (figure 14).
Exact (solid lines) and computed (dotted lines) profiles of
$u_z$
for a manufactured solution of viscoelastic flow through an axisymmetric pipe at
$\textit{Re}=3000$
,
$\textit{Wi}=35$
,
$\beta =0.97$
and
$b=6400$
(mesh resolution
$256 \times 512$
).
Appendix B. Mesh resolution and
$\textit{Sc}$
dependency of solution
To illustrate the impact of mesh resolution and
$\textit{Sc}$
on the EIT result, we depict the snapshots of the
$zz$
-component of the polymeric stress field (as this component is the most sensitive) obtained using different
$\textit{Sc}$
values and the respective mesh resolution and time step required to successfully evolve the EIT dynamics in axisymmetric pipe flow (figure 15). As the value of
$\textit{Sc}$
increases, polymeric stress diffusion becomes weaker, hence we see the emergence of even finer/sharper filaments oriented along the streamwise direction. The computation at higher
$\textit{Sc}$
becomes significantly expensive due to the requirement of higher mesh resolution and a smaller time step to successfully evolve the dynamics. However, it does not fundamentally change the essential features of EIT and hence the central conclusion of the study. Therefore, in the present study, we use
$\textit{Sc}=250$
,
$N_z \times N_r=256 \times 512$
and
$\Delta t=0.001$
, as these simulation parameters provide a reasonable balance between accuracy and computational cost.
Snapshots of the
$zz$
-component of the polymeric stress field for EIT in an axisymmetric pipe flow obtained using different mesh resolutions and
$\textit{Sc}$
at the physical parameter values
$\textit{Re}=3000$
,
$\textit{Wi}=35$
and
$L=5$
.
Appendix C. The PDF of
$\text{tr}(\boldsymbol{\alpha })/\boldsymbol{b}$
To estimate the deviation of mechanical energy from the exact value due to the finite
$b$
for the inner product used in the present study (2.7), we estimate the PDF of
$\text{tr}(\boldsymbol{\alpha })/b$
for the ensemble of polymeric stress field in EIT. We define the PDF as
\begin{equation} \text{PDF}(\tilde {p}) = \lim _{t\to \infty } \frac {\int _0^t \int _0^L \int _0^1 \delta (p(t,z,r)-\tilde {p})r\,{\rm d}r\,{\rm d}z\,{\rm d}t}{\int _0^t \int _0^L \int _0^1 r\,{\rm d}r\,{\rm d}z\,{\rm d}t}, \end{equation}
where
$p=\text{tr}(\boldsymbol{\alpha })/b$
,
$\tilde {p}$
is the value of
$p$
at which the PDF is evaluated, and
$\delta$
is the Dirac delta function. Although the maximum value of
$\text{tr}(\boldsymbol{\alpha })/b$
is
$\approx 0.9$
, the probability of a high value of
$\text{tr}(\boldsymbol{\alpha })/b$
is small (figure 16). In fact, the peaks in the PDF are at approximately
$\tilde {p}=0$
and
$\tilde {p}=0.4$
, suggesting that the inner product used in the present study is a good representative of the mechanical energy. Regardless of whether the inner product used in present study exactly captures the mechanical energy, it is an appropriate and physically motivated inner product and norm.
The PDF of
$\text{tr}(\boldsymbol{\alpha })/b$
for polymeric stress field in EIT at
$\textit{Re}=3000$
,
$\textit{Wi}=35$
,
$\beta =0.97$
and
$b=6400$
.
Appendix D. The VESPOD spectra based on different frequency resolutions
For wavenumber
$\kappa =1$
structures, the difference between the consecutive peak frequencies is small, therefore a higher frequency resolution is required to properly resolve the dominant peaks in the VESPOD spectra. We have plotted the leading VESPOD eigenvalue spectra for
$\kappa =1$
structures obtained using different frequency resolutions (figure 17). The frequency resolution
$\Delta f=0.02$
is too coarse to resolve all the peaks associated with
$\kappa =1$
structures, whereas
$\Delta f=0.005$
requires a larger dataset for statistical convergence, and also generates some additional peaks resolving weak variability within broadband structures. Therefore, we choose
$\Delta f=0.01$
for
$\kappa = 1$
structures to balance frequency resolution and statistical convergence, enabling clearer identification of the dominant energetic structures.
Leading eigenvalues of VESPOD spectra of
$\kappa =1$
structures obtained using different frequency resolutions at
$\textit{Re}=3000$
,
$\textit{Wi}=35$
and
$L=5$
. A total of 24 000 snapshots were used for the frequency resolutions
$\Delta f=0.02$
and
$\Delta f=0.01$
, whereas 36 000 snapshots were used for
$\Delta f=0.005$
to maintain statistical convergence.
a
b
c
zz
Re=3000
Wi=35
L=5
z
t
d
a
κ=1
Re=3000
Wi=35
L=5
b
a−c
ur′
d−f
uz′
g−i
zz
Tzz′
κ=1
g−i
Re=3000
Wi=35
L=5
a
f=0.08 (cw,1=0.40)
b
f=0.13 (cw,2=0.65)
c
f=0.18 (cw,3=0.90)
a
zz
Tzz′
a
b
Re=3000
Wi=35
L=5
ur′
a−c
κ=2
d−f
κ=3
g−i
κ=4
j−l
κ=5
Re=3000
Wi=35
L=5
Tzz′
a−c
κ=2
d−f
κ=3
g−i
κ=4
j−l
κ=5
Re=3000
Wi=35
L=5
a
b
c
zz
Re=3000
Wi=35
L=10
a
Re=3000
Wi=35
L=10
b
ur′
L=10
Re=3000
Wi=35
a,d
b,e
c,f
zz
a−c
Re=3000, Wi=70
d−f
Re=6000,Wi=35
a
Re=3000, Wi=70
b
Re=6000, Wi=35
a−c
ur′
d−f
Tzz′
κ=1
Re=3000
Wi=70
a
a
uz
Re=3000
Wi=35
β=0.97
b=6400
256×512
zz
Sc
Re=3000
Wi=35
L=5
tr(α)/b
Re=3000
Wi=35
β=0.97
b=6400
κ=1
Re=3000
Wi=35
L=5
Δf=0.02
Δf=0.01
Δf=0.005
Open access
