2.1. Governing equations

We consider the linear stability of steady fully developed flow of a viscoelastic fluid in a rigid circular pipe of radius $R$ as shown in figure 1. A cylindrical polar coordinate system is used with $r$, $\theta$ and $z$ denoting the radial, azimuthal and axial directions, respectively. The following scales are used for non-dimensionalizing the governing equations: radius of the pipe $R$ for lengths, maximum base-flow velocity $U_{max}$ for velocities, $R/U_{max}$ for time and $\rho U_{max}^2$ for pressure and stresses, with $\rho$ being the density of the fluid.

Figure 1.

Schematic diagram showing the geometry and the coordinate system considered.

The governing (non-dimensional) continuity and Cauchy momentum equations are given by

(2.1a,b)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v} = 0, \quad \frac{{\partial} \boldsymbol{v}}{{\partial} t} + (\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{v} ={-}\boldsymbol{\nabla} p +\frac{\beta}{Re}\nabla^2\boldsymbol{v} + \boldsymbol{\nabla}\boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}. \end{equation}

Here, $\boldsymbol {v}$ is the fluid velocity field, $p$ is the pressure field and ${\boldsymbol{\mathsf{T}}}$ is the polymeric contribution to the stress tensor, which in turn is given by the Oldroyd-B constitutive relation (Larson Reference Larson1988) as follows:

(2.2)\begin{equation} W \left(\frac{{\partial} {\boldsymbol{\mathsf{T}}}}{{\partial} t}+\left( \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} \right){\boldsymbol{\mathsf{T}}}-{\boldsymbol{\mathsf{T}}} \boldsymbol{\cdot} \left(\boldsymbol{\nabla} \boldsymbol{v}\right) -\left(\boldsymbol{\nabla} \boldsymbol{v} \right)^{\textrm{T}} \boldsymbol{\cdot}{\boldsymbol{\mathsf{T}}}\right) +{\boldsymbol{\mathsf{T}}} = \frac{(1-\beta)}{Re} \{\boldsymbol{\nabla}\boldsymbol{v}+(\boldsymbol{\nabla}\boldsymbol{v})^\textrm{T}\}. \end{equation}

The solvent to solution viscosity ratio is denoted by $\beta = \eta _s/\eta$, where the solution viscosity is $\eta = \eta _p + \eta _s$, $\eta _s$ and $\eta _p$ being the solvent and polymer viscosities, respectively; $\beta = 0$ and $1$ denote the UCM and Newtonian limits. For a fixed $\beta$, the dimensionless groups relevant to the stability of the Oldroyd-B fluid described previously are the Reynolds number $Re=\rho U_{max} R/\eta$ and the Weissenberg number $W=\lambda U_{max} /R$, which is a ratio of the polymer relaxation time $\lambda$ to the flow time scale. The Oldroyd-B model describes the stress in a dilute solution of polymer chains modelled as non-interacting Hookean dumbbells (Larson Reference Larson1988), and is invariably the first model used in the examination of elastic phenomena involving dilute polymer solutions. Consistent with the aforementioned microscopic picture, the Oldroyd-B model assumes the relaxation time to be independent of both the shear rate and the polymer concentration. As the model predicts a shear-rate-independent viscosity, the non-Newtonian (elastic) effects in this model arise from an effective tension along the streamlines (arising from flow-aligned dumbbells), which manifests as a shear-rate-independent first normal stress different in viscometric flows. This model has been used extensively, and with considerable success, in earlier investigations of inertialess elastic instabilities in flows with curved streamlines (Larson, Shaqfeh & Muller Reference Larson, Shaqfeh and Muller1990; Pakdel & McKinley Reference Pakdel and McKinley1996; Shaqfeh Reference Shaqfeh1996). The so-called Boger fluids constitute an experimental realization of this constitutive model (Boger & Nguyen Reference Boger and Nguyen1978). As discussed later in the manuscript (in § 4.4.1, where we use scaling arguments in the context of the FENE-P model to assess the role of shear thinning), whereas shear thinning can play an important role especially in flow through microtubes (Samanta et al. Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013; Chandra et al. Reference Chandra, Shankar and Das2018), the Oldroyd-B model does have the necessary ingredients to qualitatively predict the instabilities observed in experiments.

2.2. Base state

The base-state velocity profile is the classical Hagen–Poiseuille profile because the nonlinear terms in the upper-convected derivative of the polymer shear stress $T_{rz}$ are identically zero. The non-dimensional base flow velocity vector is given by

(2.3)\begin{equation} \bar{\boldsymbol{v}} = \begin{bmatrix} \bar{v}_{r}\\ 0\\ \bar{v}_{z} \end{bmatrix}= \begin{bmatrix} 0\\0\\U(r)\end{bmatrix}, \end{equation}

where $U(r)=1-r^2$ for pipe Poiseuille flow. Here, and in what follows, base state quantities are denoted by an overbar. The polymer contribution to the stress tensor in the base state is given by

(2.4)\begin{equation} \bar{{\boldsymbol{\mathsf{T}}}}=\begin{bmatrix} \bar{\tau}_{rr} & 0 & \bar{\tau}_{rz}\\ 0 & \bar{\tau}_{\theta\theta} & 0\\ \bar{\tau}_{zr} & 0 & \bar{\tau}_{zz}\end{bmatrix}= \frac{(1-\beta)}{Re} \begin{bmatrix} 0 & 0 & U'\\ 0 & 0 & 0\\ U' & 0 & 2WU'^2\end{bmatrix},\end{equation}

where, $f'\equiv \mathrm {D}f \equiv {\mathrm {d}f}/{\mathrm {d}r}$. Unlike the velocity profile, the base-state stress profile differs from that of a Newtonian fluid in having a tension along the streamlines proportional to the square of the velocity gradient.

2.3. Linear stability analysis

A temporal linear stability analysis is carried out wherein the base-state described previously is subjected to small-amplitude axisymmetric perturbations. Owing to the absence of a Squire-like theorem for pipe flow even in the simpler Newtonian case, in general, both axisymmetric and non-axisymmetric disturbances need to be considered for viscoelastic Oldroyd-B fluids. However, for the parameter regime probed in this study, we find axisymmetric disturbances alone to be unstable, and this study is therefore restricted to axisymmetric disturbances. The total velocity, pressure and stress are expressed in terms of their base-state values and perturbations as

(2.5a)\begin{gather} \boldsymbol{v} = \boldsymbol{\bar{v}} + \boldsymbol{\hat{v}}, \end{gather}

(2.5b)\begin{gather}p = \bar{p} + \hat{p}, \end{gather}

(2.5c)\begin{gather}{\boldsymbol{\mathsf{T}}} = \bar{{\boldsymbol{\mathsf{T}}}} + \hat{{\boldsymbol{\mathsf{T}}}}, \end{gather}

with $\hat {f}$ denoting the perturbation to the dynamical quantity $f$. For axisymmetric disturbances without swirl (i.e. $\hat {v}_{\theta } = 0$), the perturbation velocity and stress tensor are

(2.6a,b)\begin{equation} \hat{\boldsymbol{v}}= \begin{bmatrix} \hat{v}_{r}\\ 0\\ \hat{v}_{z} \end{bmatrix},\quad \text{and}\quad \hat{{\boldsymbol{\mathsf{T}}}}= \begin{bmatrix} \hat{\tau}_{rr} & 0 & \hat{\tau}_{rz}\\ 0 & \hat{\tau}_{\theta\theta} & 0\\ \hat{\tau}_{zr} & 0 & \hat{\tau}_{zz} \end{bmatrix}. \end{equation}

Next, the perturbation quantities above are represented in the form of Fourier modes in the axial ($z$) direction in the following manner:

(2.7)\begin{equation} \hat{f}(r,z;t) = \tilde{f}(r)\mathrm \exp\{\mathrm{i}k (z - c t)\}, \end{equation}

where $k$ is the axial wavenumber and $c = c_r + i c_i$ is the complex wave speed. The flow is temporally unstable (stable) if $c_i > 0$ $(< 0)$. Substituting (2.7) in the linearized versions of (2.1a,b)–(2.2), we obtain the following set of linearized governing equations:

(2.8)\begin{gather} \left(\mathrm{D} + \frac{1}{r}\right)\tilde v_{r}+\mathrm{i} k\tilde v _{z}=0, \end{gather}

(2.9)\begin{gather}G \tilde{v}_r ={-}\mathrm{D}\tilde{p} + \left[\left(\mathrm{D}+\frac{1}{r}\right)\tilde{\tau}_{rr}+\mathrm{i}k\tilde{\tau}_{rz}-\frac{\tilde{\tau}_{\theta\theta}}{r}\right] + \frac{\beta}{Re}\mathrm{L}\tilde{v}_r, \end{gather}

(2.10)\begin{gather}G \tilde{v}_z +U'\tilde{v}_r={-}\mathrm{i}k\tilde{p} + \left[\left(\mathrm{D}+\frac{1}{r}\right)\tilde{\tau}_{rz}+\mathrm{i}k\tilde{\tau}_{zz}\right] + \frac{\beta}{Re}\left(\mathrm{L}+\frac{1}{r^2}\right)\tilde{v}_z, \end{gather}

(2.11)\begin{gather}H\tilde{\tau}_{rr}=2\frac{(1-\beta)}{Re}\left(\mathrm{D}+W\mathrm{i}kU'\right)\tilde{v}_r, \end{gather}

(2.12)\begin{gather} H\tilde{\tau}_{rz}-WU'\tilde{\tau}_{rr} = \frac{(1-\beta)}{Re} [\{\mathrm{i}k-W(U''-U'\mathrm{D}-2\mathrm{i}kWU'^2)\} \tilde{v}_r +(\mathrm{D}+W\mathrm{i}kU')\tilde{v}_z ], \end{gather}

(2.13)\begin{gather} H\tilde{\tau}_{\theta\theta} = 2\frac{(1-\beta)}{Re} \frac{\tilde{v}_r}{r}, \end{gather}

