2.1. The FENE-P and FENE-CR models

In these models, the dilute polymer solution is again considered to be a suspension of non-interacting elastic dumbbells in a Newtonian solvent, as for the Oldroyd-B case, but with each dumbbell consisting of two beads connected by a nonlinear entropic spring. The polymer conformation is defined by the bead connector vector $\boldsymbol {Q}$ . In the FENE dumbbell, the (dimensional) force-extension relation characterising the entropic elasticity of the spring is given by $\boldsymbol {F}_c = ({H_{sp} \boldsymbol {Q}})/{1 -(|\boldsymbol {Q}|/Q_{max})^2}$ , where $H_{sp}$ is the (equilibrium) spring constant, with the force diverging in the limit $|\boldsymbol {Q}| \to Q_{max}$ . Here, $Q_{max}$ is the fully-stretched contour length of the polymer molecule. (Henceforth, the connector vector is rendered dimensionless using $\sqrt {(k_B T/H_{sp})}$ , and we continue to represent its non-dimensional analogue using $ \boldsymbol {Q}$ .) As mentioned in the Introduction, the nonlinear functional dependence on $\boldsymbol {Q}$ precludes the derivation of a closed-form constitutive relation for the stress. The configurational average involved in the stress (defined below), when applied to the underlying kinetic (Smoluchowski) equation, leads to a coupling to higher-order moments, and thence to an infinite moment hierarchy.

In order to obtain a closed-form constitutive equation, Peterlin (Bird et al. Reference Bird, Dotson and Johnson1980; Herrchen & Ottinger Reference Herrchen and Ottinger1997) replaced the $\boldsymbol {Q}$ -dependent nonlinear factor in the FENE spring force by its ‘pre-averaged’ equivalent, leading to the FENE-P model. This results in the following (dimensionless) relation between the polymeric stress $\boldsymbol {\tau } = \langle \boldsymbol {F}_c \boldsymbol {Q} \rangle$ and the conformation tensor $ \boldsymbol {C} = \langle \boldsymbol {Q} \boldsymbol {Q} \rangle$ :

(2.2) \begin{equation} \boldsymbol {\tau }=\frac {1-\beta }{W} (f(C_{kk}) \, \boldsymbol {C}-\boldsymbol {I}), \end{equation}

where $\langle \cdot \rangle$ represents the aforementioned configurational average. In (2.2), $ W = \lambda U_{max}/H$ is the Weissenberg number, $ C_{kk} = \text{tr}(\boldsymbol {C})$ , and $\boldsymbol {I}$ is the identity tensor representing the isotropic distribution of dumbbell conformations at equilibrium. The Peterlin closure function $ f(C_{kk})$ is defined as

(2.3) \begin{equation} f(C_{kk}) = \frac {1}{1-{C_{kk}}/{L^2}} , \end{equation}

where $ L^2 = Q_{max}^2/(k_B T/H_{sp})$ is the dimensionless parameter mentioned in the Introduction that characterises the maximum extensibility of polymer chains; $C_{kk} = Q_{eq}^2$ for $W = 0$ , with $Q_{eq}^2 = (3k_B T /H) L^2/(L^2 + 5)$ being the mean square extension at equilibrium. The Oldroyd-B model is obtained in the limit $L \to \infty$ ( $f \to 1$ ), when $\boldsymbol {F}_c$ becomes proportional to $\boldsymbol {Q}$ , implying a Hookean spring. It is worth noting (Alves, Oliveira & Pinho Reference Alves, Oliveira and Pinho2021; Davoodi et al. Reference Davoodi, Zografos, Oliveira and Poole2022) that there are (slightly) different variants of the FENE-P model prevalent in the literature, and the present work employs the version used by Zhang et al. (Reference Zhang, Lashgari, Zaki and Brandt2013) (see also Housiadas & Beris Reference Housiadas and Beris2003) in order to facilitate benchmarking of our numerical procedure. Importantly, the variants differ only by terms of $O(1/L^2)$ , and must therefore yield virtually identical results for $L \gtrsim O(10)$ . Herein, our focus is on the latter regime, implying that arguments with regard to shear thinning, in the Introduction, apply regardless of the particular variant used. The stability results detailed below should likewise be valid for the other FENE-P variants.

To describe the nonlinear rheology of Boger fluids (James Reference James2009), which exhibit negligible shear thinning of viscosity (owing to the presence of a highly viscous solvent), Chilcott & Rallison (Reference Chilcott and Rallison1988a ) proposed a model that leads to a shear-rate-independent viscosity, but retains the shear-rate dependence of the first normal stress coefficient. This is achieved by an (ad hoc) insertion of another factor $f$ multiplying the identity tensor in (2.2), leading to the FENE-CR model, with the stress given by

(2.4) \begin{equation} \boldsymbol {\tau } = (1-\beta ) \frac {f(C_{kk}) \, (\boldsymbol {C} - \boldsymbol {I})}{W}\, . \end{equation}

The evolution of the conformation tensor for both the aforementioned models can be expressed as

(2.5) \begin{eqnarray} (1-\beta ) \left (\frac {\partial \boldsymbol {C}}{\partial t}+(\mathbf {u}\cdot \nabla )\boldsymbol {C}-(\nabla \mathbf {u})^{\rm T} \cdot \boldsymbol {C}- \,\boldsymbol {C} \cdot (\nabla \mathbf {u})\right )= - \boldsymbol {\tau } , \end{eqnarray}

