2.1. Governing equations

We consider pressure-driven flow of an incompressible viscoelastic fluid in a channel with walls separated by a distance $2H$ (figure 1). The viscoelastic fluid is modelled using the Oldroyd-B constitutive equation (Larson Reference Larson1988), which is applicable to dilute polymer solutions wherein the polymer chains are assumed to be non-interacting, and each chain is modelled as an elastic dumbbell with beads connected by a linear infinitely extensible entropic spring. This model predicts a shear-rate independent viscosity and first normal stress difference in viscometric shearing flows. Many authors have used this model in the past to analyse instabilities in the flow of dilute polymer solutions in rectilinear (Sureshkumar & Beris Reference Sureshkumar and Beris1995b; Wilson, Renardy & Renardy Reference Wilson, Renardy and Renardy1999; Morozov & Saarloos Reference Morozov and Saarloos2007; Zhang et al. Reference Zhang, Lashgari, Zaki and Brandt2013; Garg et al. Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018), curvilinear (Shaqfeh Reference Shaqfeh1996) and cross-slot (Poole, Alves & Oliveira Reference Poole, Alves and Oliveira2007) geometries with considerable success. To render the governing equations dimensionless, we use the centreline maximum velocity of the laminar base state $U_{max}$ as the velocity scale, channel half-width $H$ as the length scale, $H/U_{max}$ as the time scale, $\eta U_{max}/H$ as the scale for the stresses and pressure and $\eta$ is the total solution viscosity. The dimensionless continuity and momentum equations are given by

(2.1)\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}=0, \end{gather}

(2.2)\begin{gather}Re \left(\frac{\partial \boldsymbol{u}}{\partial t}+(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{u} \right)={-}\boldsymbol{\nabla} {p}+ {\beta} \nabla^2 \boldsymbol{u} + \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\tau}. \end{gather}

Here, $Re = \rho U_{max} H/\eta$ is the Reynolds number based on the solution viscosity and $\beta = \eta _s/\eta$. The Oldroyd-B constitutive relation for the polymeric stress tensor, $\boldsymbol {\tau }$, in dimensionless form is given by

(2.3)\begin{equation} \boldsymbol{\tau}+{W} \left(\frac{\partial \boldsymbol{\tau}}{\partial t}+ (\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}) \boldsymbol{\tau}-(\boldsymbol{\nabla}\boldsymbol{u})^\intercal \cdot \boldsymbol{\tau}- \boldsymbol{\tau} \boldsymbol{\cdot}(\boldsymbol{\nabla} \boldsymbol{u})\right)=\left( 1-\beta \right)\left(\boldsymbol{\nabla} \boldsymbol{u}+\boldsymbol{\nabla}\boldsymbol{u}^\intercal \right). \end{equation}

Here, $W = \lambda U_{max}/H$ is the Weissenberg number and $\lambda$ is the microstructural relaxation time. The UCM model, which ignores the solvent contribution to the stress, is obtained from the Oldroyd-B model by setting $\beta = 0$, whereas the limit of a Newtonian fluid is obtained by setting $\beta = 1$.

Figure 1.

Schematic representation of the configuration consisting of pressure-driven flow in a channel of half-width $H$.

2.2. Base flow

The laminar base state whose stability is of interest here is the steady fully developed pressure-driven channel flow of an Oldroyd-B fluid, with the base-state velocity profile $U(z) = 1-z^2$ being identical to that of plane Poiseuille flow of a Newtonian fluid. However, unlike its Newtonian counterpart, the Oldroyd-B fluid exhibits a nonzero first normal stress difference ($T_{xx} - T_{zz}) = 8 (1-\beta ) W z^2$. Here, and in what follows, the velocity and stress fields corresponding to the base flow are denoted by uppercase letters.

2.3. Linearized governing equations