(2.14)\begin{gather}H\tilde{\tau}_{zz}-2WU'\tilde{\tau}_{rz} = 2\frac{(1-\beta)}{Re} \left[{-}2W^2U'U'' \tilde{v}_r +\{\mathrm{i}k+WU'\left(\mathrm{D}+2W\mathrm{i}kU'\right)\}\tilde{v}_z\right], \end{gather}

where $G = \mathrm {i}k(U-c)$, $H = 1 + WG$ and $\mathrm {L} = (\mathrm {D}^2 + {\mathrm {D}}/{r}-{1}/{r^2}-k^2)$. The no-slip boundary conditions $\tilde {v}_r = 0$ and $\tilde {v}_{z} = 0$ are applicable at $r = 1$, whereas at $r = 0$, the conditions $\tilde {v}_r = 0$ and $\tilde {v}_{z} = \text {finite}$, corresponding to regularity of axisymmetric disturbances in the vicinity of the centreline, are used (Batchelor & Gill Reference Batchelor and Gill1962; Khorrami, Malik & Ash Reference Khorrami, Malik and Ash1989).

2.4. Numerical method

We use two independent formulations to solve the viscoelastic eigenvalue problem for the wavespeed $c$. In the first, the governing equations for perturbation stresses (2.11)–(2.14) are substituted in (2.9)–(2.10) to obtain two linearized ordinary differential equations corresponding to the momentum balances in $r$- and $z$-directions in addition to (2.8), and the dependent variables in this formulation are $\tilde {v}_r$, $\tilde {v}_z$ and $\tilde {p}$. In the second formulation, we directly solve the system of linear equations (2.8)–(2.14), with $\tilde {v}_r/r$ as the dependent variable instead of $\tilde {v}_r$, with the other variables being $\tilde {v}_z$, $\tilde {p}$, $\tilde {\tau }_{\theta \theta }$, $\tilde {\tau }_{rr}$, $\tilde {\tau }_{rz}$ and $\tilde {\tau }_{zz}$. The simplified equations represent a homogeneous eigenvalue problem, and are solved using the standard spectral collocation numerical scheme based on Chebyshev polynomials (Boyd Reference Boyd2000; Trefethen Reference Trefethen2000). Results from the two different spectral approaches show excellent agreement. Further, the eigenvalues obtained from the spectral method were verified using a shooting method (Ho & Denn Reference Ho and Denn1977; Lee & Finlayson Reference Lee and Finlayson1986) implemented for the first formulation, based on an adaptive step size Runge–Kutta integrator and a Newton–Raphson procedure for determining the eigenvalue. The integration for the shooting method was carried out from a point near the centreline $r = \epsilon$ (with $\epsilon \rightarrow 0$) to the pipe wall at $r = 1$. The velocities at $r = \epsilon$ were obtained using a Frobenius series expansion (Garg & Rouleau Reference Garg and Rouleau1972) about the regular singular point $r = 0$. The shooting method gives very accurate (based on our choice of tolerance, typically $10^{-9}$) results when sufficiently close initial guesses are provided, whereas the number of polynomials $N$ required for convergence of eigenvalues in the spectral method depends mainly on the nature of the eigensolutions and the parameter values. Typically, the $N$ required for convergence of eigenvalues for finite $\beta$ is in the range $150$–$200$, whereas that for the UCM limit ($\beta \rightarrow 0$) is in the range $400$–$500$. There is no prior literature that reports the eigenspectrum for pipe flow of an Oldroyd-B fluid, and hence our numerical procedure was benchmarked in the Newtonian limit (obtained by setting $W = 0$ or $\beta = 1$). Results in this limit are available, for instance, in Schmid & Henningson (Reference Schmid and Henningson1994, Reference Schmid and Henningson2001).

3.1. Spectra at fixed $\beta$ and different $E$

Figures 2(a) and 2(b) show the eigenspectra for Newtonian pipe flow at $Re = 6000$ and $Re = 600$, respectively. The spectrum has the well-known ‘Y’-shaped structure (Mack Reference Mack1976; Schmid & Henningson Reference Schmid and Henningson2001) at $Re = 6000$, but this is only beginning to form in the spectrum at $Re = 600$. The Y-shaped structure is composed of three branches: (i) the ‘A branch’ corresponding to ‘wall modes’ with $c_r \rightarrow 0$ for $Re \gg 1$ on the top left; (ii) the ‘P branch’ that consists ‘centre modes’ with $c_r \rightarrow 1$ for $Re \gg 1$ on the top right and (iii) the ‘S branch’ that consists modes with $c_r \approx 2/3$ extending down to $c_i = -\infty$. For $Re \gg 1$, the decay rates of the least-stable centre and wall modes vary as $|c_i| \sim Re^{-1/2}$ and $|c_i| \sim Re^{-1/3}$, respectively (Meseguer & Trefethen Reference Meseguer and Trefethen2003), implying that the centre modes are the least stable at large $Re$. Although these scalings need not necessarily hold for the moderate $Re$ ($= 600$) considered in figures 3 and 4, the centre mode is nevertheless found to be the least stable. Consistent with previous studies (Schmid & Henningson Reference Schmid and Henningson1994, Reference Schmid and Henningson2001), all modes for Newtonian pipe flow are found to be stable. We discuss the nature of the least-stable mode in more detail in § 3.3.

Figure 2.

The ‘Y’-shaped eigenspectrum for Newtonian pipe flow subjected to axisymmetric disturbances for $k = 3$, and for (a) $Re = 6000$ and (b) $Re = 600$.

Figure 3.

Eigenspectra for pipe flow of an Oldroyd-B fluid at $\beta = 0.8$, $Re = 600$ and $k = 3$, and for different $E$ in the range $5 \times 10^{-4}$–$10^{-3}$: (a) $E = 10^{-4}$; (b) $E = 2\times 10^{-4}$; (c) $E = 3\times 10^{-4}$; (d) $E = 4\times 10^{-4}$; (e) $E = 5\times 10^{-4}$ and (f) $E = 10^{-3}$. The eigenspectra are obtained for $N=200$, and there is excellent convergence of the spectra for $N=200$ and $250$ (not shown). An elliptical-ring structure is prominent at the lower $E$, but is absent beyond $E = 10^{-3}$. The vertical locations of the CS1 and CS2 lines and the Newtonian spectrum for the same $Re$ and $k$ are shown for reference.

Figure 4.