where $\boldsymbol {\tau }$ on the right-hand side takes the forms given in (2.2) and (2.4) for the FENE-P and FENE-CR models, respectively. As in the Oldroyd-B model, the terms within the parentheses on the left-hand side of (2.5) represent the upper-convected derivative operating on the conformation tensor $ \boldsymbol {C}$ , corresponding to an affine deformation of the underlying dumbbells in the limit $W \rightarrow \infty$ (i.e. in the absence of the restoring spring force).

2.2. Base-state velocity and stress distributions

In what follows, we use uppercase Roman alphabets to denote base-state quantities corresponding to the velocity field, stress field and components of the conformation tensor. For the Oldroyd-B fluid, on account of the shear-rate-independent viscosity, the steady unidirectional velocity profile for pressure-driven channel flow is the same as its Newtonian counterpart, the plane Poiseuille profile, with $U(z) = (1-z^2)$ . However, unlike the Newtonian case, the Oldroyd-B fluid exhibits a first normal stress difference $(T_{xx} - T_{zz}) = 8W^2z^2$ , and second normal stress difference $(T_{zz} - T_{yy}) = 0$ . In the coordinate-specific expressions above, and those that appear below, $x$ , $y$ and $z$ are, respectively, the flow (streamwise), vorticity (spanwise) and wall-normal directions, with perturbation velocity components denoted by $u$ , $w$ and $v$ .

Figure 1.

Base-flow velocity profiles for the FENE-P fluid for $\beta = 0.9$ at different $L$ values, for $W = 60$ and $150$ .

2.2.1. The FENE-P model

Unlike the Oldroyd-B model, shear thinning leads to the FENE-P model yielding a non-parabolic velocity profile (Cruz, Pinho & Oliveira Reference Cruz, Pinho and Oliveira2005; Zhang et al. Reference Zhang, Lashgari, Zaki and Brandt2013) given by

(2.6) \begin{eqnarray} U(z)=-\frac {1}{2\beta }\frac {{\rm d}p}{{\rm d}x}(1-z^2)-\frac {1-\beta }{\beta }\frac {3}{8B}( (F_+G_-)|_{-1}^{z} + (F_-G_+)|_{-1}^{z})\, , \end{eqnarray}

where

(2.7) \begin{eqnarray} F_+=(Bz + \sqrt {B^2z^2+A^3})^{1/3},\quad F_-=(Bz - \sqrt {B^2z^2+A^3})^{1/3}, \end{eqnarray}

(2.8) \begin{eqnarray} G_+=3Bz + \sqrt {B^2z^2+A^3},\quad G_-=3Bz - \sqrt {B^2z^2+A^3}\, , \end{eqnarray}

and

(2.9) \begin{equation} A=\dfrac {L^2}{6W^2}\left (1 + \dfrac {3}{L^2}+\dfrac {1-\beta }{\beta }\right ), \quad B=\dfrac {L^2}{4 W^2}\,\dfrac {{\rm d}p}{{\rm d}x}\,\dfrac {1}{\beta } . \end{equation}

Since the velocity scale ( $U_{max}$ ) used is the maximum speed for plane Poiseuille flow of an Oldroyd-B fluid, ${\rm d}p/{\rm d}x = -2$ in (2.6). Figure 1 shows the velocity profiles at different $L$ values for two different $W$ values. As $L$ is decreased, the maximum speed increases above unity owing to shear thinning.

In general, $U(z)$ in (2.6) is a function of $W$ , $L$ and $\beta$ . However, a notable feature in the above expressions is that, except for the $3/L^2$ term in the expression for $A$ in (2.9), $U(z)$ depends on $W$ and $L$ only through the combination $W/L$ . For $L \gt 10$ , the term $3/L^2$ is negligible compared to the other two terms in $A$ , and velocity profiles for different values of $W$ and $L$ , such that $W/L$ is constant, must therefore be virtually identical. Indeed, instead of using (2.3), if we use the FENE-P variant where $f = (L^2 - 3)/(L^2 - C_{kk})$ (see Alves et al. (Reference Alves, Oliveira and Pinho2021) for the prevalent FENE-P variants), then the $3/L^2$ term is absent, and the aforementioned collapse is exact for any $L$ (Tej, Mohanty & Shankar Reference Tej, Mohanty and Shankar2024). This collapse is analogous to that of the polymer contribution to viscosity $\eta _p$ and the first normal stress coefficient $\Psi _1$ , for plane Couette flow, demonstrated by Yamani & McKinley (Reference Yamani and McKinley2023), and is illustrated in figure 2 for different $W$ and $L$ . The collapse generally holds for steady unidirectional shear flows, including pipe flow (Tej et al. Reference Tej, Mohanty and Shankar2024). However, as pointed out by Yamani & McKinley (Reference Yamani and McKinley2023), it breaks down for purely extensional flows, due to the two- or three-dimensional nature of the velocity field in such flows.

Figure 2.

Master curves for base-flow velocity profiles of a FENE-P fluid at a fixed $W/L$ obtained for different pairs $(W,L)$ at (a) $\beta=0.5$ and (b) $\beta=0.9$ .

The base-state conformation tensor is given by