A temporal linear stability analysis of the aforementioned base flow is carried out by imposing infinitesimal perturbations (denoted by primes) to the base flow: $\boldsymbol {u}=\boldsymbol {U}+\boldsymbol {u'},\ p=P+p',\ \boldsymbol {\tau }=\boldsymbol {T}+\boldsymbol {\tau '}$. As Squire's theorem is valid for plane Poiseuille flow of an Oldroyd-B fluid (Bistagnino et al. Reference Bistagnino, Boffetta, Celani, Mazzino, Puliafito and Vergassola2007), we restrict our analysis to two-dimensional perturbations, which are considered as elementary Fourier modes of the form $f'(x,z,t)=\tilde {f}(z)\exp [\textrm {i}k(x-ct)]$, where $f'$ is the relevant disturbance field, $\tilde {f}(z)$ is the eigenfunction, $k$ is the dimensionless wavenumber and the eigenvalue $c =c_r+{\rm i} c_i$ is the complex wavespeed of perturbations. If $c_i > 0$, the perturbations grow exponentially with time leading to an instability. Substituting the Fourier mode representation for perturbations in the linearized governing equations yields, with $\textrm {d}_z = \textrm {d}/\textrm {d} z$, and primes on the base velocity profile $U(z)$ denoting derivatives with respect to $z$:

(2.4)\begin{gather} \textrm{d}_z \tilde {v}(z)+{\rm i} k\tilde {u}(z)=0,\end{gather}

(2.5)\begin{gather} Re [\textrm{i} k (U-c)\tilde{u}(z)+\tilde{v}(z)U']={-}\textrm{i}k\tilde{p}(z)+\beta(\textrm{d}_z^2-k^2)\tilde{u}(z) \nonumber\\ \hspace{9.5pc} +\textrm{i}k\tilde{\tau}_{xx}(z)+\textrm{d}_z\tilde{\tau}_{xz}(z), \end{gather}

(2.6)\begin{gather} Re \, \textrm{i} k (U-c)\tilde{v}(z)={-}\textrm{d}_z\tilde{p}(z)+\beta (\textrm{d}_z^2-k^2) \tilde{v}(z) \nonumber\\ \hspace{5.5pc} + \textrm{i}k\tilde{\tau}_{xz}(z)+\textrm{d}_z\tilde{\tau}_{zz}(z), \end{gather}

(2.7)\begin{gather} {[}1+ \textrm{i} k W (U-c)]\tilde{\tau}_{xx}(z) = (1-\beta) [2 \textrm{i} k \tilde{u}(z) + 4 {\rm i} k W^2 (U')^2 \tilde{u}(z) +2 W U' \textrm{d}_z\tilde{u}(z) \nonumber\\\hspace{3.5pc} -4W^2 U' U'' \tilde{v}(z)]+2W U' \tilde{\tau}_{xz}(z), \end{gather}

(2.8)\begin{gather} {[}1+\textrm{i} k W (U-c)]\tilde{\tau}_{zz}(z) = 2(1-\beta)[ d_z \tilde{v}(z) + \textrm{i} k W U' \tilde{v}(z)], \end{gather}

(2.9)\begin{gather} {[}1+\textrm{i} k W (U-c)]\tilde{\tau}_{xz}(z) = (1-\beta)[\textrm{d}_z \tilde {u}(z) + \textrm{i} k \tilde{v}(z) +2 \textrm{i} k W^2 (U')^2\tilde{v}(z) \nonumber\\ \hspace{3.5pc}- W U'' \tilde{v}(z) ]+W U'\tilde{\tau}_{zz}(z). \end{gather}

2.4. Numerical procedure

In order to determine the complex eigenvalue ($c$), we use a spectral collocation method (Boyd Reference Boyd1999; Weideman & Reddy Reference Weideman and Reddy2000), where the dynamical variables (velocity, pressure and stress perturbations) are expanded as a finite sum of Chebyshev polynomials and substituted in the above linearized differential equations. In our spectral formulation, we discretize all of the six (2.4)–(2.9), and the resulting generalized eigenvalue problem is of the form

(2.10)\begin{equation} \boldsymbol{\mathsf{A}}\boldsymbol{x}=c\boldsymbol{\mathsf{B}}\boldsymbol{x},\end{equation}

where $\boldsymbol{\mathsf{A}}$ and $\boldsymbol{\mathsf{B}}$ are coefficient matrices, and $\boldsymbol {x} = (\tilde {u},\tilde {v},\tilde {p},\tilde {\tau }_{xx},\tilde {\tau }_{xz},\tilde {\tau }_{zz})^\intercal$ is the vector comprising of the coefficients of the spectral expansion at the collocation points. The size of the $\boldsymbol{\mathsf{A}}$ matrix is $6N \times 6N$, where $N$ is the number of Gauss–Lobatto collocation points. The generalized eigenvalue problem is solved using the ‘polyeig’ eigenvalue solver of Matlab. To filter out the spurious eigenvalues associated with the spectral method, we run our spectral code for two different values of $N$, say, $400$ and $500$, and eliminate those eigenvalues that do not satisfy a prescribed tolerance criterion. In the following discussion, we usually use $N$ between $400$ and $600$ for $k,E < 1$. However, for the highest $E$, we use $N = 900$ to obtain convergence of the unstable mode. In addition, a numerical shooting procedure (Ho & Denn Reference Ho and Denn1977; Lee & Finlayson Reference Lee and Finlayson1986b; Schmid & Henningson Reference Schmid and Henningson2001) is used for further validation by providing the results from the spectral method as initial guesses. The numerical shooting procedure involves an adaptive Runge–Kutta integrator coupled with a Newton–Raphson iterative scheme to solve for the eigenvalues. Only physically genuine modes from the spectral method converge with the shooting code. To benchmark the implementation of our numerical methodology, we compare (table 1) results from our procedure with those of Sureshkumar & Beris (Reference Sureshkumar and Beris1995b) for both UCM and Oldroyd-B fluids. The unstable eigenvalues are in good agreement for $N=129$ and $N=257$. In addition, we have benchmarked our results with those of Chaudhary et al. (Reference Chaudhary, Garg, Shankar and Subramanian2019) for the UCM case.

Table 1.

Validation of UCM ($\beta =0$) and Oldroyd-B $(\beta =0.5)$ results with those of Sureshkumar & Beris (Reference Sureshkumar and Beris1995b) for viscoelastic channel flow.

3.1. The Newtonian and Oldroyd-B spectra

We first discuss the key differences between the Oldroyd-B and Newtonian eigenspectra. Note that the Oldroyd-B eigenspectrum reduces to the Newtonian eigenspectrum when either $E = 0$ (for any $\beta$) or $\beta = 1$ (for any $E$). As mentioned in § 1, the Newtonian eigenspectrum for plane Poiseuille flow (see figure 2a), at sufficiently high $Re$, has a characteristic ‘Y-shaped’ structure. For $Re > 5772$, a wall mode belonging to the ‘A’ branch becomes unstable (Schmid & Henningson Reference Schmid and Henningson2001), this being the TS instability. The eigenspectrum at $Re = 800$, $E = 0.1$, $\beta = 0.8$ and $k = 1.5$ (figure 2b) shows that in addition to the elastic modification of the discrete modes of the Newtonian spectrum, the spectrum for the Oldroyd-B fluid has a pair of CS ‘balloons’ (Graham Reference Graham1998; Wilson et al. Reference Wilson, Renardy and Renardy1999; Chaudhary et al. Reference Chaudhary, Garg, Shankar and Subramanian2019). The exact locations of the two continuous spectra are obtained in the following manner (Wilson et al. Reference Wilson, Renardy and Renardy1999; Chokshi & Kumaran Reference Chokshi and Kumaran2009). The linearized governing equations and constitutive relations (2.4)–(2.9) can be recast into a single fourth-order differential equation for $\tilde {v}_z$, wherein the coefficient of the highest-order derivative vanishes when $[1+\textrm {i} k W (U-c)] = 0$ and $[1+\textrm {i} \beta k W (U-c)] = 0$. This yields $c_i = -1/(kW)$ and $c_i = -1/(\beta k W)$, with $c_r = U(z)$, for the two CS; the latter condition implies $c_r \in [0,1]$. The CS with $c_i = -1/(kW)$ is present even in the absence of solvent (i.e. the UCM limit), and henceforth will be referred to as ‘CS1’. The second continuous spectrum (abbreviated as CS2), characterized by modes with $c_i = -1/(\beta k W)$, is present only for non-zero $\beta$. Theoretically, both the CS are ‘lines’ in the $c_r$–$c_i$ plane with the aforementioned $c_i$. As the eigenfunctions corresponding to the CS eigenvalues are singular, these are resolved only approximately by the finite number of collocation points used in the spectral method. Thus, both the CS appear as balloons whose spread only decreases slowly with increasing $N$. In addition to the elastically modified Newtonian discrete modes and the CS balloons, new discrete modes (absent in the Newtonian spectrum) also appear, of which one of the centre modes is unstable at $E = 0.1$ (see the inset of figure 2b); all other discrete modes, including the continuation of the TS (wall) mode, remain stable for $Re = 800$. An analogous centre-mode instability for viscoelastic pipe flow (over a similar range of parameters) was first reported by Garg et al. (Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018), and has since been examined in more detail by Chaudhary et al. (Reference Chaudhary, Garg, Subramanian and Shankar2021). The presence of a centre-mode instability in both channel and pipe flows of an Oldroyd-B fluid is in direct contrast to the Newtonian scenario, where pipe flow is stable at any $Re$.

Figure 2.

Eigenspectra for plane Poiseuille flow of (a) Newtonian ($E= 0$) and (b) Oldroyd-B ($E = 0.1$) fluids at $Re = 800,\ k = 1.5$ and $\beta = 0.8$. The A, P, and S branches of the Newtonian spectrum are indicated in (a). The inset in (b) zooms over the region near the unstable eigenvalue.

3.2. Evolution of the unstable elasto-inertial centre mode

In this section, we discuss the emergence and trajectory of the elasto-inertial centre mode (henceforth labelled ECM-1) that turns unstable for large enough $E$, and a few of the least-stable discrete stable modes, by following two different paths in parameter space, both starting from the Newtonian limit: (i) increasing $E$ (from zero) at fixed $\beta$ and (ii) decreasing $\beta$ (from unity) at fixed $E$; we also examine the non-trivial effect of changing $\beta$ on the spectrum in the complementary UCM limit ($\beta \rightarrow 0$).

3.2.1. Effect of varying $E$ at fixed $\beta$

The focus of the ensuing discussion is on the aforementioned centre mode; a detailed depiction of the overall features of the elasto-inertial spectrum, as a function of both $E$ and $\beta$, is provided in appendix A. In general, when $E$ is increased from zero at fixed $\beta$, $Re$ and $k$, the two CS moves up towards the $c_i = 0$ line (the $c_r$ axis), and the continuation of the centre modes originally present in the Newtonian P-branch (termed the NCM-$i$, $i$ being the mode number), merge with CS1 (located at $c_i = -1/(k E \,Re)$); new discrete modes also emerge from below the CS1 in this process (see appendix A.1). We further demonstrate (see appendix A.2) that the distinct class of discrete modes, comprising the shear waves in the UCM spectrum ($\beta = 0$), referred to here as the ‘HFGL’ (high-frequency Gorodtsov–Leonov) class of modes (Gorodtsov & Leonov Reference Gorodtsov and Leonov1967), are strongly stabilized at finite $\beta$. Hence, the continuation of the HFGL modes are not relevant in determining the stability of plane-Poiseuille flow in the dilute and semi-dilute regimes of relevance to experimental studies. At higher $E$, all the Newtonian centre modes have merged with the CS, and the CS-modes are, therefore, the least stable. Crucially, beyond a threshold $E$, a new ‘elasto-inertial’ centre mode emerges above the CS, and this mode becomes eventually unstable at sufficiently high $E$ (see figure 21 of appendix A.1).

Figures 3(a) and 3(b) present the eigenspectra for different $E$ varying over the interval $(0.4, 1.1)$ for $\beta = 0.96$, with figure 3(b) being plotted in terms of the scaled growth rate $k W c_i$, which ensures that the locations of the two CS are fixed as $E$ is changed (for a given $\beta$). Figure 3(a) tracks the paths taken (with increasing $E$) by the first few discrete modes, whereas the continuous line in figure 3(b) represents the trajectory of the unstable elasto-inertial centre mode (ECM-1) alone obtained from the shooting method (the superposed symbols correspond to results obtained using the spectral method). The new elasto-inertial centre mode, which emerges above the CS1 at $E \approx 0.4$, becomes unstable for $0.48 < E < 1.04$, but becomes stable again for $E>1.04$, with $|c_i|$ eventually scaling as $1/E$ for large $E$, quite similar to pipe flow (Chaudhary et al. Reference Chaudhary, Garg, Subramanian and Shankar2021). However, unlike pipe flow, the $c_r$ for the unstable mode exceed unity over some ranges of $E$.

Figure 3.

Eigenspectra for $Re=650,\ k=1,\ \beta =0.96$ at different $E$. (a) The full spectrum. (b) Enlarged view of (a) near the unstable eigenvalue expressed using the scaled growth rate $kWc_i$. The continuous (blue) line showing the trajectory of ECM-1 is obtained using shooting method, whereas symbols show results from the spectral method.

Figure 4 shows the velocity eigenfunctions ($\tilde {v}_x, \tilde {v}_z$) for different $E$, corresponding to some of the unstable centre modes shown in figure 3. The $\tilde {v}_x$ eigenfunctions are symmetric about the channel centreline (and are therefore shown only over one half of the channel), in marked contrast with the TS (wall) and NCM-1 modes, which are antisymmetric about the channel centreline; note that NCM-1 refers to the first (least-stable) mode originally on the P-branch of the Newtonian spectrum. The eigenfunctions have their peak amplitudes closer to the channel centreline, but are nevertheless spread across the entire channel for the moderate $Re$ considered here, similar to the centre-mode instability in pipe flow (Chaudhary et al. Reference Chaudhary, Garg, Subramanian and Shankar2021). This latter fact, that the unstable eigenfunctions for moderate $Re$ and $E$ are not localized near the channel centreline despite the phase speed being close to the maximum velocity of the base flow, needs to be emphasized since this contradicts earlier interpretations of our original report on the centre-mode instability (Shekar et al. Reference Shekar, McMullen, Wang, McKeon and Graham2019). The contours corresponding to the velocity ($\hat {v}_x(x,z), \hat {v}_z(x,z)$) and stream-wise component of the polymeric stress ($\hat {\tau }_{xx}(x,z)$) eigenfunctions of the unstable mode are shown in figure 5, and reinforce the relatively modest confinement (in the neighbourhood of the centreline) at $Re = 650$.

Figure 4.

Velocity eigenfunctions corresponding to unstable eigenvalues in figure 3 for $Re = 650$, $k = 1$, $\beta = 0.96$ and at different $E$: (a) axial velocity, $\tilde {v}_x$; (b) wall-normal velocity, $\tilde {v}_z$. The $\tilde {v}_x$ eigenfunctions are symmetric about the channel centre, and are shown over the half-domain $0 \leq z \leq 1$. The eigenvalues for which the eigenfunctions are shown here are $E = 0.7$, $c = 0.99856712 + 0.00204187\textrm {i}$; $E =0.9$, $c = 1.00087623 + 0.00130115\textrm {i}$; $E =1.0$, $c = 1.00121782 + 2.88573410 \times 10^{-4}\textrm {i}$.

Figure 5.

Contours of $v_x, v_z$ and $T_{xx}$ for unstable (symmetric) centre mode in the $x$–$z$ plane for $Re =650$, $\beta = 0.96$, $k=1$ and $E=0.7$. The unstable eigenvalue is $c = 0.99856712 + 0.00204187\textrm {i}$.

Figure 6(a) shows the variation of $c_i$ for the TS mode (TSM), NCM-1 and ECM-1 modes with $E$. In the near-Newtonian limit ($E \rightarrow 0$), TSM is the least-stable mode followed by NCM-1, whereas ECM-1 just emerges from the CS1 for $E \approx 0.01$. For $E \sim 0.01$, the decay rates of TSM and ECM-1 cross each other, and for all higher values of $E$, ECM-1 is the least-stable/unstable mode. For $E > 0.02$, both TSM and NCM-1 collapse into CS1 (figure 6a); we discuss this feature in more detail in appendix B where we compare the relative stability of these two modes for different values of $Re$, $k$ and $\beta$. Thus, the mode ECM-1 is the least-stable discrete mode for $E > 0.01$, and, in fact, is the only discrete mode that lies above the CS for $E > 0.02$; for $E > 0.1$, ECM-1 becomes unstable (inset of figure 6a). The corresponding behaviour of the phase speeds of the three modes is shown in figure 6(b), where the phase speeds for TSM and NCM-1 increase with $E$, before eventually merging into CS1, whereas the phase speed of ECM-1 remains almost constant (close to unity) over the range of $E$ spanned.

Figure 6.

Relative stability of the first three least stable eigenmodes, namely TS mode (TSM), elastically modified Newtonian centre mode (NCM-1) and the new elasto-inertial centre mode (ECM-1) at $Re = 800,\ k = 1.5$ and $\beta =0.8$. (a) Variation of $c_i$ with $E$ (inset shows the range of $E$ for which ECM-1 is unstable). (b) Phase speed $(c_r)$ corresponding to the modes shown in (a). In the Newtonian limit $(E\rightarrow 0)$, TSM is the least stable mode that governs the stability of the flow, whereas ECM-1 emerges from CS1 at $E \sim 0.01$. However, as $E$ increases, both TSM and NCM-1 disappear into CS1 leaving behind ECM-1 as the least-stable mode for $E > 0.02$, which eventually becomes unstable at $E \approx 0.1$.

Unlike elasto-inertial wall modes (Chaudhary et al. Reference Chaudhary, Garg, Shankar and Subramanian2019), the elasto-inertial centre mode remains stable in the UCM limit ($\beta = 0$) for channel flow, and remains so for $\beta$ below a finite threshold. Figure 7(a) explores the effect of varying $E$ on ECM-1 for $\beta = 0$ and $0.5$. In the UCM limit $(\beta = 0)$, as $E$ is increased from the Newtonian limit ($E \rightarrow 0$), $|c_i|$ eventually decreases to very small values (figure 7a). However, $c_i$ remains negative even for very large $E$ and, therefore, no centre-mode instability is found in the UCM limit for channel flow. An analogous behaviour is found for $\beta = 0.5$. At higher $\beta$, ECM-1 does become unstable as explained in the following paragraph.

Figure 7.

Effect of increasing $E$ on the elasto-inertial centre mode (ECM-1) for UCM and Oldroyd-B fluids. (a) Variation of $c_i$ for $Re = 800,\ k = 1.5$ and $\beta = 0$ and $0.5$. The centre mode remains stable for $\beta <0.5$ even at very large values of $E$, in stark contrast to pipe flow which remains unstable at much lower $\beta$. (b) Scaled growth rate of ECM-1 for $Re = 2500,\ k = 0.19,\ \beta = 0.58$ and $Re = 650,\ k = 1,\ \beta = 0.96$. Regardless of the value of $\beta$, ECM-1 in channel always emerges from CS1 ($c_i =-1/(kW)$) in the limit $E \rightarrow 0$.

While discussing the evolution of the elasto-inertial centre mode (ECM-1) in pipe flow at fixed $\beta$, and for different $E$, Chaudhary et al. (Reference Chaudhary, Garg, Subramanian and Shankar2021) identified two qualitatively different trajectories of ECM-1 depending upon $\beta$: For $\beta \geq 0.85$, ECM-1 collapses into CS1 in the limit $E\rightarrow 0$, and does not seem to have any connection with the Newtonian spectrum (and with the least-stable Newtonian centre mode NCM-1, in particular). However, for $\beta <0.85$, the unstable centre mode smoothly continues to the least-stable centre mode of the Newtonian eigenspectrum (labelled NCM-1 in this study). For channel flow, in marked contrast, the unstable elasto-inertial centre mode never smoothly continues to its Newtonian counterpart with decreasing $E$, within the parameter regimes explored. This is because ECM-1 and NCM-1 are modes with opposite symmetry (as will be seen later in figure 27 of appendix B), the tangential velocity eigenfunction for NCM-1 is antisymmetric about the channel centreline, while it is symmetric for ECM-1 as already seen in figure 4, with the former emerging out of CS1 at a (non-zero) threshold $E$, and the latter collapsing into CS1 at a smaller $E$, for any fixed $\beta$. For the pipe-flow case, the smooth connection made possible by the axisymmetry of both modes. Figure 7(b) reinforces this trend by showing the scaled growth rate of the least-stable elasto-inertial centre mode for two different $\beta$, both greater than $0.5$ (namely, 0.58 and 0.96). The range of $E$ for which elasto-inertial centre mode remains unstable increases with $\beta$. For both $\beta$, ECM-1 follows a trajectory similar to the one shown in figures 3 and 6(a). Thus, the elasto-inertial centre mode, whether unstable or otherwise, is not the continuation/elastic modification of least-stable Newtonian centre mode (NCM-1) for any $\beta$ (regardless of $Re$ or $E$).

3.2.2. Effect of varying $\beta$ at fixed $E$

We considered the effect of varying $\beta$ on the elasto-inertial spectrum both from the Newtonian ($\beta = 1$) and UCM ($\beta = 0$) limits. The following discussion focuses on the former limit; in appendix A.2, we examine the latter limit and show the qualitative changes in the elasto-inertial spectrum as $\beta$ is increased from zero.

In figure 8(a), we discuss the role of gradually decreasing $\beta$ from unity at a fixed $E$ on the first few discrete modes. This figure shows the trajectories of the two leading Newtonian centre modes (labelled NCM-1 and NCM-2); although these appear to emanate from the same point for $\beta = 1$, a closer examination reveals two distinct, but closely-spaced modes in the Newtonian spectrum. In addition to these Newtonian centre modes (NCM-1 and NCM-2), four new modes emerge from CS1. These modes (labelled ECM-1 to ECM-4 in figures 8(a) and 8b) arise because of the combined effect of polymer elasticity and solvent viscosity at non-zero $Re$, and hence do not have counterparts in the Newtonian spectrum. The unstable centre mode belongs to this class (ECM-1 in figure 8b). Except for ECM-1, however, all the other elasto-inertial centre modes remain stable over the entire range of $\beta$, from the Newtonian ($\beta = 1$) to the UCM ($\beta = 0$) limit, regardless of $Re$ and $E$. Figure 8(b) depicts the trajectory of ECM-1 with decreasing $\beta$, starting from its emergence out of CS1 at $\beta \approx 0.95$, using the scaled growth rate $W kc_i$ (the continuous red line represents results from the shooting method). Similar to the trend exhibited by ECM-1 for varying $E$ (at fixed $\beta$; see § 3.2.1), wherein the instability existed only over a finite range of $E$, the mode is unstable only over a range of $\beta$ at fixed $E$ in figure 8(b), and becomes stable again below a critical $\beta$. Thus, the trajectory of the unstable centre mode with varying $\beta$, at a fixed $E$, is similar in both pipe (see figure 12 of Chaudhary et al. Reference Chaudhary, Garg, Subramanian and Shankar2021) and channel (figure 8 of the present work) flows. However, in contrast to the pipe case, the unstable centre mode in channel flow persists even for $Re \sim O(1)$ in the limit $\beta \rightarrow 1$, albeit at high $E$. We discuss this in detail in § 4.2.

Figure 8.

Modification of Newtonian centre modes (NCM-1,-2) in the viscoelastic spectrum and the emergence of new elasto-inertial centre modes (ECM-1,-2,-3,-4) as $\beta$ is decreased from unity at $Re=800,\ k=1.5,\ E=0.1$: (a) NCM-1 and -2 and ECM-3, ECM-4; (b) ECM-1 and ECM-2. All the new elasto-inertial centre modes emerge from CS1 as $\beta$ is reduced from unity. In (a), the modes NCM-1 and 2 are distinct, but closely placed, in the Newtonian limit. The continuous lines represent results from the shooting method while symbols denote results from the spectral method. For clarity, only the filtered eigenspectrum is shown (with the CS balloons being absent). The theoretical location of CS1 is shown using dotted lines.

In figure 9(a), we exclusively focus on the centre mode ECM-1 to illustrate the importance of the solvent viscous contribution in rendering this mode unstable, by showing the variation of the scaled growth rate with $\beta$; figure 9(b) shows the variation of the corresponding phase speeds. At a fixed $Re, E$ and $k$, ECM-1 emerges from CS1 ($c_i=-1/kW$) as $\beta$ is decreased from the Newtonian limit ($\beta \rightarrow 1$). At a critical $\beta$ (close to unity for higher $E$) the elasto-inertial mode becomes unstable, and the range of $\beta$ in which ECM-1 is unstable increases with decrease in $E$. However, the mode becomes stable again as $\beta$ is decreased below a threshold. Crucially, for $\beta < 0.5$, we find that the centre mode always remains stable in channel flow, for any $E$ and $Re$. The absence of an instability for $\beta < 0.5$ reinforces our predictions from the spectral analysis (in the previous section) that for the centre-mode instability, solvent viscosity is essential along with inertia and elasticity, again in agreement with the pipe flow results of Garg et al. (Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018) and Chaudhary et al. (Reference Chaudhary, Garg, Subramanian and Shankar2021). However, for pipe flow, the centre mode becomes unstable even as $\beta \approx 10^{-3}$, for sufficiently high $Re$. Intriguingly, this feature is not present in viscoelastic channel flows.

Figure 9.

Effect of variation in $\beta$ on the scaled growth rate ($k W c_i$) of unstable elasto-inertial centre mode (ECM-1) for $Re=800,\ k=1.5$. (a) In the UCM limit ($\beta \rightarrow 0$), the centre mode remains stable even at very large values of $E$, illustrating the role of solvent viscosity in the centre-mode instability in channel flow. For a fixed $E = 0.1$, as $\beta$ is decreased from unity, ECM-1 emerges from CS1 and becomes unstable over a small range of $\beta$ ($0.8$–$0.7$). The unstable range of $\beta$ shifts towards $\beta \rightarrow 1$ for $E = 0.6$. (b) The corresponding variation of $c_r$ with $\beta$.

4.1. Scaled neutral curves

Figure 10 is strongly suggestive of a collapse of neutral curves, especially for the smaller $E$, on suitable rescaling. Figures 11(a) and 11(b) show a collapse of the different small-$E$ neutral loops onto a single master curve in the $Re E^{3/2}$–$k E^{1/2}$ plane, for the $\beta$ chosen in the aforementioned figures, implying that the threshold Reynolds number diverges as $E^{-3/2}$ as one approaches the Newtonian limit $E = 0$. In figures 11(c) and 11(d), the phase speeds along the neutral curve exhibit a similar collapse when plotted as $(1-c_r)/E$ versus $kE^{1/2}$, suggesting that $(1-c_r) \sim O(E)$ along the neutral curve. A similar collapse was also reported for pipe flow (Garg et al. Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018; Chaudhary et al. Reference Chaudhary, Garg, Subramanian and Shankar2021). An alternate route to the Newtonian limit, that of $\beta$ approaching unity for a fixed $E$, also appears to yield a collapse of the neutral curves when plotted in terms of $Re[E(1-\beta )]^{3/2}$ and $k[E(1-\beta )]^{1/2}$, in the limit $[E(1-\beta )]\ll 1$ (figure 12a). However, this collapse is not as perfect as that obtained previously for small $E$, even in the limit $\beta \rightarrow 1$. In particular, the upper branch of the $Re$–$k$ curves collapses very well for $\beta \approx 0.99$, but the collapse is not perfect in the lower branches and near the minimum of the neutral curves. Figure 12(b) shows the rescaled critical Reynolds number, $Re_c E^{3/2}$, and the corresponding rescaled critical wavenumber, $k_c E^{1/2}$, as a function of $(1-\beta )$. This plot suggests that $Re_c$ and $k_c$ begin to approach the scalings $Re_c \propto (E(1-\beta ))^{-3/2}$, $k_c \propto (E(1-\beta ))^{-1/2}$ only for $\beta \approx 0.99$.

Figure 12.

Collapse in the limit $(1-\beta ) \ll 1$ and $E$ fixed. In (a), neutral stability curves at different $E$ and $\beta$ plotted in terms of the scaled Reynolds number $Re[E(1-\beta )^{-3/2}]$ and wavenumber $k[E(1-\beta )]^{-1/2}$. For $\beta \rightarrow 1$, the rescaled neutral curves exhibit a data collapse. In (b), rescaled critical parameters at different $E$ and $\beta$ plotted as $Re_c E^{3/2}$, $k_c E^{1/2}$ vs. $(1-\beta )$ fall on lines of slopes $-3/2$ and $-1/2$ respectively, indicating again that, $Re_c \propto [E(1-\beta )]^{-3/2}$ and $k_c \propto [E(1-\beta )]^{-1/2}$.

4.2. Critical parameters and scalings

The critical parameters ($Re_c, k_c$ and $c_{rc}$) are plotted as a function of $E(1-\beta )$ in figure 13. The variation of $Re_c$ (figure 13a) is non-monotonic with $E(1-\beta )$, with $Re_c$ scaling as $(E(1-\beta ))^{-3/2}$ for $E(1-\beta ) \ll 1$, but showing a nearly vertical rise beyond a threshold $E$, denoted $E_{min}$, in a manner very similar to pipe flow (Garg et al. Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018; Chaudhary et al. Reference Chaudhary, Garg, Subramanian and Shankar2021). A similar non-monotonic behaviour of $Re_c$ with $E$ has been obtained for elasto-inertial wall mode instabilities in plane Poiseuille flow of Oldroyd-B (Sadanandan & Sureshkumar Reference Sadanandan and Sureshkumar2002; Brandi et al. Reference Brandi, Mendonça and Souza2019) and FENE-P (Zhang et al. Reference Zhang, Lashgari, Zaki and Brandt2013) fluids. However, because wall modes in channel flow are strongly stabilized by solvent viscous effects, the minima in $Re_c$–$E$ curves shift towards higher $Re_c$ with increase in $\beta$ for a fixed $E$ (see, for example, figure 1(a) of Sadanandan & Sureshkumar Reference Sadanandan and Sureshkumar2002). In stark contrast, for the unstable centre modes (figure 13a), the $Re_c$ shift towards lower values as $\beta$ approaches unity, thereby illustrating the contrasting roles played by solvent viscous effects on the centre- and wall-mode instabilities. Figure 13(b) further reinforces the effect of $\beta$ by showing the variation of the minimum $Re_c$ (obtained from figure 13a) and the corresponding $E_{min}$ with $(1-\beta )$. Unlike pipe flow, where the centre-mode instability ceases to exist below $Re_c \approx 60$, the instability in channel flow persists down to $Re_c \sim O(1)$ for $\beta \rightarrow 1$, albeit at very high $E$. Figure 13(c) shows that the critical wavenumber scales as $k_c \propto [E(1-\beta )]^{-1/2}$ for $E (1-\beta ) \ll 1$, whereas figure 13(d) shows that the critical phase speed scales as $( 1-c_r ) \propto [E(1-\beta )]$, both similar to pipe flow.

Figure 13.

Variation of critical parameters with $E (1-\beta )$: (a) the critical Reynolds number scales as $Re_c \propto [E(1-\beta )]^{-3/2}$ for $E(1-\beta ) \ll 1$; (b) the minimum $Re_c$ in (a) and the corresponding $E_{min}$; (c) critical wavenumber $k_c \propto [E(1-\beta )]^{-1/2}$; and (d) phase speed, $(1- c_r) \propto [E(1-\beta )]$. As shown in (b), the centre-mode instability persists in channel flow up to $Re \approx 5$ for very high $E \sim 10^4$ and for $\beta \approx 0.99$.

Similar to the collapse of the neutral curves for $E \ll 1$, a collapse is also exhibited by the eigenfunctions when plotted using a suitably rescaled wall-normal coordinate for $Re \gg 1$, $E \ll 1$. In this regard, there are two possible asymptotic regimes: one in which ($k, \beta$) are fixed and $Re$ and $E$ are varied so as to remain in the unstable region, and the other in which $\beta$ is fixed, and the eigenfunctions are tracked along different sets of critical parameters ($Re_c, k_c$) for different $E$. For the latter case, figure 14 shows that the tangential and normal velocity eigenfunctions are increasingly localized in the vicinity of the channel centreline, within a boundary layer of thickness of $O(Re^{-1/3 })$; the $Re$-dependence of this boundary layer thickness may be obtained using a scaling analysis, as outlined in Chaudhary et al. (Reference Chaudhary, Garg, Subramanian and Shankar2021). Instead, if one considers a fixed $k$, and the limit $Re,W \rightarrow \infty$, such that the ratio $W/Re^{1/2} \sim O(1)$ in order to be in the unstable region, the eigenfunctions become localized in a boundary layer of thickness of $\textit{O}(Re^{-1/4})$ in the vicinity of channel centreline.

Figure 14.

The collapse of stream-wise and wall-normal eigenfunctions corresponding to $Re_c$ and $k_c$ (at $\beta =0.8$ and different $E$) when plotted against the rescaled wall-normal coordinate scaled using the viscous layer thickness of $\textit{O}(Re^{-1/3})$: (a) axial and (b) wall-normal.

4.3. Effect of solvent viscosity on critical parameters

The centre-mode instability in pipe Poiseuille flow discussed in our earlier works (Garg et al. Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018; Chaudhary et al. Reference Chaudhary, Garg, Subramanian and Shankar2021), rather counter-intuitively, required the presence of solvent viscous effects, with the flow being stable in the UCM limit. Nevertheless, the pipe-flow instability does continue to exist for very low $\beta \sim 0.001$, with $Re_c$ exhibiting a weak divergence for $\beta \rightarrow 0$. In marked contrast, a finite solvent viscous threshold is required for the channel flow instability, with the instability ceasing to exist below $\beta \approx 0.5$ at $E =0.01$ (figure 15a). We have further verified that this is, in fact, the lowest $\beta$ for which the instability is present for any $E$. Figure 15(a) also shows a non-monotonic behaviour of $Re_c$ with $\beta$, at fixed $E \sim O(1)$, rather similar to the variation of $Re_c$ with $E$ (at fixed $\beta$). In the limit of $\beta \rightarrow 1$, $Re_c$ does diverge for channel flow, in a manner similar to that seen in pipe flow (see figure 5 of Garg et al. Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018). The divergence of $Re_c$ for $\beta \rightarrow 1$ appears, at first sight, to contradict the results shown in figure 13(b), where $Re_c$ decreases in the same limit. There is no inconsistency, however, because the parameters kept constant differ in the two cases. In figure 15(a), $E$ is fixed at $0.1$, whereas in figure 13(b), $E$ is allowed to vary, and increases to very high values for $\beta \rightarrow 1$. The eigenfunctions at the lowest $\beta$ for which the centre-mode instability is present are shown in figure 16. Interestingly, the eigenfunctions at $\beta = 0.6$ (and $Re = 800$, $k = 1.5$) are qualitatively similar to the eigenfunctions at a much higher $\beta = 0.96$ (and $Re = 650$, $k = 1$) shown in figure 4, suggesting that the shape of the centre-mode eigenfunctions is rather robust over the entire unstable range of $\beta$.

Figure 15.

Variation of (a) $Re_c$ and (b) critical wavenumber $k_c$ as a function of the viscosity ratio $\beta$ at fixed $E$. The minimum $\beta$ required to sustain the centre-mode instability in channel flow is ${\approx }0.5$.

Figure 16.

Normalized eigenfunctions for (a) stream-wise $\tilde {v}_x$ and (b) wall-normal $\tilde {v}_z$ perturbation velocities near the lowest value of $\beta$ for which the centre-mode instability exists in viscoelastic channel flow. Data shown for the eigenvalue $c = 0.99778 + 5.78112 \times 10^{-5} \textrm {i}$ at $Re =800,\ k=0.6,\ E=0.1,\ \beta =0.6$.

4.4. Role of diffusion on the centre-mode instability

In this section, we explore the role of stress diffusion on the centre-mode instability. The underlying microscopic origin of stress diffusion is the Brownian (translational) diffusion of the polymer molecules, with a diffusivity $D \sim 10^{-12}$ m$^2$ s$^{-1}$, and a corresponding Schmidt number $Sc = \nu /D \sim 10^6$, with $\nu$ being the kinematic viscosity of the polymer solution. To this end, the Oldroyd-B constitutive equation is now augmented with a stress diffusion term, whose importance, in dimensionless terms, is characterized by $D \lambda /H^2$ (Chaudhary et al. Reference Chaudhary, Garg, Subramanian and Shankar2021). Although many older (Sureshkumar & Beris Reference Sureshkumar and Beris1995a; Sureshkumar et al. Reference Sureshkumar, Beris and Handler1997) and a few recent (Lopez et al. Reference Lopez, Choueiri and Hof2019) DNS studies have incorporated an artificially large diffusion coefficient with $Sc \sim O(1)$, the work of Sid et al. (Reference Sid, Terrapon and Dubief2018) has demonstrated that the two-dimensional EIT structures are suppressed for $Sc < 9$. It therefore behooves us to examine whether stress diffusion has a similar effect on the centre-mode instability analysed in this study, especially because of our premise that the centre-mode instability is the mechanism underlying the onset of EIT. Based on the $D$ given previously, a typical relaxation time $\lambda \sim 10^{-3}$ s, and with channel half-width $H \sim 1$ mm, the dimensionless diffusivity $D \lambda /H^2 \sim 10^{-9}$. Note that, with the stress diffusion term included, boundary conditions need to be prescribed for the polymeric stress. Following earlier efforts (Sureshkumar & Beris Reference Sureshkumar and Beris1995a), these are obtained by using the constitutive equation without diffusion at the two boundaries. Figure 17 shows the threshold $Re$ for the centre-mode instability as a function of $D \lambda /H^2$, for fixed sets of $E$, $\beta$ and $k$. For $D\lambda /H^2 \rightarrow 0$, the threshold $Re$ for instability for the model with stress diffusion approaches that of the Oldroyd-B model without diffusion; importantly, $Re_c$ remains virtually unaltered for the aforementioned estimate of $D\lambda /H^2 \sim O(10^{-9})$. However, similar to pipe flow (Chaudhary et al. Reference Chaudhary, Garg, Subramanian and Shankar2021), $Re_c$ increases steeply for $D\lambda /H^2$ greater than a threshold that is a function of $E$ and $\beta$. For $(\beta ,E,k) \equiv (0.8, 0.16, 1)$, this threshold is $O(10^{-3})$, corresponding to $Sc = E/(D \lambda /H^2) \sim 100$ for $E = 0.1$. This stabilization of the linear centre-mode instability beyond a threshold stress diffusivity is broadly consistent with the disappearance of the span-wise structures in the fully nonlinear simulations of Sid et al. (Reference Sid, Terrapon and Dubief2018) discussed previously.

Figure 17.

The effect of stress diffusion coefficient $D \lambda /H^2$ on the threshold $Re$ required for centre-mode instability at different $E$, $\beta$ and $k$.

4.5. Comparison with experiments

We compare our theoretical predictions with the experiments of Srinivas & Kumaran (Reference Srinivas and Kumaran2017), who studied the flow of 30 and 50 ppm polyacrylamide (PAAm) solutions (molecular weight $5 \times 10^6$) through a rectangular microchannel with cross-sectional dimensions $160\ \mathrm {\mu }\textrm {m} \times 1.5$ mm, resulting in a cross-sectional aspect ratio of $9.375$. Srinivas & Kumaran (Reference Srinivas and Kumaran2017) characterized the transition using the increase in the standard deviation of velocity fluctuations, as inferred from PIV. On account of the large cross-sectional aspect ratio, the flow in the channel is expected to be quasi-two-dimensional, and have a parabolic (plane-Poiseuille) form with the characteristic scale in the gradient direction being the shorter of the two cross-sectional dimensions ($160\ \mathrm {\mu }$m). Therefore, for the purposes of comparison with our predictions, we identify the experimental channel half-width (see figure 1) to be $H = 80\ \mathrm {\mu }$m. We estimated the elasticity numbers ($E \equiv \lambda \nu /\rho H^2$, $\lambda$ being the longest relaxation time of polymer, whereas $\nu$ and $H$ are the kinematic viscosity of the solution and channel half-width), respectively, for these experiments using the Zimm relaxation time of the polymer solution. The longest relaxation time from the Zimm model is given by $\lambda _{Zimm} = \eta _s R_g^3/(k_B T)$ ((4.83) of Doi & Edwards Reference Doi and Edwards1986), where $R_g$ is the radius of gyration of a polymer molecule, $T$ is the absolute temperature, $k_B$ is the Boltzmann constant and $\eta _s$ is the solvent viscosity. It is important to note here that the Zimm relaxation time of $\lambda _{Zimm} = 30$ ms used in Srinivas & Kumaran (Reference Srinivas and Kumaran2017) is overestimated by a factor of around $20$ owing to the use of end-to-end distance in the expression given previously. However, if $R_g = 0.184 \times 10^{-6}$ m corresponding to a molecular weight of $5 \times 10^6$ g mol$^{-1}$ is used (Kulicke, Kniewske & Klein Reference Kulicke, Kniewske and Klein1982), we obtain $\lambda _{Zimm} = 1.4$ ms, for $\eta _s = 0.001$ Pa s. Using these estimates, we obtain $E = 0.22$ for $H = 80\ \mathrm {\mu }$m corresponding to the flow conditions of Srinivas & Kumaran (Reference Srinivas and Kumaran2017). The $Re_c$ from our stability analysis are in reasonable agreement with the threshold $Re_t$ inferred from experiments (table 2).

Table 2.

Comparison of present theoretical predictions for $Re_c$ with the experimentally inferred transition Reynolds number $Re_t$ of Srinivas & Kumaran (Reference Srinivas and Kumaran2017) for the flow of PAAm solutions in rectangular microchannels. Note that the $Re$ defined in Srinivas & Kumaran (Reference Srinivas and Kumaran2017) was based on the average velocity and channel-full width, and hence must be multiplied by $3/4$ to obtain the $Re$ used in the present work, which is based on the maximum base-flow velocity and channel half-width. Here, $C_p$ denotes the concentration of the polymer solutions used.

A point, made on more than one occasion in this paper, is that viscoelastic channel flow continues to be linearly unstable even at $Re \sim O(1)$, provided the elasticity number is sufficiently large. In figure 13(b), $Re_c$ dips down to about $5$ at an $E$ of $O(200)$ (with $\beta = 0.99$). In this regard, it is worth mentioning the recent experiments of Steinberg and co-workers (Varshney & Steinberg Reference Varshney and Steinberg2017, Reference Varshney and Steinberg2018a,Reference Varshney and Steinbergb), which demonstrate the feasibility of achieving very high $E$ with dilute polymer solutions. The experiments involve a channel flow set-up, although the focus is entirely different; the authors analyse elasticity-induced transitions in the free-shear flow set up between a pair of cylindrical obstacles embedded in the imposed pressure-driven flow. Importantly, the experiments access $W$ in excess of $10^3$ with $Re$ still being substantially smaller than unity. Note that the polymer concentration in the previous experiments is quite low ($c = 80$ ppm, with the overlap concentration $c^* \approx 200$ ppm), and shear thinning effects are therefore negligible. In contrast, there have been other reports of instabilities (Bodiguel et al. Reference Bodiguel, Beaumont, Machado, Martinie, Kellay and Colin2015; Poole Reference Poole2016; Picaut et al. Reference Picaut, Ronsin, Caroli and Baumberger2017; Chandra et al. Reference Chandra, Mangal, Das and Shankar2019) in channel/tube flows of highly shear-thinning concentrated solutions ($\beta < 0.2$), but these observations cannot be explained by the centre-mode instability which is absent for $\beta \leq 0.5$.