Eigenspectra for the Oldroyd-B fluid for different $E$ in the range $0.002$–$0.1$, at $\beta = 0.8$, $Re = 600$, $N = 200$ and $k = 3$: (a) $E = 0.002$; (b) $E = 0.003$; (c) $E = 0.006$; (d) $E = 0.01$; (e) $E = 0.05$ and (f) $E = 0.1$. The spectra, shown for a narrower range of $c_r$ and $c_i$ compared with figure 3, demonstrate how the discrete centre modes (labelled $2$, $3$ and $4$) merge into and emerge out (labelled $2'$, $3'$ and $4'$) of the CS as $E$ is increased. The least-stable centre mode (labelled 1) always stays above the CS, and eventually becomes unstable at $E = 0.1$. The vertical locations of the CS lines and the Newtonian spectrum at the same $Re$ and $k$ are shown for reference.

The spectra for pipe flow of an Oldroyd-B fluid reduce to the Newtonian one either when $E \rightarrow 0$ (at fixed $\beta$) or when $\beta \rightarrow 1$ (at fixed $E$). We therefore first examine the effect of viscoelasticity as $E$ is increased from zero, with $\beta$ fixed at $0.8$. Figures 3 and 4 show the viscoelastic eigenspectra for $Re =600$, with $E$ ranging from $5 \times 10^{-4}$ to $0.1$. The values of $\beta$ and $Re$ are chosen so they are close to the experimental conditions of Samanta et al. (Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013) and our earlier theoretical work (Garg et al. Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018). With increasing $E$, figures 3 and 4 show the classical Y-shaped structure of the Newtonian spectrum to be altered by elasticity in a singular manner. There are important differences between the two spectra even for the smallest $E$, the most prominent of these being the appearance of two CS for the viscoelastic case, similar to viscoelastic plane shear flows (Renardy & Renardy Reference Renardy and Renardy1986; Sureshkumar & Beris Reference Sureshkumar and Beris1995; Graham Reference Graham1998; Wilson, Renardy & Renardy Reference Wilson, Renardy and Renardy1999; Grillet et al. Reference Grillet, Bogaerds, Peters and Baaijens2002; Chokshi & Kumaran Reference Chokshi and Kumaran2009). It is now well understood (Graham Reference Graham1998; Chaudhary et al. Reference Chaudhary, Garg, Shankar and Subramanian2019; Subramanian, Reddy & Roy Reference Subramanian, Reddy and Roy2020) that the CS arise from the local nature of the constitutive model for the polymeric stress, and disappear when non-local diffusive effects are incorporated in the constitutive relation (see § 4.3). The eigenvalues corresponding to the CS are obtained by setting to zero the coefficient of the highest-order derivative in the differential equation governing the stability. This coefficient turns out to be the product $[1 + \mathrm {i}kW(U-c)] [1 + \beta \mathrm {i}kW(U-c)]$, which leads to a pair of horizontal ‘lines’ in the $c_r$–$c_i$ plane with $c_i = -1/(kW)$ and $c_i = -1/(\beta k W)$, and with $0 \leq c_r \leq 1$. Henceforth, these two CS are respectively abbreviated as ‘CS1’ and ‘CS2’ respectively, with CS1 being present even in the limit of a UCM fluid, and CS2 being present only when there is a solvent contribution ($\beta \neq 0$), receding to $c_i = -\infty$ in the limit $\beta \rightarrow 0$.

Figure 3 explores the spectra for the smallest $E$ (ranging from $10^{-4}$ to $10^{-3}$), the range of $c_r$ and $c_i$ being chosen so as to provide a larger view of the spectra. Here, in addition to the modified Y-shaped structure of the Newtonian spectra and the two CS lines, there exist a class of modes that form a ‘ring’ that surrounds the CS at small $E$ of $O(\sim 10^{-4})$. Similar to CS1 and CS2 , all modes belonging to the ring structure are stable for the range of $E$ explored. For small $E$, the modes on the ring appear to be symmetrically distributed (figures 3a to 3c) about the S branch, forming an approximate ellipse. As $E$ is increased, the modes move towards the CS with the ring decreasing in size. For $E = 4\times 10^{-4}$, these modes move closer, intermingling with the other modes that emerge from the CS (figure 3d), and the ring structure is now fully distorted. At still higher $E\sim 10^{-3}$ (see figures 3e and 3f), the modes originally on the ring collapse, wrapping around the CS in an irregular manner. To understand the origin of the ring structure, it is relevant to recall a prominent feature of the viscoelastic spectra (at non-zero $Re$) in the UCM limit ($\beta = 0$): an infinite sequence of discrete modes corresponding to damped shear waves in a viscoelastic fluid (discussed in § 3.2), and are referred to as the high-frequency Gorodtsov–Leonov (‘HFGL’) modes (Gorodtsov & Leonov Reference Gorodtsov and Leonov1967; Kumar & Shankar Reference Kumar and Shankar2005; Chaudhary et al. Reference Chaudhary, Garg, Shankar and Subramanian2019). This sequence corresponds to $c_i = -1/(2kW)$, and extends to infinity in either direction parallel to the $c_r$ axis. As we demonstrate in figure 10, at any finite $\beta$, the infinite-in-extent HFGL line curves downwards, eventually meeting the S branch, and thereby leading to the aforementioned ring structure for sufficiently small $E$.

In figure 4, we explore the spectra for larger $E$, in the range $0.002$–$0.1$, with the ranges of $c_r$ and $c_i$ being chosen to provide a more magnified view of the spectra. As $E$ is increased, the vertical locations of the two CS move up towards $c_i = 0$, and in the process, the discrete ‘elastic’ centre modes (labelled 2, 3 and 4 in figure 4a; and that lie above the CS) disappear into the CS. As $E$ is increased further, new discrete elastic centre modes (now shown with the labels $2'$, $3'$ and $4'$ in figures 4c and 4d, which lie below the CS) emerge out of the CS. The labelling of the modes that emerge below the CS are for the purposes of reference only, and there is no connection between these modes (with primes) and those (without primes) that disappeared into the CS. This was ascertained from the absence of any resemblance between the eigenfunctions of the modes that disappear into and reappear from the CS. A more detailed account of the evolution of the centre modes as $E$ is varied is provided in figure 7. Importantly, the least-stable centre mode (labelled 1) does not merge into the CS, and always stays above it. As $E$ is increased to $0.1$, this centre mode becomes unstable, and corresponds to the instability first reported by Garg et al. (Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018). This scenario of the unstable centre mode being a smooth continuation of its stable Newtonian counterpart is, however, sensitive to $\beta$, and we show (and provide further details in § 3.1.1) that there exist other parameter regimes where the centre mode that eventually becomes unstable emerges out the CS with increasing $E$, and there is no connection to the least-stable Newtonian centre mode. Other new stable centre modes (with $c_r \rightarrow 1$; see figures 4d–4f) and wall modes (with $c_r \rightarrow 0$, and even negative; see figures 4e and 4f), which have no Newtonian counterparts, appear below the CS with increasing $E$. Increase in $E$ has a stabilizing effect on these modes. The aforementioned annihilation and creation of discrete modes with increase in $E$ occurs because both the CS are branch cuts (Wilson et al. Reference Wilson, Renardy and Renardy1999) for Poiseuille flow. Note that, for plane Couette flow, only CS2 is a branch cut. It is well known that discrete eigenmodes can appear or disappear out of the branch cut as parameters are varied, and this aspect is discussed further in § 3.1.1 in the specific context of the centre mode.

Figure 5 shows the spectra in the near-Newtonian limit of $\beta = 0.96$ and for $E$ ranging over the interval $(0.4,4)$, overlaid in a single plot, in order to demonstrate the variation of not just the (eventually) unstable centre mode, but also of the other stable modes. For the higher $E$ considered in figure 5, the two CSs lie very close to $c_i = 0$ (and to each other for the chosen $\beta$), and the modes in the Newtonian P branch have therefore already disappeared into the CS, with new modes emerging from below. Thus, the trajectories of the modes shown in figure 5(a) are for the modes that start off below the CS. The enlarged version in figure 5(b) shows the spectra in terms of the scaled growth rate $kWc_i$, which fixes the vertical location of both the CS (for fixed $\beta$), and allows one to focus on the trajectory of the unstable centre mode with varying $E$. The continuous curve indicating the trajectory of the centre mode, as $E$ is varied, is obtained using the shooting method with much finer increments in $E$. This figure shows that the centre mode first emerges out of the CS, in the form of a bump in the continuous spectrum balloon, at $E \approx 0.6$, and becomes unstable as $E$ is increased to $0.712$. The centre mode remains unstable for $0.712 \leq E \leq 2.5$, but becomes stable for $E > 2.5$, with $|c_i|$ eventually scaling as $1/E$ for large $E$. Thus, figures 4 and 5(b) show that there are two qualitatively different trajectories of the unstable centre mode with increasing $E$. For the lower $\beta$ $(=0.8)$, the centre mode appears as a smooth continuation of the least-stable Newtonian centre mode, whereas for $\beta = 0.96$, it emerges from the continuous spectrum, with no obvious connection to the Newtonian spectrum. This aspect is discussed in more detail in § 3.1.1.

Figure 5.

(a) Eigenspectra for different values of $E$ for $\beta =0.96$, $Re=500$ and $k=1$. (b) Enlarged version of the region in (a) near the unstable centre mode $c_r = 1$. The scaled growth rate $kWc_i$ fixes the vertical location of both the CS (for $\beta = 0.96$, CS1 and CS2 lie very close to each other). The continuous line for the trajectory of the unstable centre mode is obtained using the shooting method.

Figures 6(a)–6(d) show the velocity ($v_r$ and $v_z)$ and stress ($T_{rz}$ and $T_{zz}$) eigenfunctions, for different $E$, corresponding to a few of the unstable centre modes shown in figure 5(b). The velocity and $T_{rz}$ eigenfunctions are largely insensitive to variations in $E$, but the $T_{zz}$ eigenfunction shows a distinct and sharp peak for the smaller $E$ ($= 0.9$) near the radial location where the phase speed of the disturbances equals the local base flow velocity. Although the amplitudes of the axial velocity eigenfunctions in figure 6 are larger near the central core region of the pipe, the disturbance fields are nevertheless spread across the entire pipe cross-section for the parameters considered. As shown in § 4.2, only for sufficiently large $Re$ ($>$1000) does the localization of the velocity eigenfunctions near the centre become prominent. It is worth emphasizing this feature here because recent studies (Shekar et al. Reference Shekar, McMullen, Wang, McKeon and Graham2019) have inaccurately characterized the centre mode instability, analysed in Garg et al. (Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018) and the present work, as always being localized in the vicinity of the centreline regardless of $Re$.

Figure 6.

Velocity and polymer stress eigenfunctions corresponding to the unstable centre modes (in figure 5b) at different $E$ for $\beta = 0.96$, $Re = 500$ and $k=1$: (a) radial velocity; (b) axial velocity; (c) $rz$ polymer stress; and (d) $zz$ polymer stress.

3.1.1. The origin of the centre mode at fixed $\beta$ and varying $E$

The origin of the centre mode is more clearly demonstrated in figure 7(a) through the variation of $c_i$ with $E$ for the first four least-stable modes from the Newtonian P branch, obtained using the shooting method. For $\beta = 0.8$, consistent with the spectra in figure 4, the least-stable Newtonian centre mode always lies above the CS (figure 7a), smoothly continuing with increasing $E$, eventually becoming unstable for $E \approx 0.1$ (shown later in inset (A) of figure 9a). However, the other (more) stable Newtonian centre modes (labelled $2$, $3$ and $4$ in figure 4a) vanish into CS1 as $E$ is increased, and new modes appear out of CS1 with further increase in $E$, subsequently suffering a second jump across the CS2 line. The modes that emerge out of CS2 were those identified as $2'$, $3'$ and $4'$ in the spectra in figure 4. This feature is also evident in the variation of the phase speeds with $E$ in figure 7(b). An analogous phenomenon was reported by Chokshi & Kumaran (Reference Chokshi and Kumaran2009) for the least-stable wall mode in plane Couette flow of an Oldroyd-B fluid.

Figure 7.

Effect of increasing $E$ on the first four least-stable centre modes of Newtonian origin for $\beta = 0.8$, $Re = 600$ and $k = 3$. EI: elasto-inertial. (a) Growth rate versus $E$, with the inset presenting a magnified view of the jumps suffered by the individual modes as they cross CS1 ($c_i = -1/(kW)$) and CS2 ($c_i = -1/(\beta kW)$). (b) Phase speeds corresponding to the modes shown in (a).

In figure 8(a), we examine the effect of increasing $E$ in the UCM and near-UCM limits at fixed $Re$ and $k$ (note that, regardless of the value of $\beta$, $E = 0$ corresponds to the Newtonian limit). For $\beta = 0$, the decay rate of the least-stable Newtonian centre mode decreases with increasing $E$, even to the point of reducing to $\sim 2 \times 10^{-4}$ at $E \approx 8 \times 10^{-3}$ (about 1/100th of the decay rate in the Newtonian limit), but the mode remains stable. As elastic effects are responsible for the unstable centre mode, it might be expected that this instability should persist even in the absence of solvent contribution to the stress. The eigenspectra for pipe flow of a UCM fluid were computed for a vast range of parameters $0.5 < k < 3$, $100<Re<20\,000$ and $0<E<1$. Unlike the spectrum for plane channel flow of a UCM fluid (Sureshkumar & Beris Reference Sureshkumar and Beris1995; Chaudhary et al. Reference Chaudhary, Garg, Shankar and Subramanian2019), only stable modes were obtained for pipe flow of a UCM fluid subjected to axisymmetric disturbances. Thus, as originally stated in Garg et al. (Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018), the centre mode instability in viscoelastic pipe flow requires the combined effects of both the polymer elasticity and solvent viscous effects, in addition to fluid inertia.

Figure 8.

Effect of increasing $E$ on the least-stable Newtonian centre mode for UCM and Oldroyd-B fluids. (a) Plots of $c_i$ for the least-stable centre mode at $Re = 6000$ for $\beta = 0$ and $0.2$, with the inset (A) showing the enlarged region near the unstable range of the mode and inset (B) showing the corresponding phase speeds. (b) Scaled growth rate of elasto-inertial modes at $Re = 6000$, $\beta = 0.92$ and $Re = 500$, $\beta = 0.96$, and the inset showing the corresponding phase speeds.

For $\beta = 0.2$ (see figure 8a), however, the least-stable Newtonian centre mode does become unstable for $0.01 < E < 0.011$. For both $\beta$, the centre mode trajectory is similar to that shown in figure 7, in that it remains above the CS over the range of $E$ examined. The corresponding phase speeds (inset (B) of figure 8a), for both $\beta = 0$ and $0.2$, show a weak non-monotonic behaviour with $E$, although $c_r \leq 1$ for all $E$. The contrasting behaviour for $\beta$ close to unity (representing dilute solutions) is shown in figure 8(b). The main figure shows the variation of the scaled growth rate $kWc_i$ with $E$ for two different sets of (near-unity) $\beta$ and $Re$. The instability occurs at significantly larger values of $E \sim O(1)$, in contrast to figure 8(a), and the unstable range of $E$ is also larger. In contrast to the trend for $\beta = 0.2$, the continuation of the Newtonian centre modes for both $\beta = 0.92$ and $0.96$ collapses into the CS at smaller $E$ (not shown). It is the trajectories of the new discrete modes, that emerge from the CS at slightly larger $E$, and that become unstable for $E \sim O(1)$, that are shown in figure 8(b). The corresponding phase speeds for $\beta = 0.92$ and $0.96$ are shown in the inset of figure 8(b).

In figures 7 and 8, we show two different trajectories for the centre mode, as a function of $E$, depending on $\beta$. In order to clarify the change in the nature of the centre mode trajectory, from a continuous variation of $c_i$ with increasing $E$ at smaller $\beta$, to a discontinuous variation for near-unity $\beta$, figure 9(a) shows the behaviour of the centre mode for $\beta = 0.8$, $0.82$ and $0.85$. The centre mode trajectory remains above CS1 until instability, for both $\beta = 0.8$ and $0.82$, whereas for $\beta = 0.85$, the centre mode disappears into the CS at $E \approx 0.009$ (inset (A) of figure 9a). Thus, in this case, there exists a range $0.009\lessapprox E \lessapprox 0.024$ where the centre mode does not exist. This range, which extends from the point of encounter of this mode with CS1 to the point of emergence of the new mode from CS1 at higher $E$, varies with increasing $\beta$. Evidently, the critical $\beta$, below which the centre mode is a smooth continuation of the least-stable Newtonian centre mode, lies somewhere between $0.82$ and $0.85$ (for $Re = 600$ and $k = 3$). Note that, despite the discontinuous transition in terms of the collapse into the CSs, the interval of instability in $E$ varies smoothly with increasing $\beta$ (the inset (B) in figure 9a). Figure 9(b) shows the corresponding phase speeds, and the enlarged region in the inset shows that the trend for $c_r$ versus $E$ curves is more or less same for $\beta = 0.8$, $0.82$ and $0.85$, except for $\beta =0.85$, where the absence of the mode in the interval $0.009\lessapprox E \lessapprox 0.024$ leads to a gap in the $c_r$ curve.

Figure 9.

Effect of increasing $E$ on the viscoelastic centre mode for fixed values of $\beta = 0.8$, $0.82$ and $0.85$, $Re = 600$ and $k = 3$. (a) Growth rates, with inset (A) showing the enlarged region over the range of $E$ for which the centre mode is discontinuous (marked by vertical dotted lines) owing to CS1 with $c_i=-1/(kW)$. The locations of CS2 (with $c_i = -1/(\beta k W)$) are not shown owing to the closely separated values of $\beta$ used in this figure. Inset (B) shows the region where the centre mode is unstable for all the three values of $\beta = 0.8$, $0.82$ and $0.85$. (b) Phase speeds corresponding to the modes shown in (a). Inset in (b) shows the enlarged region near $E \rightarrow 0$.

3.2. Spectra at fixed $E$ and different $\beta$

We next examine the viscoelastic eigenspectra as $\beta$ is increased from zero, at fixed $E$. We begin with figure 10 which illustrates the effect of increasing $\beta$, starting from the UCM spectrum, at $E = 10^{-4}$. A moderate $Re$ ($=600$) is chosen in order to keep the spectral features relatively simple, requiring only a modest resolution (the number of collocation points $N$), and thereby allowing us to focus on the large-scale features. Figure 10 illustrates the singular feature of the bending down of the HFGL line for non-zero $\beta$. The bending down can be interpreted as a (very strong) stabilization of these modes owing to the solvent viscosity. For the larger $\beta$ ($0.5$ and $0.8$), the bending is ‘complete’, leading to the ring-like structure within the range of $c_i$ examined; this then clarifies the origin of the structure seen before in figure 3. Figure 11 shows spectra at a higher $Re$ ($=6000$), for different $\beta$, and with $E = 0.01$. The spectrum for $\beta = 0$ (figure 11a) now has a more intricate structure, necessitating an enlarged view of the phase speed interval $(0,1)$. The features of the high-$Re$ UCM channel-flow spectrum were first explained in Chaudhary et al. (Reference Chaudhary, Garg, Shankar and Subramanian2019), and include CS1 (which appears as a balloon owing to the finite resolution), the HFGL and additional discrete modes with $c_r \in [0,1]$ that lie on either side of the HFGL line, roughly along the contours of an ‘hourglass’. These features of the UCM pipe-flow spectrum, in figure 11(a), are analogous to the channel flow case described previously.

Figure 10.

Eigenspectra for $E = 10^{-4}$, $Re = 600$ and $k = 3$ demonstrating the bending down of the HFGL modes with increasing $\beta$, leading to the appearance of a ring-like structure for the higher $\beta$. The vertical location of CS1 is independent of $\beta$, whereas that of CS2 varies with $\beta$.

Figure 11.

Unfiltered viscoelastic eigenspectra for $E=0.01$, $Re=6000$, $k=2$ and different $\beta$: (a) $\beta = 0$, UCM; (b) $\beta = 1\times 10^{-4}$; (c) $\beta = 5\times 10^{-4}$; (d) $\beta = 5\times 10^{-3}$; (e) $\beta = 5\times 10^{-2}$; and (f) $\beta = 0.2$. The decay rate of the least-stable centre mode in the UCM limit decreases with increase in $\beta$ and the centre mode eventually becomes unstable at $\beta = 0.2$.

Figures 11(b)–11(d) show the spectra for $\beta$ in the range $10^{-4}$–$5\times 10^{-3}$. Figure 11(b) shows that even the smallest $\beta$ has a profound effect on the HFGL modes. In contrast to the UCM spectrum at $Re = 600$ (figure 11a), where the bending of HFGL line became evident only for $c_r$ well outside the base-state interval, the bending down of the HFGL modes is evident at $Re = 6000$ even for $c_r \in (0,1)$; see figures 11(b) and 11(c). The bent HFGL line has all but disappeared as $\beta$ is increased to $10^{-3}$ (figure 11c), again demonstrating that the HFGL modes are rapidly damped by small amounts of solvent viscosity. Owing to this drastic stabilization even at rather small $\beta$, the HFGL modes in the original UCM spectrum become irrelevant to the parametric regimes (corresponding to relatively dilute solutions, with $\beta \sim 0.6$ and higher) explored later in this study. Further, the ‘density’ of stable modes present in the hourglass structure in the UCM limit also decreases rapidly as $\beta$ is increased from zero, with the hourglass structure virtually absent for $\beta = 0.05$. Most importantly, although almost all other modes in the hourglass structure of the UCM spectrum are rapidly stabilized with increasing $\beta$ (figures 11d–11f), the least-stable centre mode (with $c_r\approx 1$ and $c_i \rightarrow 0$) is rather unaffected by the small increase in $\beta$. In fact, as shown in figures 11(e) and 11(f), for the largest $\beta$ shown ($\beta = 0.2$), the centre mode becomes unstable. Thus, as originally stated in figure 8(a), it appears that all three effects, namely, elasticity, solvent viscous stresses and fluid inertia are important ingredients for the instability of the centre mode in viscoelastic pipe flow.

In figure 12, we show the eigenspectra (overlaid) as $\beta$ is reduced from unity, again at fixed $E$, $Re$ and $k$; note that the $\beta$ shown are all higher than the threshold value for collapse into the CS (the analogue of that identified in figure 9a, but for $Re = 6000$). Figure 12(a) is for $\beta$ close enough to unity that the centre mode has not emerged out of the CS yet (the other stable modes, with $c_r \rightarrow 0$, are not shown). Thus, the trends in this figure pertain to all other (least-stable) modes on the P branch. Figure 12(a) shows no discernible trend in the behaviour of the P branch modes with changing $\beta$. For instance, as $\beta$ is decreased, the least-stable Newtonian mode moves in the clockwise sense in $(c_r,c_i)$-plane. In contrast, the mode LSCM2 smoothly continues from a Newtonian mode at the junction of the ‘APS’ structure present at $\beta = 0$. The remaining modes are, however, smooth continuations of the modes of the Newtonian P branch, but these move in the counter-clockwise sense with decreasing $\beta$. Eigenspectra for smaller $\beta$ in the interval $0.85\leq \beta \leq 0.98$ are shown in figure 12(b), the focus being on the centre mode. The centre mode first emerges at $\beta = 0.96$, and becomes unstable for $\beta \in [0.88,0.945]$. The smooth (blue) curve, passing through the spectral centre mode eigenvalues, shows the trajectory of the centre mode with decreasing $\beta$, obtained using the shooting method. Thus, at a fixed $E$, the unstable centre mode always emerges out of the CS as $\beta$ is decreased from unity.

Figure 12.

Viscoelastic pipe-flow eigenspectra for $E=0.15$, $Re=6000$, $k=1$ and for varying $\beta$: (a) $(1-\beta )=0$ to $0.02$; and (b) $(1-\beta )=0.02$ to $0.15$. Panel (a) shows the region in the vicinity of the ${\rm P}$ branch. Panel (b) focuses on the trajectory of the centre mode that becomes unstable for $\beta \in (0.88,0.945)$.

The velocity and stress eigenfunctions corresponding to some of the unstable centre modes in figure 12(b) are shown in figure 13. For the higher $Re$ ($= 6000$), the axial velocity eigenfunctions are more localized near the centre (as $\beta$ approaches unity), compared with those for $Re = 600$ shown in figure 6(b). The axial stress $T_{zz}$ (figure 13d) shows a sharp peak as $\beta$ approaches unity at fixed $E$, similar to the feature that was seen earlier (figure 6d), albeit with decreasing $E$ at a fixed $\beta$ for $Re = 600$. For the values of $\beta$ examined here, both the axial and radial eigenfunctions exhibit a rather smooth variation with $r$, unlike the rapid, oscillatory variation (not shown) characteristic of wall modes ($c_r \rightarrow 0$) for $\beta \rightarrow 0$. The latter are analogous to wall modes in viscoelastic channel flow whose structures was examined in detail by Chaudhary et al. (Reference Chaudhary, Garg, Shankar and Subramanian2019) (see figure 20 therein); the overall similarity of the pipe and channel flow UCM spectra was discussed previously.

Figure 13.

Velocity and polymer stress eigenfunctions corresponding to the unstable centre modes (in figure 12b) at different $\beta$ for $E = 0.15$, $Re = 6000$ and $k=1$: (a) radial velocity; (b) axial velocity; (c) $T_{rz}$ stress; and (d) $T_{zz}$ stress.

3.2.1. The origin of the centre mode at fixed $E$ and varying $\beta$

While the discussion pertaining to figure 12(b) showed that centre mode emerges out of the CS as $\beta$ is decreased from unity, in figure 14(a), we address the question of what happens as $\beta$ decreases to zero (the UCM limit). This figure shows the variation of the scaled growth rate, $kWc_i$, of the centre mode with varying $\beta$ at fixed $Re$, $E$ and $k$. As $\beta$ decreases from unity, the centre mode emerges out of CS1 (when $kWc_i = -1$) at a critical $\beta$, and becomes unstable as $\beta$ decreases further. The critical $\beta$ corresponding to the emergence of the centre mode is closer to unity for higher $E$. The range of unstable $\beta$ also approaches unity for larger $E$, while also narrowing down in extent, with a concomitant increase in the growth rate. A similar narrowing down occurs when $E$ approaches the lower threshold for the instability, for the chosen $\beta$ (the blue curve in figure 14a). Figure 14(b) shows that the corresponding $c_r$ remains close to (and less than) unity for the entire range of $\beta$.

Figure 14.

(a) Scaled growth rate and (b) phase speed, for the centre mode, with varying $\beta$, for different fixed sets of parameters $(E,Re, k) = (0.01, 6000, 2), (0.02, 6000, 2), (0.05, 6000, 2), (0.15, 6000, 1)$ and $(1, 500, 1)$.

In figure 15, we show the three possible behaviours, within the parameter regimes explored, for the trajectory of the least-stable centre mode as $\beta$ is varied from the UCM to the Newtonian limit. For the smallest elasticities (e.g. $E = 0.005$ in figure 15a), when the two CS are highly stable, and well outside the range of $c_i$ shown, the centre mode, while remaining stable, smoothly continues all the way from the UCM limit ($\beta = 0$) to the Newtonian ($\beta = 1$) limit without suffering any discontinuities or abrupt endings. For moderate elasticities (e.g. $E = 0.015$ shown in figure 15b), the $c_i$ versus $\beta$ curve for the least-stable centre mode starts from the Newtonian end ($\beta = 1$), but abruptly ends as it encounters CS2 from below. On the other hand, the least-stable centre mode in the UCM limit continues to finite $\beta$, abruptly ending at the location of its encounter with CS1 from above. Corresponding phase speeds for $E = 0.015$ are shown in figure 15(c), with the inset showing an enlarged view near $\beta =1$, where the variation of the phase speed $c_r$ with $\beta$ is quite sharp. For the chosen parameters, the centre mode still remains stable for all $\beta$. Finally, for higher elasticity (e.g. $E = 0.15$), the $c_i$ versus $\beta$ curve for the least-stable mode from the Newtonian end continues all the way up to the UCM limit without suffering discontinuities as shown in figure 15(d), ending up as a centre mode in the UCM spectrum. Inset (A) shows a magnified view of the sharp variation of the $c_i$ curve near $\beta =0$. The least-stable centre mode in the UCM limit behaves similar to the previous case of $E=0.015$, with an abrupt ending as it collapses onto CS1 from above, the only difference now being that the mode is unstable for a small range of $\beta$ (owing to the higher $E$); inset (B) provides the enlarged view of the unstable range of $\beta$. The corresponding phase speeds for $E=0.15$ are shown in figure 15(e) with inset (A) showing the enlarged view of the non-monotonic behaviour near the Newtonian limit. Note that the Newtonian centre mode does not suffer a jump despite crossing the CS2 curve ($c_i = -1/(\beta kW)$) in figure 15(d); this is only an apparent crossing since, as shown in figure 15(e), its phase speed exceeds unity, and it therefore ‘goes around’ CS2 with decreasing $\beta$. In contrast, the discontinuities in the centre mode trajectory, in figure 15(b), occur because $0 < c_r < 1$.

Figure 15.

The three possible centre mode trajectories with variation in $\beta$ for $Re=600$, $k=3$ at various fixed values of the elasticity number: $E = 0.005$ in (a), $E = 0.015$ in (b,c), and $E = 0.15$ in (d,e). Insets $(A)$ in (d) and (e) show the enlarged regions near $\beta =0$ and $\beta = 1$, respectively. Inset $(B)$ in panels (d) show the enlarged views of the unstable range of $\beta$. (a) Growth rate and phase speed; (b) growth rate; (c) phase speed; (d) growth rate; and (e) phase speed.

3.3. Centre versus wall modes in viscoelastic pipe and channel flows

In the results presented so far, we have characterized the behaviour of the elasto-inertial centre mode as a function of $E$ and $\beta$. Although this mode may either be directly related to a Newtonian centre mode (for $\beta$ below a threshold), or be disconnected from the Newtonian spectrum (for $\beta$ above), the interpretation is nevertheless that the EIT observed in recent experiments (Samanta et al. Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013; Chandra et al. Reference Chandra, Shankar and Das2018; Choueiri et al. Reference Choueiri, Lopez and Hof2018) is the outcome of a linear instability associated with this centre mode. In sharp contrast to this picture, in a recent effort, Shekar et al. (Reference Shekar, McMullen, Wang, McKeon and Graham2019) have argued based on DNS and a singular-value decomposition analysis that EIT in channel flow might instead be closely related to the elastically modified TS mode. As is well known, the TS mode is the least-stable wall mode in the Newtonian limit, and this remains true for the range of elasticities considered by the authors. Thus, the premise of Shekar et al. (Reference Shekar, McMullen, Wang, McKeon and Graham2019) continues to be along the lines of a subcritical bifurcation to EIT, similar in spirit to the earlier efforts of Meulenbroek et al. (Reference Meulenbroek, Storm, Bertola, Wagner, Bonn and van Saarloos2003) and Morozov & van Saarloos (Reference Morozov and van Saarloos2005, Reference Morozov and van Saarloos2007) in the inertialess limit, and to the work of Stone et al. (Reference Stone, Waleffe and Graham2002), Stone & Graham (Reference Stone and Graham2003), Stone et al. (Reference Stone, Roy, Larson, Waleffe and Graham2004) and Li & Graham (Reference Li and Graham2007) based on an elastic modification of 3D ECS structures. The main difference is that the bifurcation ascribed by Shekar et al. (Reference Shekar, McMullen, Wang, McKeon and Graham2019) is supposedly to a finite amplitude 2D mode, with EIT-like dynamics. The authors reported results for $Re = 1500$ (where the Newtonian flow is turbulent), $\beta = 0.97$, and for $0 < W < 50$. It is worth noting that, for these parameters, the elastically modified ECSs originally examined by Graham and co-workers (Li & Graham Reference Li and Graham2007) also exist, although Shekar et al. (Reference Shekar, McMullen, Wang, McKeon and Graham2019) restrict themselves to two-dimensional initial conditions.

While the present study is restricted to linear (modal) stability of pipe flow of an Oldroyd-B fluid, it is nevertheless instructive to compare the viscoelastic pipe and channel flow spectra in order to assess the relative importance of centre and wall modes in these geometries. Such an assessment would help set the template (in terms of the relevant linear modes, both discrete and continuous) for a nonlinear bifurcation analysis. We show representative eigenspectra for pipe flow (in figure 16) for a range of $E$ that subsumes the range ($0 < E < 0.013$) considered by Shekar et al. (Reference Shekar, McMullen, Wang, McKeon and Graham2019), for $\beta = 0.97$, $Re = 1500$ and $k=0.4{\rm \pi}$. In each panel, the corresponding Newtonian spectrum is also shown for comparison (as open blue circles). For $0.0005 <E < 0.05$, with increasing $E$, discrete modes collapse into the CS, and new ones emerge from below, similar to what was shown in figure 3. The centre mode emerges from CS1 at $E = 0.46$, (see inset of figure 16d), becoming unstable at $E \approx 0.5$ (figure 16e). Importantly, for the parameters considered in figure 16, the centre mode always remains the least-stable or least-unstable mode. This feature remains true even for other regimes investigated in this study ($Re \in 100$–$2000$). This is unlike Newtonian channel flow, where there is a range of parameters where the wall mode (i.e. the TS mode) is the least stable (or unstable), and this remains true for small but finite $E$.

Figure 16.

Viscoelastic pipe-flow eigenspectra at $\beta = 0.97$ for different $E$: (a) $0.0005$, (b) $0.005$, (c) $0.05$, (d) $0.46$ and (e) $0.5$, compared with the Newtonian eigenspectrum ($E =0$), for $Re = 1500$ and $k=0.4{\rm \pi}$.

An important feature of Newtonian pipe flow is the absence of a critical-layer singularity (Drazin & Reid Reference Drazin and Reid1981; Schmid & Henningson Reference Schmid and Henningson2001) for axisymmetric disturbances, as a result of which there is no axisymmetric analogue of the two-dimensional TS instability. This difference between pipe and channel flows appears to persist even in the presence of elasticity. In figure 17, we show, via contour plots, the spatial structure of the least-stable centre and wall modes marked in figure 16(b). Further, and in sharp contrast to viscoelastic channel flow, where the elastically modified TS mode was shown to have the $T_{xx}$ (the stream-wise component of the normal stress) eigenfunction strongly localized in the critical layer (see figure 2 of Shekar et al. Reference Shekar, McMullen, Wang, McKeon and Graham2019), neither the least-stable centre nor the wall mode in pipe flow exhibits a comparably strong localization of $T_{zz}$; in fact, the extent of localization is more stronger for the centre mode. For these reasons, the connection between the (stable) TS wall mode to the elasto-inertial structures suggested by Shekar et al. (Reference Shekar, McMullen, Wang, McKeon and Graham2019) (in the context of viscoelastic channel flow) is not applicable for viscoelastic pipe flow. This aspect will be discussed in more detail in a future communication (Khalid et al. Reference Khalid, Chaudhary, Garg, Shankar and Subramanian2021), where we show that, even for viscoelastic channel flows, the parameter regime relevant to the proposed TS-mode-based subcritical mechanism of Shekar et al. (Reference Shekar, McMullen, Wang, McKeon and Graham2019) is somewhat restricted. It is worth emphasizing that all of the experiments on viscoelastic transition (with the exception of Srinivas & Kumaran Reference Srinivas and Kumaran2017) pertain to the pipe geometry. Further, and importantly, recent simulations in both the channel (Samanta et al. Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013; Sid et al. Reference Sid, Terrapon and Dubief2018) and pipe (Lopez et al. Reference Lopez, Choueiri and Hof2019) geometries have found analogous (span-wise oriented) coherent structures, suggesting a common underlying mechanism for elasto-inertial transition.

Figure 17.

Contours of the radial ($\hat {v}_r$), axial ($\hat {v}_x$) velocities and axial normal stress ($\hat {\tau }_{zz}$) in the $r$–$z$ plane for the least-stable wall and centre modes (marked in figure 16b) in pipe flow for $E = 0.005, \beta = 0.97, Re =1500$ and $k = 0.4{\rm \pi}$: (a) least-stable wall mode $c = 0.561844626 - 0.233038878\mathrm {i}$ from the A branch of the eigenspectrum; and (b) for the least-stable centre mode $c = 0.935239154-0.064928230\mathrm {i}$ from the P branch of the eigenspectrum. The location of the critical layer is shown using white dashed lines.

4.1. Collapse of neutral curves

The qualitatively similar character of the neutral curves at different $E$ in figure 18 is strongly suggestive of a collapse upon suitable rescaling of both $Re$ and $k$ with the elasticity number $E$. Figure 21 shows that such a collapse is indeed possible for sufficiently small $E$, when $Re$ is rescaled as $Re E^{3/2}$ and $k$ as $k E^{1/2}$. These scalings are found to be valid for fixed $\beta$, although the nature of the collapsed curve does depend on $\beta$ (as evident from figures 21a and 21c). Similarly, as shown in figures 21(b) and 21(d), the curves for the rescaled phase speed $(1-c_r)/E$, plotted as a function of $kE^{1/2}$, again exhibit a collapse, implying that $(1-c_r)$ is $O(E)$ for $E \ll 1$.

Figure 21.

Collapse of neutral curves for different $E$ in the $Re$–$k$ plane (a,c) and collapse of the corresponding phase speeds (b,d) for two different $\beta$: (a) rescaled neutral curves for $\beta = 0.65$; (b) $(1-c_r)$ collapse for $\beta =0.65$; (c) rescaled neutral curves for $\beta = 0.9$ and (d) $(1-c_r)$ collapse for $\beta =0.9$.

While the collapse obtained above is for a fixed $\beta$ and for $E \ll 1$, a further collapse is obtained in the dual limit $E (1-\beta ) \ll 1$, $(1-\beta ) \ll 1$, when the neutral curves are plotted in terms of $Re[(1-\beta )E]^{3/2}$ and $k [(1-\beta )E]^{1/2}$ as shown in figure 22, implying that the threshold $Re$ and $k$ scale as $Re \propto [(1-\beta )E]^{-3/2}$ and $k \propto [(1-\beta )E]^{-1/2}$, respectively, in this limit. The rescaled neutral curves in figure 22 begin to collapse onto a single one only for $\beta > 0.9$, the collapse being perfect for the lower branch, but less so for the upper ones. Thus, the role of the solvent viscosity appears to be ‘universal’ only as far as the lower branch is concerned. Importantly, however, since the critical $Re$ occurs on the lower branches of the neutral curves, the transition to the elasto-inertial turbulent state is governed by the combination $E(1-\beta )$ for $E (1-\beta ) \ll 1$, $(1-\beta ) \ll 1$. It is worth noting that the nearly vertical nature of the upper branch implies that the instability appears to exist in the limit of $Re \rightarrow \infty$, with $E$ fixed. An axisymmetric version of the ‘elastic Rayleigh’ equation (the elastic analogue of the classical Rayleigh equation; see Rallison & Hinch Reference Rallison and Hinch1995; Subramanian et al. Reference Subramanian, Reddy and Roy2020), which also has $E(1-\beta )$ as the governing parameter, is known to govern the linearized dynamics of perturbations in this limit, and involves a balance of inertial and elastic forces in the fluid. There is, however, no instability associated with the elastic Rayleigh equation for plane- (Kaffel & Renardy Reference Kaffel and Renardy2010) and pipe-Poiseuille (Chaudhary, Shankar & Subramanian Reference Chaudhary, Shankar and Subramanian2020) flows, and the lack of collapse of the (near-vertical) upper branches, and the implied instability for $Re \rightarrow \infty$, in figure 22, betrays therefore the singular nature of the inviscid elastic limit, with viscous effects playing a likely role even as $Re \rightarrow \infty$.

Figure 22.

Neutral curves for different $\beta$ and $E$ plotted in terms of $Re[(1-\beta )E]^{3/2}$ versus $k[(1-\beta )E]^{1/2}$: the rescaled neutral curves collapse for $\beta \rightarrow 1$.

4.2. Critical parameters and scalings

Figures 23(a), 23(c) and 23(d) show the variation of critical parameters $Re_c$, $k_c$ and $c_{r,c}$ with $E (1-\beta )$ for different $\beta$. Irrespective of $\beta$, the critical parameters conform to scaling laws for small $E(1-\beta )$; thus, $Re_c \propto (E(1-\beta ))^{-3/2}$, $k_c \sim (E(1-\beta ))^{-1/2}$ and $(1-c_{r,c}) \sim (1-\beta ) E$. Further, the curves for $\beta > 0.9$ collapse onto a universal curve in this limit (as was expected from the findings of the previous section in the dual limit $E(1-\beta ) \ll 1$, $(1-\beta ) \ll 1$. The collapse breaks down for $E(1-\beta ) > 0.05$, with the breakdown occurring at the point where the original neutral curves in the $Re$–$k$ plane start shifting upwards (after becoming two-lobed), with the lower lobe shrinking in size with increasing $E$ (figure 18). As $E (1-\beta )$ is increased, $Re_c$ reaches a minimum value and beyond a threshold value of $E (1-\beta )$, it increases rather sharply indicating the flow to be stable beyond this threshold. However, this threshold shifts to higher $E(1-\beta )$ as $\beta \rightarrow 1$, and $Re_c$ therefore continues to decrease for $\beta \rightarrow 1$, with the lowest $Re_c$ found being as small as $63$ (albeit for $E \sim 10$). The latter suggests that pipe flow of strongly elastic dilute polymer solutions can become unstable at an $Re$ much lower than that for their Newtonian counterparts. Figure 23(b) shows that $Re_c$ in figure 23(a) decreases approximately in a linear manner with $\beta$, although there appears to be an eventual deviation from linearity for the highest $\beta = 0.99$ analysed in this study. The aforementioned deviation from linearity suggests the approach of $Re_c$ to a finite lower bound regardless of $E$ or $\beta$, and that this lower bound is attained with $E(1-\beta )$ being finite. However, note that the corresponding $E$ diverges as $1/(1-\beta )$ for $\beta \rightarrow 1$, implying that the flow only becomes unstable for a very high $W$ in this limit.

Figure 23.

Variation of critical parameters with $E (1-\beta )$ for $\beta$ ranging from $0.4$ to $0.99$. (a) Plot of $Re_c$ versus $E(1-\beta )$, (b) the minima of $Re_c$ of (a) decreases approximately in a linear manner with $\beta$, but appears to approach a finite value as $\beta \rightarrow 1$, whereas the corresponding $E$ diverges as $\beta \rightarrow 1$, (c) $k_c$ versus $E(1-\beta )$ and (d) $(1-c_{r,c})$ versus $E(1-\beta )$. Here $Re_c$ and $k_c$ follow the scalings $Re_c\propto [E(1-\beta )]^{-3/2}$ and $k_c\propto [E(1-\beta )]^{-1/2}$, respectively, below a critical value of $E(1-\beta )$.

The scalings for the parameters ($Re$, $c_r$, $k$) with $E$ for $E\ll 1$, found previously, may also be justified using a scaling analysis for the boundary layer near the centreline, as briefly outlined in Garg et al. (Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018). In the limit $Re \gg 1$, $E\ll 1$, there is a ‘core’ region around the centreline with (dimensionless) extent $\delta \ll 1$, where inertial, elastic and viscous stresses are equally important. The scalings for $Re, k, \delta$ and $c_r$ in terms of $E$ can be derived by rescaling (2.8)–(2.14) in the region near the centreline as follows. The radial coordinate $r$ near the centreline can be expressed as $r=\delta \xi$, with $\xi \sim O(1)$. For $E\ll 1$, our numerical results show that $k_c$ becomes large for $E \ll 1$, and so we set $k \sim \delta ^{-1}$ in the analysis, as suggested by the continuity equation. The eigenvalue $c$ approaches unity in the said limit, and as $r \rightarrow 0$, $U\sim 1$, $(U-c) \sim \delta ^2$ and we therefore expand $c$ as $c=1+\delta ^2c^{(1)}$. The derivatives near the centreline get rescaled as ${\mathrm {d}}/{\mathrm {d}r}=({1}/{\delta })({\mathrm {d}}/{\mathrm {d}\xi }) \equiv \delta ^{-1} {D_1}$.

The base-flow profile becomes $U = 1-r^2 \equiv 1-\delta ^2\xi ^2$, and (2.8)–(2.14) take the following forms near the centreline:

(4.7)\begin{equation} \delta^{{-}1}(\mathrm{D_1}+\xi^{{-}1})\tilde{v}_r+\mathrm{i}k\tilde{v}_z=0, \end{equation}

(4.8)\begin{align} -\mathrm{i}k\delta^2(c^{(1)}+\xi^2)\tilde{v}_r &={-}\delta^{{-}1}\mathrm{D_1}\tilde{p}+\{\delta^{{-}1} (\mathrm{D_1}+\xi^{{-}1})\tilde{\tau}_{rr} \nonumber\\ &\quad +\mathrm{i}k\tilde{\tau}_{rz}-\delta^{{-}1}\xi^{{-}1}\tilde{\tau}_{\theta\theta}\} + \beta Re^{- 1}\delta^{{-}2}\mathrm{L_1}\tilde{v}_r, \end{align}

(4.9)\begin{align} -\mathrm{i}k\delta^2(c^{(1)}+\xi^2)\tilde{v}_z-2\delta\xi\tilde{v}_r &={-}\mathrm{i}k\tilde{p}+[\delta^ {{-}1}(\mathrm{D_1}+\xi^{{-}1})\tilde{\tau}_{rz} \nonumber\\ &\quad +\mathrm{i}k\tilde{\tau}_{zz}] +\beta Re^{{-}1}\delta^{{-}2}(\mathrm{L_1}+\xi^{{-}2})\tilde{v}_z, \end{align}

(4.10)\begin{equation} \{1- \mathrm{i}kW\delta^2(c^{(1)}+\xi^2)\} \tilde{\tau}_{rr} = 2(1-\beta)Re^{{-}1}(\delta^{- 1}\mathrm{D_1}-2W\mathrm{i}k\delta\xi)\tilde{v}_r, \end{equation}

(4.11)\begin{align} & \{1- \mathrm{i}kW\delta^2(c^{(1)}+\xi^2)\} \tilde{\tau}_{rz} +2W\delta\xi\tilde{\tau}_{rr} = (1- \beta)Re^{{-}1} [\{\mathrm{i}k \nonumber\\ &\quad -2W(1-\xi\mathrm{D_1}+4W\mathrm{i}k\delta^2\xi^2)\}\tilde{v}_r+(\delta^{{-}1}\mathrm{D_1}-2W \mathrm{i}k\delta\xi)\tilde{v}_z], \end{align}

(4.12)\begin{equation} \{1- \mathrm{i}kW\delta^2(c^{(1)}+\xi^2)\} \tilde{\tau}_{\theta\theta} = 2(1-\beta)Re^{{-}1}\delta^ {{-}1}\xi^{{-}1}\tilde{v}_r, \end{equation}

(4.13)\begin{align} & \{1- \mathrm{i}kW\delta^2(c^{(1)}+\xi^2)\} \tilde{\tau}_{zz} +4W\delta\xi\tilde{\tau}_{rz} = 2(1- \beta)Re^{{-}1} [{-}8W^2\delta\xi\tilde{v}_r \nonumber\\ &\quad +\{\mathrm{i}k-2W\delta\xi(\delta^{{-}1}\mathrm{D_1}-4W\mathrm{i}k\delta\xi)\}\tilde{v}_z], \end{align}

where $\mathrm {L_1}=(\mathrm {D_1}^2+\xi ^{-1}\mathrm {D_1}-\xi ^{-2}-k^2\delta ^{2})$. In the $r$-momentum equation, a balance of inertial stresses and solvent viscous stresses gives $\delta \sim Re^{-1/3}$. The left-hand side of the linearized constitutive equations reveal that in order for the elastic and viscous contributions to be of the same order, we require $W \sim \delta ^{-1}$ or, after using $\delta \sim Re^{-1/3}$, $Re \sim E^{-3/2}$. We thus obtain the following scaling relationships for $Re \gg 1$, $E \ll 1$, along the neutral curve:

(4.14a–d)\begin{equation} \delta \sim Re^{{-}1/3},\quad k\sim Re^{1/3},\quad Re\sim E^{{-}3/2},\quad \text{and}\quad (1-c)\sim Re^{{-}2/3}. \end{equation}

As shown in figure 24, the eigenfunctions for different $Re$ and $k$ along the lower branch of the neutral curve do exhibit a collapse when plotted as a function of the boundary layer coordinate $\xi$. It is also possible to obtain the following scalings for a fixed $k$ by rescaling the (2.8)–(2.14):

(4.15a–d)\begin{equation} \delta\sim Re^{{-}1/4},\quad Re\sim E^{{-}2},\quad (1-c)\sim Re^{{-}1/2},\quad \text{and}\quad \tilde{v}_z \sim Re^{1/4}\tilde{v}_r. \end{equation}

These scalings are illustrated by the collapse of the $\tilde {v}_r$ and $\tilde {v}_z$ eigenfunctions, corresponding to the unstable mode, as shown in figure 25 for a few selected pairs ($Re, W$) that are large enough to justify the limit $(Re, W)\rightarrow \infty$ such that $Re \sim E^{-2}$.

Figure 24.

Collapse of eigenfunctions at different $Re$ and $k$ along the lower branch of the neutral curve for $\beta = 0.4$: (a) $\tilde {v}_z$ and (b) $\tilde {v}_r$. The eigenvalues (with $c_i = 0$) for which the rescaled eigenfunctions are shown are $c = 0.996707, 0.995912, 0.995131, 0.994367$ and $0.993621$, respectively.

Figure 25.

Collapse of eigenfunctions at different $Re$ and $E$, but for a fixed $k = 1$ for (a) axial velocity $\tilde {v}_z$ and (b) scaled radial velocity $Re^{1/4} \tilde {v}_r$, on scaled radial axis in the limit $Re\rightarrow \infty$ and $W\rightarrow \infty$ for a fixed $W/Re^{1/2}$ and $\beta =0.5$, corresponding to the eigenvalues $c=0.994842+0.000112i,0.996321+0.000103i,0.996985+0.000093i$ and $0.997383+0.000085i$, respectively.

Figure 26(a) shows that the critical Reynolds number ($Re_c$) and critical wavenumbers ($k_c$) for different values of $E$ and $(1-\beta )E$ follow the scalings $Re_c\propto [(1-\beta )E]^{-3/2}$ and $k_c\propto [(1-\beta )E]^{-1/2}$ for $E(1-\beta ) \ll 1$. We had earlier reported (Garg et al. Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018) that $Re_c$ diverges weakly as $\beta ^{-1/4}$ for $\beta \rightarrow 0$, based on results extending down to a $\beta$ of $0.025$. However, new results for lower values of $\beta$ (down to $10^{-3}$) in figure 26(b) show that $Re_c$ does not diverge as $\beta ^{-1/4}$, but appears instead to diverge more weakly, or perhaps even asymptote to a constant. Thus far, we have not found any unstable mode in the UCM limit for the corresponding $Re$ and $E$. However, note that the structure of the centre mode changes qualitatively for the smallest $\beta$, characterized by the onset of small-scale oscillations (Chaudhary et al. Reference Chaudhary, Garg, Shankar and Subramanian2019), and our efforts thus far prevent us from discriminating between $Re_c$ approaching a constant with regards to a weak divergence in the said limit.

Figure 26.

(a) Rescaled critical parameters $Re_c E^{3/2}$ and $k_c E^{1/2}$ versus $(1-\beta )$ fall on straight lines of slopes $-3/2$ and $-1/2$, respectively, for smaller $(1-\beta )$ implying $Re_c\propto [(1-\beta )E]^{-3/2}$ and $k_c\propto [(1-\beta )E]^{-1/2}$ in the limit $\beta \rightarrow 1$. (b) Variation of critical parameters $Re_c$ and $k_c$ with $\beta$ for $E=0.01$.

4.3. Role of stress diffusion on the unstable centre mode

As discussed in the introduction, artificial stress diffusion is often used for regularization in DNS studies of viscoelastic flows (Sureshkumar & Beris Reference Sureshkumar and Beris1995; Sureshkumar et al. Reference Sureshkumar, Beris and Handler1997; Lopez et al. Reference Lopez, Choueiri and Hof2019). Recently, it has been shown that this additional diffusivity can qualitatively affect the stress dynamics (Gupta & Vincenzi Reference Gupta and Vincenzi2019), even to the extent of suppressing signatures associated with EIT (Sid et al. Reference Sid, Terrapon and Dubief2018). In this section, therefore, we briefly examine the effect of stress diffusion on the onset of the centre mode instability. The constitutive equation for the polymeric stress, (2.2), is now augmented with the stress diffusion term:

(4.16)\begin{align} W \left( \frac{\partial \boldsymbol{T}}{\partial t} + (\boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla}) \boldsymbol{T} -\boldsymbol{T}\boldsymbol{\cdot} ( \boldsymbol{\nabla} \boldsymbol{v} ) - ( \boldsymbol{\nabla} \boldsymbol{v} )^{\textrm{T}} \boldsymbol{\cdot} \boldsymbol{T} \right) + \boldsymbol{T} + \frac{D \lambda}{R^2} \nabla^2 \boldsymbol{T} = \frac{1-\beta}{Re} \{\boldsymbol{\nabla} \boldsymbol{v} +( \boldsymbol{\nabla} \boldsymbol{v})^{\textrm{T}}\},\end{align}

where $D$ is the stress diffusivity. El-Kareh & Leal (Reference El-Kareh and Leal1989) showed that the stress diffusion term owes its origin to the translational diffusion of the polymer molecules and estimated the diffusivity $D$ to be $O(10^{-12})\ \textrm {m}^2\ \textrm {s}^{-1}$. Using relaxation times $\lambda \sim 10^{-3}$ s reported in Chandra et al. (Reference Chandra, Shankar and Das2018) for polymer concentrations ${\sim }500$ ppm, and for tube diameters ${\sim }0.1$–$1$ mm, the dimensionless diffusivity $D \lambda /R^2 \sim 10^{-9}$–$10^{-7}$. It is useful to represent diffusive effects using a Schmidt number $Sc \equiv \nu /D = E /(D \lambda /R^2)$, with $Sc \rightarrow \infty$ representing the absence of diffusion. The linearized equations for viscoelastic pipe flow using (4.16) were solved using a spectral method. Additional boundary conditions are now required for the stress components. At the pipe wall ($r=1$), the stress equation is imposed without the diffusivity, whereas a regularity condition for the stress is imposed at the centreline (Beris & Dimitropoulos Reference Beris and Dimitropoulos1999; Lopez et al. Reference Lopez, Choueiri and Hof2019). A finite stress diffusivity regularizes the continuous spectrum modes and leads to an additional family of stable diffusive modes. The decay rate of this family increases with increasing $(D \lambda /R^2)$. However, the discrete modes existing in the absence of stress diffusion are only weakly perturbed for small values of the diffusivity $(D \lambda /R^2 \rightarrow 0)$.

Figure 27 shows that the $Re$ for onset of the centre mode instability increases with increasing ($D \lambda /R^2$), implying that the stress diffusivity has a stabilizing effect; for $D \lambda /R^2 \rightarrow 0$, the onset becomes independent of the diffusivity, approaching the values shown in figure 27. The threshold diffusivity for stabilization is seen to depend on $E$ and $\beta$. Importantly, the instability continues to exist for the experimentally relevant values of $D \lambda /R^2 \sim 10^{-9}$–$10^{-7}$, but would be suppressed at much larger values of $D \lambda /R^2 \sim 10^{-4}$–$10^{-2}$ (or, equivalently, $Sc < 1000$, for $E \sim 0.1$) used in earlier DNS studies (Sureshkumar & Beris Reference Sureshkumar and Beris1995; Sureshkumar et al. Reference Sureshkumar, Beris and Handler1997; Lopez et al. Reference Lopez, Choueiri and Hof2019). Thus, consistent with the results of Sid et al. (Reference Sid, Terrapon and Dubief2018), the results of figure 27 reinforce the importance of using simulation techniques, which avoid an artificially enhanced diffusivity, to access the axisymmetric structures associated with the centre-mode instability.

Figure 27.

The effect of stress diffusion (characterized by $D \lambda /R^2$) on the threshold $Re$ required for onset of instability at different $E, \beta$ and $k$.

4.4. Comparison with recent experimental studies and DNS

4.4.1. Comparison with experiments

We have replotted, in figure 28, the results of Samanta et al. (Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013) for the transition Reynolds number $Re_t$ as a function of $E(1-\beta )$, based on the reported viscosities and relaxation times of the different polymer solutions used in the experiments. The present theoretical results yield similar critical Reynolds numbers $Re_c$ only at much higher values of $E(1-\beta )$. Samanta et al. (Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013) estimated the relaxation time using the CaBER technique (Anna & McKinley Reference Anna and McKinley2001), in which the flow is extensional, and the polymer chains are highly stretched. However, the CaBER procedure is known to have some disadvantages in the estimation of relaxation time for polymers in low-viscosity solvents owing to inertia being neglected in the filament thinning dynamics. The CaBER relaxation time also exhibits a significant concentration dependence even below the nominal overlap concentration (Clasen et al. Reference Clasen, Plog, Kulicke, Owens, Macosko, Scriven, Verani and McKinley2006). The data for $Re_t$ from the experiments of Chandra et al. (Reference Chandra, Shankar and Das2018), also plotted in figure 28, shows good agreement with the theoretical $Re_c$; in that, both threshold $Re$ are of the same order of magnitude for comparable $E(1-\beta )$. Chandra et al. (Reference Chandra, Shankar and Das2018) used small-amplitude oscillatory strain experiments to infer the relaxation times; in contrast to CaBER, the polymer chains are not greatly perturbed about their equilibrium conformations. Although the threshold $Re$ from Chandra et al. (Reference Chandra, Shankar and Das2018) are comparable to theory, the latter predicts $Re_c \sim E^{-3/2}$ along the lower branch of the theoretical envelope, and $Re_t \sim E^{-1/2}$ in Chandra et al. (Reference Chandra, Shankar and Das2018). This difference in the scaling exponents could be due to shear thinning in the experiments, which can also significantly alter the parabolic nature of the base velocity profile. These effects are not accounted for in the Oldroyd-B model used in this study. The following scaling analysis examines the role of shear thinning on the scaling exponent characterizing the $Re_c$ versus $E$ behaviour for small $E$. We begin by noting the limiting behaviour of viscosity and relaxation time, for large $W$, for the more realistic FENE-P model, where shear thinning arises on account of the chains being finitely extensible (Bird, Dotson & Johnson Reference Bird, Dotson and Johnson1980):

(4.17a,b)\begin{equation} \eta_1 \sim \eta (\dot{\gamma} \lambda)^{{-}2/3},\quad \text{and}\quad \lambda_1 = \lambda (\dot{\gamma} \lambda)^{{-}4/3}, \end{equation}

where $\eta$ and $\lambda$ are viscosity and relaxation time at zero shear rate ($\dot {\gamma } = U/R$). The effective Reynolds number $Re_1$ and Weissenberg number $W_1$, evaluated using the shear-rate dependent viscosity and relaxation time, are given in terms of those involving the corresponding zero-shear-rate quantities, as

(4.18a,b)\begin{equation} Re_1 = \frac{\rho UR}{\eta_1} = E^{2/3}Re^{5/3},\quad \text{and}\quad W_1 = \frac{\lambda_1 U}{R}, \end{equation}

and the effective elasticity number $E_1$ becomes

(4.19)\begin{equation} E_1 = \frac{W_1}{Re_1}= \frac{\lambda_1 \eta_1}{\rho R^2} = E^{{-}1}Re^{{-}2}. \end{equation}

We now postulate that the scaling for $Re_c$ determined above for an Oldroyd-B fluid, is valid for a FENE-P fluid as well, but with $Re_1$ and $E_1$ replacing $Re$ and $E$ in order to account for the shear-rate dependence of viscosity and relaxation time. Using (4.18a,b)–(4.19) in the theoretical scaling $Re_{1,c} \propto E_{1}^{-3/2}$ gives $Re_c \propto E^{-5/8}$; the scaling exponent now being closer to that ($-1/2$) observed in experiments (Chandra et al. Reference Chandra, Shankar and Das2018, Reference Chandra, Shankar and Das2020). A similar argument has been used earlier to successfully account the effect of shear thinning on the onset of inertialess elastic instability in Taylor–Couette flow (Larson, Muller & Shaqfeh Reference Larson, Muller and Shaqfeh1994).

Figure 28.

Comparison of the present theoretical predictions with experimental results of Samanta et al. (Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013) and Chandra et al. (Reference Chandra, Shankar and Das2018).

4.4.2. Comparison with DNS of Lopez et al. (2019)

The recent DNS study by Lopez et al. (Reference Lopez, Choueiri and Hof2019) on viscoelastic pipe flow used the FENE-P model and showed that at a fixed $Re = 3500$, the flow fully relaminarizes as $W$ is increased, and at even larger $W$, the laminar state again becomes unstable, with the post-instability friction factor approaching the MDR asymptote. In figure 29, we show the neutral stability curve in the $W$–$k$ plane corresponding to the centre mode instability for $Re = 3500$ and $\beta = 0.9$ (parameters corresponding to the DNS of Lopez et al. Reference Lopez, Choueiri and Hof2019), according to which the flow is unstable in the range $176.9 <W<4783.6$. The closed loop in the $W$–$k$ plane, at a fixed $Re = 3500$, arises because the centre mode instability is absent both in the low- and high-$W$ limits, as can be inferred from the corresponding neutral curves in the $Re$–$k$ plane shown in figure 18(b). The range of $W$ corresponding to the $W$–$k$ loop where the linear instability of the centre mode is found is significantly higher than the range ($16<W<80$) over which EIT was observed in the simulations of Lopez et al. (Reference Lopez, Choueiri and Hof2019). However, the Oldroyd-B model used here does not account for shear thinning effects inherent in the FENE-P model used in their simulations. More work is thus needed to address these discrepancies between the theoretical predictions and experimental and/or DNS studies.

Figure 29.

Neutral stability curve in the $W$–$k$ plane at fixed $Re = 3500$ and $\beta = 0.9$; the region inside the loop is unstable.