(2.10) \begin{equation} {\boldsymbol {C}} = \begin{bmatrix} \dfrac {1 + 2 (W U^{\prime})^2}{\bar {f}^3} &\quad 0 &\quad \dfrac {W U^{\prime}}{\bar {f}^2} \\ 0 &\quad \dfrac {1}{\bar {f}} &\quad 0 \\ \dfrac {W U^{\prime}}{\bar {f}^2} &\quad 0 &\quad \dfrac {1}{\bar {f}} \end{bmatrix} , \end{equation}

where $U^{\prime} = {\rm d}U/{\rm d}z$ , and $\bar {f}$ (with overbar denoting the base-state) satisfies the cubic relation $\bar {f}^3 -\bar {f}^2 -({2 W^2 (U^{\prime})^2})/{L^2} = 0$ . Of the three roots of this cubic, only one is physically admissible (i.e. real-valued). The base-state polymeric stress tensor, obtained from (2.2), is

(2.11) \begin{equation} {\boldsymbol {T}} = (1-\beta ) \begin{bmatrix} \dfrac {2 W U^{\prime2}}{\bar {f}^2} &\quad 0 &\quad \dfrac {U^{\prime}}{\bar {f}} \\ 0 &\quad 0 &\quad 0 \\ \dfrac {U^{\prime}}{\bar {f}} &\quad 0 &\quad 0 \end{bmatrix}. \end{equation}

Since the expression for $\bar {f}$ contains $U^{\prime}$ , its presence in $T_{xz}$ and $T_{xx}$ is the reason for the shear-rate dependence of both viscosity and the first normal stress coefficient.

2.2.2. FENE-CR model

The base-state velocity profile in the FENE-CR model is identical to the plane Poiseuille flow profile due to the shear-rate independent viscosity. The conformation and polymeric stress tensors in the base state are given by

(2.12) \begin{equation} \boldsymbol {C} = \begin{bmatrix} 1+2\left (\dfrac {WU^{\prime}}{\overline {f}} \right )^2 &\quad 0 &\quad \dfrac {WU^{\prime}}{\overline {f}} \\ 0 &\quad 1 &\quad 0 \\ \dfrac {WU^{\prime}}{\overline {f}} &\quad 0 &\quad 1 \end{bmatrix} ,\quad {\boldsymbol {T}}= (1-\beta ) \begin{bmatrix} 2\left (\dfrac {W}{\overline {f}} \right )U^{\prime} &\quad 0 &\quad U^{\prime} \\ 0 &\quad 0 &\quad 0 \\ U^{\prime} &\quad 0 &\quad 0 \end{bmatrix}\, . \end{equation}

Here, $\bar {f}$ satisfies a quadratic equation whose physically valid solution is

(2.13) \begin{equation} \overline {f} = \dfrac {L^2 + \sqrt {L^4 + 8(WU^{\prime})^2(L^2 - 3)}}{2(L^2 - 3)}\,. \end{equation}

While $\bar {f}$ above involves $U^{\prime}$ similar to the FENE-P model, the absence of the $1/\bar {f}$ factor in $T_{xz}$ in (2.12), leads to the viscosity being independent of the shear rate. Its presence in $T_{xx}$ , however, results in a shear-rate-dependent first normal stress coefficient.

2.3. Linearised equations for small perturbations

A temporal linear stability analysis is carried out using the two models above, by writing the different field variables as the sum of their base-state values (denoted by uppercase alphabets) and superposed infinitesimal perturbations (lowercase alphabets with primes): $\mathbf {u}=\mathbf {U}+\mathbf {u}^{\prime}$ , $p=P+p^{\prime}$ , $\boldsymbol {\tau }=\mathbf {T}+\boldsymbol {\tau }^{\prime}$ $\mathbf {c}=\mathbf {C} + \boldsymbol {c}^{\prime}$ . After substitution of these forms, the governing equations are linearised about the base state; similar to our earlier effort (Khalid et al. Reference Khalid, Chaudhary, Garg, Shankar and Subramanian2021a ,Reference Khalid, Shankar and Subramanian b ), we consider only two-dimensional perturbations, varying only in the streamwise and wall-normal directions. The use of two-dimensional perturbations is justified for Oldroyd-B fluids (within a modal stability analysis), owing to the existence of a Squire’s theorem, which states that corresponding to each unstable three-dimensional disturbance at given $Re$ , $W$ , there is a two-dimensional unstable mode at a lower $Re$ , $W$ , with a larger growth rate (Bistagnino et al. Reference Bistagnino, Boffetta, Celani, Mazzino, Puliafito and Vergassola2007). We are not aware of an extension of the Squire’s theorem for the FENE-P and FENE-CR models, but nevertheless assume two-dimensional disturbances in these cases. The two-dimensional perturbations are taken to be individual Fourier modes of the form $ f^{\prime}(x,z,t)=\tilde {f}(z)\exp[{\rm i}k(x-ct)]$ , where $f^{\prime}$ is any of the perturbation fields, $ \tilde {f}(z)$ is the eigenfunction, $ k$ is the real streamwise wavenumber, and $ c= c_r + {\rm i}c_i$ is the complex eigenvalue, with $c_r$ the wave speed, and $c_i$ the growth rate, and the flow being temporally unstable when $ c_i \gt 0$ . The linearised governing equations for both models are provided in Appendices A and B.

The linearised system is solved using a spectral collocation method, as explained in our earlier effort (Khalid et al. Reference Khalid, Chaudhary, Garg, Shankar and Subramanian2021a ). To benchmark our numerical implementation for the FENE-P model, we have compared results from our procedure with those of Zhang et al. (Reference Zhang, Lashgari, Zaki and Brandt2013); the eigenvalues are in good agreement (table 1). In addition, we have benchmarked the recent results of Buza et al. (Reference Buza, Page and Kerswell2022) (personal communication) both at finite $ Re$ and in the creeping-flow ( $ Re =0$ ) limit. Our FENE-CR results were benchmarked with those given in table 6 of Sureshkumar & Beris (Reference Sureshkumar and Beris1995); see table 2.

Table 1.

Validation of our FENE-P results with those of Zhang et al. (Reference Zhang, Lashgari, Zaki and Brandt2013) for viscoelastic channel flow at $ Re =3960,\ \beta =0.5,\ W=3.96,\ k=1.15,\ L = 60$ .

Table 2.

Validation of our FENE-CR results with those of Sureshkumar & Beris (Reference Sureshkumar and Beris1995): critical Reynolds number $Re_c$ and critical frequency $\nu _c$ at different $L$ .

3.1. Eigenspectra for the elastoinertial centre mode

Figure 3 shows the spectra obtained for both Oldroyd-B (figures 3 a,c) and FENE-P fluids (figures 3 b,d), the latter for $L = 200$ and $100$ , for $E = 0.1$ . In figures 3(a) and 3(c), the Oldroyd-B limit is obtained from the FENE-P fluid by taking $L = 10^4$ ; we show subsequently (in figure 4) that this value of $L$ is sufficient to capture the Oldroyd-B spectrum accurately, the spectrum being virtually invariant to a further increase in $L$ . As is well known, viscoelastic shearing flows, including the channel flow configuration examined here, are characterised by the presence of singular continuous spectra, whose (stable) eigenvalues lie on a line segment in the $c_r$ – $c_i$ plane. However, owing to the singular nature of the associated eigenfunctions, the eigenvalues belonging to the continuous spectra are not fully converged for any finite $N$ , with the theoretical line segment only being resolved as a balloon in the spectral method, the balloon extending to $c_i \gt 0$ if $N$ is not large enough. These eigenvalues, which are easily filtered out based on a convergence criterion, are therefore not considered in the ensuing discussion.

Figure 3.

Elastoinertial spectra for viscoelastic channel flow (a) in the Oldroyd-B limit ( $L = 10^4$ ), and (b) for finite extensibilities using the FENE-P model ( $L = 100,\ U_{max,F} = 1.09451$ and $L = 200,\ U_{max,F} = 1.05297$ ). Data for $Re = 800$ , $E = 0.1$ , $\beta = 0.8$ , $k = 1.5$ . (c,d) Zoomed-in regions of (a,b), respectively, near the elastoinertial centre mode (ECM-1).

Figure 4.

Dual role of finite extensibility on the centre mode illustrated using an enlarged view of the elastoinertial spectra, near the ECM obtained using the FENE-P model. (a) $Re = 150$ , $E = 0.4$ , $\beta = 0.9$ , $k = 1.4$ . The ECM is unstable in the Oldroyd-B limit and is monotonically stabilised by a decrease in $L$ . (b) $Re = 150$ , $E = 1$ , $\beta = 0.9$ , $k = 1$ . The ECM is stable in the Oldroyd-B limit, and is destabilised by an initial decrease in $L$ before being stabilised again for sufficiently small $L$ . The encircled eigenvalues are the ones for which eigenfunctions are shown in figure 5. (c) $Re = 150$ , $E = 1$ , $\beta = 0.9$ , $k = 0.4$ . The ECM is unstable in the Oldroyd-B limit. A decrease in $L$ leads to an initial stabilisation, a subsequent destabilisation, and an eventual stabilisation. (d) $c_i$ vs. $L$ for the adjacent panels. Summary of the results by showing the variation of $c_i$ with $L$ in (a,b,c).

Figures 3(a) and 3(c), respectively, show the overall spectra for the Oldroyd-B and FENE-P fluids, while figures 3(b) and 3(d) show the enlarged region near the unstable centre mode. Figure 3(c) shows an unstable centre mode with phase speed close to, but a little less than unity, for the Oldroyd-B fluid, the latter value being the base-state maximum. For the FENE-P fluid, figure 3(d) shows that the centre mode is unstable for $L = 200$ , while being stable for $L = 100$ . The phase speeds of these modes are slightly higher than unity, owing to the maximum velocity for the Oldroyd-B fluid being used as the velocity scale. In figures 3(b) and 3(d), $W = 80$ , with $W/L = 0.8$ and $0.4$ for $L = 100$ and $200$ , respectively. For such values of $W/L$ , figure 3(c) of Yamani & McKinley (Reference Yamani and McKinley2023) shows that the polymer contribution to the viscosity $\eta _p$ starts to decrease from its zero-shear value, and for a fixed pressure gradient, the (dimensionless) base-state maximum ( $U_{max,F}$ ) for the FENE-P fluid therefore increases beyond that for an Oldroyd-B fluid (i.e. $U_{max,F} \gt 1$ ; see caption to figure 3). Because the phase speed of the centre-mode instability is very close to $U_{max,F}$ , we find the $c_r$ values to be higher than unity for the FENE-P case. Another noteworthy feature in figures 3(c) and 3(d) is that while the centre mode is unstable in the Oldroyd-B limit, it is only marginally unstable for $L = 200$ , and is stable for $L = 100$ , suggesting a stabilising effect of finite extensibility at the chosen (small) $E$ .

We explore the above trend in detail, via a plot of $c_r$ versus $c_i$ for different $L$ , in figure 4(a). The figure shows, for a small $E = 0.4$ (slightly higher than $E = 0.1$ in figure 3), that decreasing $L$ from the Oldroyd-B limit down to $L = 40$ has a monotonically stabilising effect. This stabilisation can be rationalised by noting that for $W=60$ considered in figure 4(a), $W/L$ increases from $0$ (Oldroyd-B limit) to $1.5$ ( $L=40$ ). As the onset of shear thinning happens for $W/L \sim O(1)$ , one expects the shear thinning of $\Psi _1$ , and the concomitant decrease of the effective Weissenberg number $W_{eff} = \Psi _1 \dot {\gamma }^2/\tau _{xz}$ , to have a stabilising effect. The centre mode is seen to become stable (corresponding to $c_i$ crossing zero) for $L \leqslant 180$ , with $W/L \geqslant 1/3$ .

Figure 5.

Velocity and stress perturbation eigenfunctions corresponding to the ECM in figure 4, for $Re = 150$ , $E = 1.0$ , $\beta = 0.9$ , $k = 1.0$ , for different $L$ . All eigenfunctions are shown over the half-channel domain ( $0 \leqslant z \leqslant 1$ ) on account of symmetry. The corresponding eigenvalues are: $L = 40$ ( $U_{max} = 1.08012$ ), $c = 1.07575 - 0.00118415{\rm i}$ , $z_c = 0.066$ ; $L = 100$ ( $U_{max} = 1.06042$ ), $c = 1.06141 + 0.00174635{\rm i}$ , $z_c = 0.00315$ ; Oldroyd-B ( $U_{max} = 1$ ), $c = 1.00097 - 0.00224793{\rm i}$ . Here, $z_c$ denotes the critical-layer locations for the stable eigenmodes.

In sharp contrast, figure 4(b) considers a parameter regime where the centre mode is stable in the Oldroyd-B limit; the key difference relative to figure 4(a) is the increased value of $E\,(= 1)$ . The centre mode now becomes unstable with decreasing $L$ . This counterintuitive destabilising role of finite extensibility is in direct contradiction with the argument of weakened elastic effects due to shear thinning, used above to rationalise the stabilisation at smaller $E$ . This prediction is nevertheless consistent with the results of Buza et al. (Reference Buza, Page and Kerswell2022), who also demonstrated the instability of the FENE-P centre mode in parameter regimes where it is stable in the Oldroyd-B limit. The destabilising role of finite extensibility is present over a limited range of $L$ values ( $45 \lt L \lt 150$ ), and there is an eventual stabilisation for $L \approx 40$ and lower, corresponding to $W/L \gt 3.5$ . Thus as $L$ is decreased to sufficiently small values (at fixed $W$ ), the eventual significant increase in $W/L$ and the implied approach to the high-shear-rate asymptote renders the argument, based on stabilisation due to finite-extensibility-induced shear thinning, consistent. Crucially, the value of $L$ for the polyacrylamide used in the experiments of Choueiri et al. (Reference Choueiri, Lopez, Varshney, Sankar and Hof2021) falls within the above unstable range of $L$ values: using $M_w = 5 \times 10^6$ g mol $^{-1}$ , $M_m \approx 70$ g mol $^{-1}$ and $N_m \sim 5$ in $L = \sqrt {M_w/(M_m N_m)}$ yields $L \approx 120$ . Thus the predicted destabilisation due to finite extensibility is likely relevant to the observed onset of EIT in the said experiments.

In figure 4(c), corresponding to the same $E$ as in figure 4(b) but a smaller $k$ , rather remarkably, we see a behaviour with respect to $L$ that is doubly non-monotonic. That is, the centre mode is unstable in the Oldroyd-B limit, but stabilises initially with a decrease in $L$ . It again becomes unstable for $200 \lt L \lt 400$ , before eventually becoming stable for $L \lt 200$ . While we do not have a physical argument to explain the rather puzzling destabilising role played by finite extensibility, for the $O(1)$ values of $W/L$ in figures 4(b) and 4(c), we attempt to isolate the probable reasons by considering the FENE-CR fluid in § 4. Figure 4(d) summarises the contrasting behaviour of $c_i$ , as a function of $E$ , seen in figures 4(a–c); the inset shows the enlarged region where $c_i \gt 0$ (near $L \approx 200$ ) for $k = 0.4$ .

The absolute values of the eigenfunctions for the velocity ( $\tilde {v}_z$ , $\tilde {v}_x$ ) and stress ( $\tilde {\tau }_{xz}$ , $\tilde {\tau }_{xx}$ ) perturbations, for different $L$ values, are shown in figure 5; for reasons related to symmetry (or anti-symmetry) of the field variables, the eigenfunctions are shown only for one-half of the channel. The three values of $L$ correspond to the encircled eigenvalues in figure 4(b), and are such that the flow is stable in the Oldroyd-B limit, unstable for $L = 100$ , and stable again for $L = 40$ . Despite this qualitative difference in the nature of stability, the velocity eigenfunctions are qualitatively similar for both Oldroyd-B and FENE-P models. It is only in the stress eigenfunctions, the $\tilde {\tau }_{xx}$ eigenfunction in particular, that one begins to see qualitative differences, with the peak of the FENE-P eigenfunctions being more prominent and localised near the channel centre for the smaller $L$ . This localisation points to the development of a critical layer for the stable eigenmodes. Apart from the critical-layer signature, however, there is no clear distinguishing feature between the two FENE-P stress eigenfunctions that might explain the difference in stability.

Figure 6.

Elastoinertial spectra for $(W,L) \equiv (150,80)$ and $(300,160)$ such that $W/L$ is fixed: data for $Re = 150$ , $\beta = 0.9$ : (a) overall spectra, (b) zoomed-in spectra. Despite $W/L$ being the same, the two spectra are different. For each pair, the spectra are shown for two $N$ values in order to ensure convergence of the discrete eigenvalue.

It was argued in § 2.1 that for a fixed $\beta$ , the base-state velocity profiles are identical for a given $W/L$ . In figure 6, we show the spectra for two $(W,L)$ pairs, such that $W/L$ is fixed. Despite the base-state profiles being identical for the two cases, their eigenspectra are not, implying that stability is not governed by the combination $W/L$ , but by $W$ and $L$ independently; this is consistent with the form of linearised equations in Appendix A. Indeed, the mode that is unstable for $W = 150,\ L = 80$ becomes stable for $W = 300,\ L = 160$ . The lack of collapse of spectra with the same $W/L$ is because despite the base-state flow being a unidirectional shearing flow, the perturbed flow is not.

Figure 7.

Neutral stability curves in the $Re$ – $k$ plane for different $E$ , for $\beta = 0.8$ : (a) Oldroyd-B and (b) FENE-P ( $L = 70$ ). The role of finite extensibility is to suppress the instability at higher $Re$ values, leading to a closed unstable region (an ‘island’) in the $Re{-}k$ plane. The threshold (minimum) $Re$ is largely unaffected by a decrease in $L$ for the range of $E$ examined.

3.2. Neutral stability curves for the elastoinertial centre mode

The neutral stability curves in the $Re$ – $k$ plane, for fixed values of $E$ and $\beta$ , are shown in figure 7 for both Oldroyd-B and FENE-P models. For the Oldroyd-B case, the neutral curves take the form of unstable tongues in the $Re$ – $k$ plane (Khalid et al. Reference Khalid, Chaudhary, Garg, Shankar and Subramanian2021a ), with the unstable region extending all the way down to $k = 0$ , and with the upper and lower branches of each tongue conforming to the scaling $Re \propto k^{-1}$ for $k\ll 1$ . Within the FENE-P framework, as $k$ is decreased, the associated increase in $Re$ implies an increase in $W$ (since $E$ is fixed), thereby making shear thinning important. As already discussed in the previous subsection, the latter leads to a weakening of the elastic forces and an eventual stabilisation. Thus for any finite $L$ , one generally expects the unstable neutral tongues, present in the Oldroyd-B limit, to close out at sufficiently large $Re$ , leading to an unstable island instead. Further, the instability must disappear when the unstable island becomes vanishingly small, as will be seen to happen for sufficiently small $L$ . This feature has been recognised earlier in the context of hoop-stress-driven instabilities in curvilinear flows (Byars et al. Reference Byars, Öztekin, Brown and McKinley1994; McKinley et al. Reference McKinley, Öztekin, Byars and Brown1995; Castillo Sánchez et al. Reference Castillo Sánchez, Jovanović, Kumar, Morozov, Shankar, Subramanian and Wilson2022). For the centre-mode instability, the existence of unstable islands, as opposed to tongues, in the $Re$ – $k$ plane is confirmed in figure 7(b). This, of course, implies the existence of both upper and lower $W$ -thresholds for instability at a fixed $k$ . Further, for the chosen range of $E$ , the lower threshold $W$ required for instability is not significantly different from that for the Oldroyd-B model in figure 7(a), despite the relatively modest value of $L\,(= 70)$ used.

Figure 8.

A comparison of neutral curves in the $Re{-}k$ plane, corresponding to the centre-mode instability, obtained using the Oldroyd-B and FENE-P models, for (a) $\beta = 0.8$ , $E = 0.28$ , (b) $\beta = 0.8$ , $E = 0.42$ , and (c)  $\beta = 0.9$ , $E = 1.05$ ( $L = 100$ ). Finite extensibility has a stabilising effect in (a); in contrast, in (b) and (c), finite extensibility destabilises the centre mode, leading to instability even for $k$ values where the Oldroyd-B fluid remains stable.

In figure 8, we demonstrate the contrasting effects of decreasing $L$ , on the neutral tongues in the $Re$ – $k$ plane present in the Oldroyd-B limit, with a specific focus on the threshold $Re$ ( $Re_c$ , interpreted as the minimum of the neutral islands in the $Re$ – $k$ plane) required for the centre-mode instability. Figure 8(a) shows the role of decreasing $L$ for a relatively small $E$ ( $= 0.28$ ): both $L = 70$ and $L=80$ exhibit a closed unstable loop in the $Re$ – $k$ plane, for reasons mentioned above. As $L$ is decreased from infinity down to $70$ , $Re_c$ gradually increases, with a concomitant decrease in the size of the island; this trend is largely consistent with the stabilising role of shear thinning, for small $E$ , seen earlier in figure 4(a). However, as $E$ is increased to $0.42$ (figure 8 b), the counterintuitive destabilising role of finite extensibility becomes evident, with $Re_c$ for the FENE-P model becoming lower than that in the Oldroyd-B limit. Interestingly, the range of wavenumbers over which the flow is unstable also shifts towards higher values with decreasing $L$ . Figure 8(c) illustrates the stronger destabilisation due to finite extensibility at a higher $\beta = 0.9$ and $E = 1.05$ . Here, in contrast to figures 8(a) and 8(b), the unstable region, while still an island, is significantly larger in its spatial extent than the original Oldroyd-B tongue. The decrease in $Re_c$ is also substantial: from $Re_c \approx 300$ (for the Oldroyd-B fluid) to $Re_c \approx 50$ for $L=50$ .

Figure 9.

Variation of critical parameters ( $Re_c$ and $k_c$ ) with $E$ for Oldroyd-B and FENE-P models at $\beta = 0.8$ : (a) critical Reynolds number $Re_c$ ; (b) critical wavenumber $k_c$ .

Figures 9(a) and 10(a) show the variation of the critical Reynolds number $Re_c$ , and figures 9(b) and 10(b) the corresponding wavenumber $k_c$ , with $E$ , for a FENE-P fluid with $L = 100$ , and for two different values of $\beta$ . Results obtained using the Oldroyd-B model are shown for the purposes of comparison. For the lower $\beta \ (=0.8)$ , figure 9(a) shows that $Re_c \propto E^{-3/2}$ for $E \ll 1$ , for both Oldroyd-B and FENE-P fluids, with the only difference being the slightly higher $Re_c$ (at a given $E$ ) for the latter. This slight increase could again be attributed to the expected stabilising role of shear thinning, owing to the higher $W$ values along this branch. The $E^{-3/2}$ scaling above can be rationalised (Chaudhary et al. Reference Chaudhary, Garg, Subramanian and Shankar2021) using a boundary-layer-type scaling argument that posits that at high $Re$ , the disturbance is confined to a thin region near the channel centreline, with elastic, viscous and inertial effects being equally dominant. The $Re_c \propto E^{-3/2}$ scaling implies that $W_c \propto E^{-1/2}$ for $E \ll 1$ ; thus at sufficiently small $E$ values (i.e. $E \lt L^{-2}$ ), $W_c/L$ will be large enough for the centre mode to be stabilised, and one expects a departure from the $E^{-3/2}$ scaling. However, in figure 9(a), for the lowest $E$ examined for the FENE-P model ( $E \sim 0.05$ ), $W_c \sim 40$ and $W/L \sim 0.4$ , implying that the stabilising influence of finite extensibility is not yet significant. It is possible that for the chosen $L$ , at sufficiently small $E$ values (lower than those examined in figure 9 a), the centre mode is stabilised in the FENE-P model due to shear thinning of $\Psi _1$ .

For larger $E$ , in figure 9(a), $Re_c$ increases sharply, appearing to diverge at a finite $E\ (=0.4)$ for the Oldroyd-B case; the increase is more gradual for the FENE-P fluid. Note that the latter increase is on account of the unstable island shifting upward in the $Re$ – $k$ plane (data not shown), while simultaneously decreasing in size. Given the eventual small size of the shrinking island, it is difficult to determine whether the FENE-P increase is also indicative of a finite- $E$ divergence. Nevertheless, as seen in figure 9(a), the instability does appear to be absent for both models at sufficiently high $E$ . Figure 10(a) shows a similar trend for $\beta = 0.9$ . The divergence of the Oldroyd-B threshold is now postponed until $E\approx 1$ , while the increase in $Re_c$ for the FENE-P model is more gradual than for $\beta = 0.8$ . Figures 9(b) and 10(b) show that the threshold wavenumber mirrors the variation in $Re_c$ : $k_c$ plummets to zero at a finite $E$ for the Oldroyd-B case, while exhibiting a slower decrease for FENE-P, with this decrease being more gradual for $\beta = 0.9$ .

Figure 10.

Variation of critical parameters ( $Re_c$ and $k_c$ ) with $E$ for the Oldroyd-B and FENE-P models at $\beta = 0.9$ : (a) critical Reynolds number $Re_c$ ; (b) critical wavenumber $k_c$ .

Figure 11.

Variation of critical parameters ( $Re_c$ and $k_c$ ) with $E$ for the Oldroyd-B and FENE-P models at $\beta = 0.97$ : (a) critical Reynolds number $Re_c$ ; (b) critical wavenumber $k_c$ . The instability persists in the limit $E \rightarrow \infty$ , $Re \rightarrow 0$ .

For the Oldroyd-B fluid, the above trend continues until $\beta = 0.990552$ , after which there is a continuous decrease of $Re_c$ with $E$ , all the way up to $E \rightarrow \infty$ , implying a smooth crossover from the elastoinertial centre mode to a purely elastic centre-mode instability. Figure 11(a) shows that for a $\beta \ (= 0.97)$ lower than the critical value above, while the Oldroyd-B results still show the instability to be absent beyond a threshold $E$ , the FENE-P results show a continuation of the elastoinertial centre-mode instability to the elastic one instead. As seen in figure 11(b), $k_c$ now remains $O(1)$ even as $E \rightarrow \infty$ in the FENE-P case. In the next subsection, we proceed to examine in detail this purely elastic centre-mode instability, in FENE-P channel flow, for $\beta \geqslant 0.97$ .

3.3. The elastic centre mode: $Re = 0$ ( $E \rightarrow \infty$ )

Figure 12.

Elastic spectrum for pressure-driven channel flow of (a) FENE-P ( $L = 100$ , $U_{max} = 1.01254$ ) and (b) Oldroyd-B fluids, for $\beta = 0.98$ , $k = 2$ , $W = 200$ . (c,d) Zoomed-in regions of (a) and (b), respectively, near the centre mode (ECM).

In this subsection, we focus exclusively on the creeping flow limit and examine parameter regimes corresponding to the elastic centre-mode instability, within the framework of the FENE-P model. The elastic spectra for pressure-driven channel flow for $\beta = 0.98$ , for both FENE-P and Oldroyd-B fluids, are shown in figures 12(a–d); as for the elastoinertial case in figure 3, figures 12(a) and 12(b) show the overall spectrum, while figures 12(c) and 12(d) focus on the region close to the centre mode. For the chosen $\beta$ , the earlier effort by Khalid et al. (Reference Khalid, Shankar and Subramanian2021b ) has shown that the centre mode is stable for the Oldroyd-B fluid (the unstable interval being $0.990552 \lt \beta \lt 1$ ). In contrast, the elastic spectrum for the FENE-P fluid shows an unstable centre mode, as expected from the results shown in figure 11(a) in § 3.2; this also reinforces the trend of finite-extensibility-induced destabilisation at higher (but finite) $E$ values already found in § 3.1. This trend is explored more systematically in figure 13(a), wherein the elastic centre mode is tracked as $L$ is decreased from the Oldroyd-B limit ( $L = \infty$ ) down to $L = 40$ , again for $\beta = 0.98$ . The centre mode is unstable for $45 \lt L \lt 110$ , corresponding to $0.225 \lt W/L \lt 0.55$ . Figure 13(b) shows that in contrast to the sensitive change in the stability characteristics, there is little variation in the base-flow velocity profiles over the said range of $L$ values. The velocity and stress eigenfunctions for the Oldroyd-B and FENE-P cases are shown in figure 14. Similar to the elastoinertial case discussed earlier (figure 5), one finds the velocity eigenfunctions to be qualitatively similar for both FENE-P and Oldroyd-B fluids, with significant differences discernible only in the $\tilde {\tau }_{xx}$ and $\tilde {\tau }_{xz}$ eigenfunctions, and with differences in the former being more pronounced.

Figure 13.

(a) Elastic spectra: the enlarged spectrum (near the centre mode) showing the variation in the centre mode as $L$ is varied at $\beta = 0.98$ , $k = 2$ , $W=200$ and $Re = 0$ . (b) The base flow velocity profiles for the same $\beta$ and $W$ , and for different $L$ values.

Figure 14.

Velocity, stress and pressure eigenfunctions corresponding to the elastic centre mode in figure 13 for $\beta = 0.98$ , $k = 2$ and $W=200$ , in the Oldroyd-B limit and for $L=100$ . The corresponding eigenvalues are: $L = 100$ ( $U_{max} = 1.01254$ ), $c = 1.0120157 + 0.000102{\rm i}$ ; Oldroyd-B ( $U_{max} = 1$ ), $c = 0.999277 - 0.000850{\rm i}$ .

Figure 15.

The elastic centre-mode instability within the FENE-P framework. (a) The neutral curves in the $W$ – $k$ plane for different $\beta$ values for $L = 100$ . (b) The variation of the critical Weissenberg number $W_c$ with $(1-\beta )$ for different $L$ values. (c) An expanded view of the region of (b) for $40 \leqslant L \leqslant 200$ .

The neutral curves for the elastic centre-mode instability, in the $W$ – $k$ plane, are shown for different $\beta$ values (for $L = 100$ ) in figure 15(a). Similar to the $Re$ – $k$ neutral curves for the elastoinertial instability discussed earlier, unstable regions take the form of islands, but now in the $W$ – $k$ plane. The reason for the island, as opposed to the tongue for the Oldroyd-B case (see Khalid et al. Reference Khalid, Shankar and Subramanian2021b ), is again the same as in the elastoinertial case – shear thinning. The variation of the critical Weissenberg number $W_c$ (minima of the neutral curves in the $W$ – $k$ plane) with $(1-\beta )$ is shown for different $L$ values in figure 15(b), along with the result obtained earlier using the Oldroyd-B model by Khalid et al. (Reference Khalid, Shankar and Subramanian2021b ). There is a systematic decrease of $W_c$ with decreasing $L$ , with the critical $W_c$ required for instability coming down to $W_c \lt 200$ for $L \approx 100$ ; eventually, for small enough $L\ (\lesssim 40)$ , the $W_c$ curve starts to move up again, as evident from the magnified view in figure 15(c). The lowest critical value $W_c \approx 163$ for $\beta = 0.975$ is consistent with the value reported in the study of Buza et al. (Reference Buza, Page and Kerswell2022). Both figures 15(a) and 15(b) also serve to show that the elastic instability exists over a wider range of $\beta$ ( $0.97 \lesssim \beta \lt 1$ ), compared to the Oldroyd-B model ( $0.990552 \leqslant \beta \lt 1$ ). In both cases, $W_c$ eventually increases for $\beta$ sufficiently close to unity, owing to the absence of the elastic instability in the Newtonian limit ( $\beta = 1$ ).