2.1. Governing equations

We consider steady, fully developed pressure-driven flow of a viscoelastic fluid between concentric cylindrical surfaces with radii $R_1$ and $R_2$ ( $R_1 \lt R_2$ ), corresponding to a gap width $d = R_2 - R_1$ ; the gap-width ratio is defined as $\epsilon = d/R_1$ . A cylindrical coordinate system with $r$ , $\theta$ , $z$ denoting the radial, azimuthal and axial coordinates, respectively, shown in figure 1, is used with the pressure gradient imposed along the azimuth. The scales used for rendering dimensionless the physical quantities, associated with both Oldroyd-B and FENE-P models, are summarised in table 1. The dimensionless continuity and inertialess Cauchy momentum equations take the form

(2.1) \begin{equation} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{v}=0 , \end{equation}

(2.2) \begin{equation} 0 =-\boldsymbol{\nabla }\!{p}+ {\beta } \,{\nabla} ^2 \, \boldsymbol{v} + {(1-\beta )} \, \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{\tau } , \end{equation}

with $\boldsymbol{v}$ and $p$ being the velocity and pressure fields, and $\boldsymbol{\tau }$ the polymeric stress tensor. In (2.2), $\beta = \eta _s/\eta$ is the solvent to solution viscosity ratio, where $\eta = \eta _s + \eta _p$ is the total (zero-shear) solution viscosity with $\eta _p$ and $\eta _s$ being the zero-shear polymer and solvent viscosities, respectively.

In the Oldroyd-B model (Larson Reference Larson1988), the polymer solution is modelled as a non-interacting suspension of Hookean dumbbells and the polymeric stress $\boldsymbol{\tau }$ is governed by

(2.3) \begin{eqnarray} \boldsymbol{\tau }+\textit{Wi} \left(\frac {\partial \boldsymbol{\tau }}{\partial t}+ \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{\tau }-\boldsymbol{\nabla }\boldsymbol{v}^\intercal \boldsymbol{\cdot }\boldsymbol{\tau }- \,\boldsymbol{\tau }\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v} \right)= \big (\boldsymbol{\nabla }\boldsymbol{v}+\boldsymbol{\nabla }\boldsymbol{v}^\intercal \big ). \end{eqnarray}

Here, $\textit{Wi} = \lambda U_m/d$ is the Weissenberg number. In the FENE-P model, the dumbbells are finitely extensible and the polymeric contribution to the stress tensor is given by

(2.4) \begin{equation} \boldsymbol{\tau } = \frac {f \boldsymbol{C} - \boldsymbol{I}}{\textit{Wi}}. \end{equation}

Here, $\boldsymbol{C}$ denotes the dimensionless conformation tensor, while $\boldsymbol{I}$ is the identity tensor representing the isotropic distribution of dumbbell conformations in the absence of an imposed flow. The Peterlin closure function $f$ (Bird et al. Reference Bird, Dotson and Johnson1980; Herrchen & Ottinger Reference Herrchen and Ottinger1997) is defined as

(2.5) \begin{equation} f = \frac {L^2 - 3}{L^2 -tr(\boldsymbol{C})}, \end{equation}

where $L$ represents the dimensionless extensibility parameter, de’s maximum extension to its equilibrium root-mean-square value; the Oldroyd-B model being attained in the limit of infinite extensibility viz., $L \rightarrow \infty$ , when $f \rightarrow 1$ . The conformation tensor is governed by the following evolution equation:

(2.6) \begin{equation} \frac {\partial \boldsymbol{C}}{\partial t}+ \boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{C}-\boldsymbol{\nabla }\boldsymbol{v}^\intercal \boldsymbol{\cdot }\boldsymbol{C}- \,\boldsymbol{C}\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{v} = -\boldsymbol{\tau }. \end{equation}

No-slip boundary conditions at the inner and outer cylinder walls apply to all velocity components. Since stress diffusion is neglected in the constitutive equations, no boundary conditions need be specified for either the polymeric stress or conformation tensors.

2.2. Base state for the Oldroyd-B Model

The base state velocity profile for Dean flow of an Oldroyd-B fluid is identical to its Newtonian counterpart (Diprima Reference Diprima1959; Joo & Shaqfeh Reference Joo and Shaqfeh1992) and, therefore, independent of $\textit{Wi}$ . The non-dimensional velocity profile, after using the coordinate transformation $r^{*} = R_1(1 + \epsilon y)$ , is given by

(2.7) \begin{equation} \overline {\boldsymbol{v}} = \left [\!\begin{array}{ccccc} \overline v_{r} \\[3pt] \overline v_{\theta } \\[3pt] \overline v_{z} \\ \end{array}\!\right ] = \left [\!\begin{array}{ccccc} 0 \\[3pt] U(y) \\[3pt] 0 \\ \end{array}\!\right ] \!, \end{equation}

(2.8) \begin{equation} U(y) = \frac { -K \left [ (1 + \epsilon y) - \frac {1}{(1 + \epsilon y)}\right ] + \left [(1 + \epsilon y) \ln (1 + \epsilon y) \right ] }{K \left [\frac { \ln (1 + \epsilon )}{\epsilon }\right ] - \frac {1}{4} (2 + \epsilon )} , \end{equation}

where $K = {(1 + \epsilon )^2 \ln (1+ \epsilon )}/{(1 + \epsilon )^2 - 1}$ and the average velocity that has been used as the scale in (2.8) is given by $U_m = ( {( {\partial p}/{\partial \theta }) R_1}/{2 \eta }) (K [ { \ln (1 + \epsilon )}/{\epsilon } ] - ({1}/{4}) (2 + \epsilon ) )$ . The polymer contribution to the base state stress tensor is given by

(2.9) \begin{equation} \boldsymbol{\overline \tau } = \left [\!\!\begin{array}{ccccc} \overline \tau _{rr} & \overline \tau _{r \theta } & \overline \tau _{rz} & \\[3pt] \overline \tau _{\theta r} & \overline \tau _{\theta \theta } & \overline \tau _{\theta z} & \\[3pt] \overline \tau _{zr } & \overline \tau _{z \theta } & \overline \tau _{zz} & \\ \end{array}\!\!\!\!\!\right ] = \left [\!\begin{array}{ccccc} 0 & U' - \frac {\epsilon }{1 + \epsilon y}U & 0 & \\[11pt] U' - \frac {\epsilon }{1 + \epsilon y}U & 2 \textit{Wi} \left(U' - \frac {\epsilon }{1 + \epsilon y}U\right)^2 & 0 & \\[11pt] 0 & 0 & 0 & \\ \end{array}\!\!\right ] \!, \end{equation}

where $ \zeta ' \equiv D\zeta \equiv{{\rm d}\zeta }/{{\rm d}y}$ .

2.3. Base state for the FENE-P model

The base state velocity field, conformation tensor and the polymer stress tensor in this case are given by

(2.10) \begin{equation} \overline {\boldsymbol{v}} = \left [\! \begin{array}{ccccc} \overline v_{r} \\[3pt] \overline v_{\theta } \\[3pt] \overline v_{z} \\ \end{array} \!\right ] = \left [\! \begin{array}{ccccc} 0 \\[3pt] U(y) \\[3pt] 0 \\ \end{array}\!\right ] \!, \end{equation}

(2.11) \begin{equation} \boldsymbol{\overline C} = \left [\!\! \begin{array}{ccccc} \overline C_{rr} & \overline C_{r\theta } & \overline C_{rz} & \\[5pt] \overline C_{\theta r} & \overline C_{\theta \theta } & \overline C_{\theta z} & \\[5pt] \overline C_{zr} & \overline C_{z\theta } & \overline C_{zz} & \\ \end{array}\!\!\!\!\!\!\right ] = \left [\! \begin{array}{ccccc} \frac {1}{\bar f} & \frac {\textit{Wi}\, \overline C_{rr} \left(U'- \frac {U \epsilon }{1+ \epsilon y}\right)}{\bar f} & 0 & \\[11pt] \frac {\textit{Wi}\, \overline C_{rr} \left(U'- \frac {U \epsilon }{1+ \epsilon y}\right)}{\bar f} & \frac {1 + 2 \textit{Wi}\, \overline C_{r\theta }\left( U' -\frac {U \epsilon }{1+ \epsilon y}\right)}{\bar f} & 0 & \\[11pt] 0 & 0 & \frac {1}{\bar f} \\ \end{array}\!\!\!\right ]\!, \end{equation}

(2.12) \begin{equation} \boldsymbol{\overline \tau } = \left [\!\!\begin{array}{ccccc} \overline \tau _{rr} & \overline \tau _{r \theta } & \overline \tau _{rz} & \\[5pt] \overline \tau _{\theta r} & \overline \tau _{\theta \theta } & \overline \tau _{\theta z} & \\[5pt] \overline \tau _{zr } & \overline \tau _{z \theta } & \overline \tau _{zz} & \\ \end{array}\!\!\!\!\!\right ] = \left [ \begin{array}{ccccc} 0 & \overline C_{rr} \left(U'- \frac {U \epsilon }{1+ \epsilon y}\right) & 0 & \\[11pt] \overline C_{rr} \left(U'- \frac {U \epsilon }{1+ \epsilon y}\right) & 2 \overline C_{r\theta } \left(U'- \frac {U \epsilon }{1+ \epsilon y}\right) & 0 & \\[11pt] 0 & 0 & 0 \\ \end{array}\!\!\!\!\right ]\!. \end{equation}

By combining (2.5) in the base state with $\overline \tau _{r \theta } =\, \overline C_{rr} (U'- ( {U \epsilon }/{1+ \epsilon y})) = ({U' - ( {U \epsilon }/ {1 + \epsilon y})}/{\bar f})$ and eliminating $\bar f$ , we obtain

(2.13) \begin{equation} \tau _{r \theta }^3 + \frac {L^2}{2 \textit{Wi}^2} \tau _{r \theta } - \frac {L^2}{2 \textit{Wi}^2} \left (U' - \frac {U \epsilon }{1+ \epsilon y}\right ) = 0 . \end{equation}

This equation for $\tau _{r \theta }$ is coupled to $U(y)$ through the $\theta$ -component of (2.2), which provides a second relation between $\tau _{r \theta }$ and $U(y)$ . The two equations are simultaneously solved, numerically, to obtain the base state fields; more details regarding the numerical procedure can be found from Tej et al. (Reference Tej, Mohanty and Shankar2024). The base state velocity profile does depend on $\textit{Wi}$ and $L$ , although only via the combination $\textit{Wi}/L$ for $L \gg 1$ . Figure 2(a) illustrates the collapse of the profiles for different ( $\textit{Wi}$ , $L$ ) pairs, with $\textit{Wi}/L$ fixed (and for $L\gg 1$ ). This collapse, first demonstrated by Yamani & McKinley (Reference Yamani and McKinley2023) for simple shear flow, was later generalised by Tej et al. (Reference Tej, Mohanty and Shankar2024) to steady unidirectional rectilinear and curvilinear shear flows of a FENE-P fluid. Figure 2(b) shows velocity profiles for different $\textit{Wi}/L$ , for a fixed $\beta$ . Interestingly, the maximum velocity, scaled by the average, is seen to exhibit a non-monotonic dependence on $\textit{Wi}/L$ (see inset of figure 2 b). For small $\textit{Wi}/L$ , the profile expectedly assumes the Newtonian parabolic form and $U_{\textit{max}}$ equals $({3}/{2}) \, U_{avg}$ ; for sufficiently large $\textit{Wi}/L$ , shear thinning ensures that the solvent contribution is dominant and $U_{\textit{max}}$ again approaches $( {3}/{2}) \, U_{avg}$ . For intermediate $\textit{Wi}/L$ values, shear thinning leads to flatter profiles and a smaller $\beta$ -dependent $U_{\textit{max}}$ . The aforementioned effect of changing $\beta$ at a fixed $\textit{Wi}/L$ and $\epsilon$ is shown in figure 2(c), where shear thinning again leads to $U_{\textit{max}}/U_{avg} \lt 1.5$ , for the smaller $\beta$ values; the profile reverts to the Newtonian form for $\beta \rightarrow 1$ . Finally, the role of varying $\epsilon$ is shown in figure 2(d) where, for $\epsilon \ll 1$ (narrow-gap limit), the maximum velocity occurs at the centreline; this location shifting towards the inner cylinder for finite gap-width ratios.

Figure 2.

Dean flow velocity profiles for a FENE-P fluid: (a) $\beta = 0.98$ and $\epsilon = 0.1$ for two $(\textit{Wi}, L)$ pairs, with $\textit{Wi}/L = 10$ ; (b) $\beta = 0.5$ and $\epsilon = 0.1$ with varying $\textit{Wi}/L$ . The inset shows the variation of the ratio of maximum to average velocities with $\textit{Wi}/L$ ; (c) master curves for $\textit{Wi}/L = 10$ and $\epsilon = 0.1$ for varying $\beta$ ; (d) master curves for $\textit{Wi}/L = 10$ and $\beta = 0.9$ for varying $\epsilon$ .

2.4. Linear stability analysis

A temporal linear stability analysis is carried out wherein the aforementioned base states ((2.7)–(2.9) for the Oldroyd-B model and (2.10)–(2.12) for the FENE-P model) are subjected to small-amplitude perturbations. We consider both axisymmetric and non-axisymmetric disturbances in this study. The total velocity, pressure and stress are each expressed as the sum of a base-state contribution ( $\overline {\chi }$ ) and an imposed perturbation ( $\hat {\chi }$ ):

(2.14) \begin{equation} \chi = \overline {\chi } + \hat \chi , \end{equation}

where the variable $\chi$ represents any of the field variables. Next, the perturbation fields are assumed to be of the Fourier mode form:

(2.15) \begin{equation} \hat {\chi }(y,z,\theta ,t) = \skew 2\tilde {\chi }(y) \exp \left [\textrm{i} \left (\alpha z + n \theta - \frac {\omega }{\textit{Wi}} t \right )\right ]\!. \end{equation}

Here, $\alpha$ and $n$ are the axial and azimuthal wavenumbers, respectively, and $\omega = \omega _r + i \omega _i$ is the complex frequency. Since the frequency $\omega$ has been rendered dimensionless using the relaxation time and time $t$ using the convective time scale $d/U_m$ , then $\omega$ appears via the above-mentioned combination $\omega /\textit{Wi}$ . Axisymmetric disturbances correspond to $n = 0$ , while $\alpha = 0$ corresponds to two-dimensional perturbations varying in the $r$ – $\theta$ plane; three-dimensional disturbances correspond to $\alpha \neq 0$ and $n \neq 0$ . The flow is temporally unstable if $\omega _i \gt 0$ . Substituting (2.15) in the linearised versions of (2.1), (2.2) and (2.3) leads to the set of governing equations, given in Appendices A and B for the Oldroyd-B and FENE-P models, respectively.

2.5. Numerical methods

We solved the linearised governing equations using a spectral method (Trefethen Reference Trefethen2000; Weideman & Reddy Reference Weideman and Reddy2000) and a shooting procedure. In the spectral method, we expand the dynamical variables as a finite sum of Chebyshev polynomials and substitute this expansion in the linearised governing equations – (A1)–(A10) for the Oldroyd-B model or (B1)–(B4) for the FENE-P model – to obtain a generalised eigenvalue problem as follows:

(2.16) \begin{equation} \boldsymbol{A} \boldsymbol{x} = \omega \boldsymbol{B} \boldsymbol{x} . \end{equation}

Here, $\boldsymbol{A}$ and $\boldsymbol{B}$ are coefficient matrices, with $\boldsymbol{x} = (\tilde v _r , \tilde v _\theta , \tilde v _z, \tilde p, \tilde \tau _{rr}, \tilde \tau _{r\theta }, \tilde \tau _{rz}, \tilde \tau _{\theta \theta }, \tilde \tau _{\theta z}, \tilde \tau _{zz} )^T$ for the Oldroyd-B model and $\boldsymbol{x} = (\tilde v _r , \tilde v _\theta , \tilde v _z, \tilde p, \tilde C _{rr}, \tilde C_{r\theta }, \tilde C_{rz}, \tilde C_{\theta \theta }, \tilde C _{\theta z}, \tilde C_{zz} )^T$ for the FENE-P model. The size of the $\boldsymbol{A}$ matrix is $10N \times 10N$ , where $N$ is the number of Gauss–Lobatto collocation points. The generalised eigenvalue problem is solved using the ‘eig’ solver of Matlab. Depending on the parameter regime, the number of collocation points $N$ required to capture physically genuine eigenvalues differs. Convergence in this method is achieved by comparing the eigenspectra for two different values of $N$ and filtering them using a prescribed tolerance. The shooting method is based on a Runge–Kutta integrator with a Gram–Schmidt orthonormalisation procedure (Ho & Denn Reference Ho and Denn1977) and a Newton–Raphson iteration for the eigenvalue. This procedure is implemented only for the Oldroyd-B model. It is known that the spectral method can give unphysical (spurious) eigenvalues (Schmid & Henningson Reference Schmid and Henningson2001), but the shooting method gives genuine eigenvalues, given a good initial guess. Thus, the genuineness of the eigenvalues, obtained from the spectral method, is confirmed by providing them as initial guesses for the shooting method. To benchmark the implementation of our numerical methodology in the narrow-gap limit, we compare (table 2 and figure 3) results from our procedure with those of Joo & Shaqfeh (Reference Joo and Shaqfeh1992), finding excellent agreement.

Table 2.

Validation of UCM $(\beta = 0)$ results with those of Joo & Shaqfeh (Reference Joo and Shaqfeh1992).

Table 2. Long description

The table presents a comparison of Dean flow parameters, specifically alpha_c and epsilon to the power of 1/2 times Wi_c, with results from Joo and Shaqfeh and present findings. The table has three rows and three columns. The columns are labeled ‘Dean flow’, ‘alpha_c’, and ‘epsilon^1/2 Wi_c’. The rows are labeled ‘Dean flow’, ‘Joo and Shaqfeh’, and ‘Present’. Row 1: Dean flow, alpha_c, epsilon^1/2 Wi_c. Row 2: Joo and Shaqfeh, 6.6 ± 0.1, 4.06 ± 0.02. Row 3: Present, 6.56, 4.06.

Figure 3.

Benchmarking of neutral curves obtained from our numerical (spectral-cum-shooting) procedure (for $\textit{Re} = 0$ , $ n = 0$ ) with those of Joo & Shaqfeh (Reference Joo and Shaqfeh1992) for Dean flow of an Oldroyd-B fluid, for different $\beta$ . Continuous lines show results from the present work, while discrete points represent the data of Joo & Shaqfeh (Reference Joo and Shaqfeh1992).

Figure 4.

Eigenspectrum showing (a) stationary ( $\alpha = 7$ ) and (b) propagating ( $\alpha = 14$ ) modes of instability via the HSM1 mode in Dean flow of an Oldroyd-B fluid at $\textit{Re} = 0$ , $n = 0$ , $\beta = 0.98$ , $\epsilon = 0.1$ and $\textit{Wi} = 25$ .

3.1. Axisymmetric disturbances ( $n = 0$ )

Figure 4 shows the axisymmetric elastic spectrum for Dean flow of an Oldroyd-B fluid at two different axial wavenumbers ( $\alpha$ ). For both $\alpha$ , there appear to be a large number of discrete modes accumulating towards the point $(\omega _r,\omega _i) \equiv (0,-1)$ along multiple branches; the latter point is shown here to be a part of the continuous spectrum (CS) for $n = 0$ . The numerical eigenvalues are likely a large but finite representation of infinite sequences, pertaining to the analytical eigenvalue problem, that converge to the aforementioned point CS. In this regard, we have verified that progressively more discrete eigenvalues are resolved in the neighbourhood of the CS, with no indication of the sequence terminating with increasing $N$ . The spectra in figures 4(a) and 4(b) also include unstable modes – either stationary ( $\omega _r = 0$ ; see figure 4 a) or propagating ( $\omega _r \neq 0$ ; see figure 4 b) ones; in the latter case, they appear symmetrically on either side of the imaginary axis. For the stationary scenario, the least stable/most unstable of the discrete modes is the classical axisymmetric HSM, originally identified by Joo & Shaqfeh (Reference Joo and Shaqfeh1991, Reference Joo and Shaqfeh1994), and is labelled ‘HSM1’ henceforth. A more detailed discussion of the parameter regimes in which HSM1 is stationary or propagating is provided in Appendix C. It will be seen therein that the stationary-to-propagating transition can happen with an increase in $\alpha$ for a fixed $\textit{Wi}$ or for a fixed $\alpha$ and with decreasing $\textit{Wi}$ .

Figure 5.

Schematic of the continuous spectra branches in Dean flow of a FENE-P fluid with $\textit{Wi}/L \sim O(1)$ and $\beta \rightarrow 1$ : (a) $ n = 0$ ; (b) $ n \neq 0$ .

Now, it is well known that the eigenspectra for viscoelastic flows also contain (stable) continuous spectra of eigenvalues (Wilson et al. Reference Wilson, Renardy and Renardy1999; Roy et al. Reference Roy, Garg, Reddy and Subramanian2022). These are obtained by setting the coefficient of the highest order derivative of the governing linearised equation for the wall-normal velocity perturbation to zero. For rectilinear flows such as plane Couette and Poiseuille flows of an Oldroyd-B fluid, this leads to two distinct CS with decay rates $\omega _i = -1$ and $\omega _i = -1/\beta$ . In general, these CS are horizontal line segments in the $\omega _r$ – $\omega _i$ plane whose length is proportional to $k \Delta V$ , with $\Delta V$ being the base range of velocities and $k$ the streamwise wavenumber. Along similar lines, for Dean flow, the horizontal extent of the CS is proportional to $n \Delta {V}$ . For Dean flow subjected to axisymmetric disturbances ( $n = 0$ ) as in figure 4, the CS, within the Oldroyd-B framework, therefore degenerate to a pair of points on the imaginary axis with $\omega _i = -1$ and $-1/\beta$ ; the the two points approaching each other in the limit $\beta \rightarrow 1$ ( $\beta = 0.98$ in figure 4).

We now turn to the effect of $L$ on the spectra. As $L$ is decreased, owing to shear thinning of rheological properties, the nature of CS for both rectilinear and curvilinear shearing flows, including Dean flow in particular, is significantly more complex compared with that for the above-mentioned Oldroyd-B fluid. In the interest of brevity, we provide here only the key features of the CS that are relevant to the present work; a detailed discussion can be found from Mohanty et al. (Reference Mohanty, Tej, Subramanian and Shankar2026). A schematic of the CS for Dean flow, pertinent to the dilute-solution limit ( $\beta \gtrsim 0.95$ ) and for $\textit{Wi}/L \sim O(1)$ (so shear thinning effects remain significant), is shown in figure 5(a) for $n = 0$ . As $L$ is decreased from the Oldroyd-B limit, the ‘point CS’ on the imaginary axis, in the above-mentioned Oldroyd-B limit, turn into vertical line segments. The two line segments are shown to be distinct in the schematic, but may overlap depending on the particular choice of parameters. The reason underlying the extended CS, even for $n = 0$ , is due to the variation of the base-state Peterlin function $\bar {f}$ as a function of the wall-normal coordinate, which leads to a range of relaxation times and, thence, decay rates. Concurrently, a pair of symmetrically placed ‘wing’-like CS emerge on either side of the imaginary axis (Mohanty et al. Reference Mohanty, Tej, Subramanian and Shankar2026). The horizontal extent of the wings and the vertical line segments (belonging to the CS) increases with decreasing $L$ . The computed numerical spectra illustrating these features are shown in figure 6 for different $L$ at a fixed $\textit{Wi}$ .

Apart from the CS, the spectra in figure 6 also contain discrete modes, which are not depicted in the figure 5(a) schematic. In fact, figure 6(a–h) shows that decreasing $L$ has a strong influence on the configuration of discrete modes. In contrast to the previous Oldroyd-B case, where there appeared to be multiple infinite sequences of discrete modes accumulating at the (point) CS, a decrease in $L$ leads to a progressive reduction in the number of modes. When $L$ is reduced to $398$ , the most unstable mode changes from a stationary one to a pair of propagating ones; this change is due to a coalescence of HSM1 and the next (least stable) mode, followed by a bifurcation. The propagating pair moves down towards the CS as $L$ is decreased further; this shear-thinning-induced stabilisation is discussed further in § 5.1. Eventually, for $L = 100$ , there are no discrete modes, either stable or unstable, and the spectrum comprises solely modes belonging to the CS for $\textit{Wi} \gtrsim 16$ (see figure 6 h).

Figure 6.

Eigenspectra for Dean flow of FENE-P fluids at different $L$ ; data for $n = 0$ , $\epsilon = 0.1$ , $\alpha = 7$ , $\beta = 0.98$ , $\textit{Wi} = 20$ and $\textit{Re} = 0$ : (a) $ L = 10^4$ (Oldroyd-B); (b) $ L = 5000$ ; (c) $ L = 1000$ . (d) $ L = 500$ ; (e) $ L = 398$ ; (f) $ L = 390$ ; ( g) $ L = 150$ and (h) $ L = 100$ .

Figure 7.

Eigenspectrum for non-axisymmetric ( $n = 1$ ) disturbances showing both HSM1 and HSM2 in Dean flow of Oldroyd-B ( $\textit{Wi}=20$ ) and FENE-P ( $L = 200,\textit{Wi}=400$ ) fluids. Data for $\alpha = 7$ , $\beta = 0.98$ and $ \epsilon = 0.1$ .

Figure 7. Long description

Two scatter plots depict the eigenspectrum for non-axisymmetric disturbances in Dean flow of Oldroyd-B and FENE-P fluids. Panel A: The scatter plot shows the eigenspectrum for Oldroyd-B fluid with Wi equal to 20. The horizontal axis is labeled ωr and ranges from -0.5 to 3.5. The vertical axis is labeled ωi and ranges from -1.4 to 0.2. The plot includes data points for N equal to 180 and N equal to 200, represented by magenta circles and green asterisks, respectively. Key features include HSM1, HSM2, and CS, each labeled with arrows pointing to their respective locations on the plot. Panel B: The scatter plot shows the eigenspectrum for FENE-P fluid with L equal to 200 and Wi equal to 400. The horizontal axis is labeled ωr and ranges from -10 to 70. The vertical axis is labeled ωi and ranges from -50 to 20. The plot includes data points for N equal to 180 and N equal to 200, represented by magenta circles and green asterisks, respectively. Key features include HSM1, HSM2, and CS, each labeled with arrows pointing to their respective locations on the plot.

3.2. Non-axisymmetric disturbances ( $n \neq 0$ )

The schematic of the CS for $n \neq 0$ is shown in figure 5(b) for FENE-P fluids with $\textit{Wi}/L \sim O(1)$ ; again, the details of how these are obtained theoretically are discussed by Mohanty et al. (Reference Mohanty, Tej, Subramanian and Shankar2026). There are four distinct CS (for $L \gtrsim 50$ ) in figure 5(b), with each of these comprising two branches. The branches arise due to the absence of a mid-plane symmetry in Dean flow and degenerate in the case of plane-Poiseuille flow. The actual numerical spectra for a non-zero $n$ ( $n = 1$ ) for both the Oldroyd-B and FENE-P cases, are shown in figures 7(a) and 7(b). In figure 7(a), we present the spectra for the same set of parameters as for the axisymmetric case ( $n = 0$ ) in figure 6(a). For $n = 1$ , the theoretical CS are no longer points even in the Oldroyd-B fluid, but are horizontal line segments located again at $\omega _i = -1$ and $-1/\beta$ ; these are only resolved as balloons in figure 7(a).

The constitution of the discrete modes in the spectra for $n = 1$ also differs qualitatively from those for $n = 0$ . In particular, the discrete modes in figure 6(a) (for $n = 0$ ) with negative $\omega _r$ are virtually absent in the spectrum for $n = 1$ (see figure 7 a). The absence is understandable since a negative $\omega _i$ implies a propagation in the direction opposite to the base state shear flow, which is an unlikely possibility in the absence of inertia. We have verified (data not shown) that as $n$ is increased continuously from $0$ to $1$ , there is also an overall reduction in the number of discrete modes; this reduction arising because many of these discrete modes disappear into the CS. In figure 7(a), the classical hoop-stress mode (HSM1) is the mode that dictates instability in the Oldroyd-B limit. This mode could be stationary or propagating for axisymmetric disturbances, but is propagating for non-axisymmetric disturbances.

Amidst the collection of other discrete modes, we also identify another mode (labelled ‘HSM2’), which remains stable at the $\textit{Wi}$ chosen in figure 7(a). One reason behind singling out HSM2, amongst other discrete modes, maybe seen in figure 7(b), where it is unstable with a growth rate that is nearly equal to that of HSM1, for $L = 200$ , albeit at a higher $\textit{Wi} = 400$ . Figure 8(a–d) shows the complete spectra for the Oldroyd-B fluid at different $\textit{Wi}$ , with both HSM1 and HSM2 identified, while figure 9 tracks the loci of HSM1 and HSM2 in the complex- $\omega$ plane with increasing $\textit{Wi}$ . The inset in this figure shows the variation of the growth rate with $\textit{Wi}$ for both modes. While an increase in $\textit{Wi}$ has a monotonic destabilising influence on HSM1, it has a non-monotonic effect on HSM2, with this mode being unstable only over an intermediate range of $\textit{Wi}$ , and eventually stabilising at sufficiently high $\textit{Wi}$ ; the growth rates of HSM2 are much smaller compared with those of HSM1.

Returning to figure 7(b), the spectrum (for $L = 200$ ) also demonstrates the presence of multiple non-trivial CS branches, for $n \neq 0$ , which nevertheless conform to the schematic shown in figure 5(b). While figure 7(b) is for $L = 200$ , HSM2 is further destabilised with decreasing $L$ , leading to it becoming the critical mode for sufficiently small $L$ , in certain regions of the parameter space. This is in contrast to the Oldroyd-B limit, where HSM1 is always the critical mode. The detailed behaviour of HSM1 and HSM2 growth rates, for the FENE-P case, is discussed in § 5.2. It is nevertheless worth a brief mention here, since the unexpected destabilising effect of finite elasticity first seen for the CM by Buza et al. (Reference Buza, Page and Kerswell2022b ), is also present for HSM2.

Figure 8.

Eigenspectra showing HSM1 and HSM2 in Dean flow of an Oldroyd-B fluid at $\textit{Re} = 0$ , $\alpha = 7$ , $\beta = 0.98$ , $ \epsilon = 0.1, n = 1$ for different $\textit{Wi}$ : (a) $\textit{Wi} = 20$ ; (b) $\textit{Wi} = 60$ ; (c) $\textit{Wi} = 200$ ; (d) $\textit{Wi} = 400$ .

Figure 9.

Variation of HSM1 and HSM2 growth rates with increasing $\textit{Wi}$ . Note the difference in scales of the left and right $y$ -axes (for HSM1 and HSM2, respectively) on account of the vast difference in the growth rate magnitudes. The arrows indicate the direction of increasing $\textit{Wi}$ . The inset shows the variation of $\omega _i$ with $\textit{Wi}$ for both the modes. Data for $\textit{Re} = 0$ , $\alpha = 7$ , $\beta = 0.98$ , $ \epsilon = 0.1, n = 1$ .

Figure 9. Long description

A line graph with two data lines representing HSM1 and HSM2 growth rates. The horizontal axis is labeled with omega subscript r and ranges from 10^-1 to 10^2. The left vertical axis is labeled with omega subscript i and ranges from 0 to 250, while the right vertical axis is labeled with omega subscript i and ranges from -0.6 to 1.2. The black line represents HSM1, and the red line represents HSM2. Arrows indicate the direction of increasing Wi. The inset shows a smaller line graph with the same axes labels and ranges, depicting the variation of omega subscript i with Wi for both modes. The data lines in the inset show different trends for HSM1 and HSM2 as Wi increases.

We end our discussion with figure 10, which presents the elastic spectrum for $\beta = 0.994$ and $\epsilon = 0.1$ ; the spectrum now includes an unstable CM. The contrast in the choice of parameters, relative to the spectra in figures 4 and 6, is on account of the restricted nature of the unstable region in parameter space, in the limit of $\epsilon \rightarrow 0$ , corresponding to plane Poiseuille flow – the centre-mode instability only exists for $\beta \gt 0.99052$ in this limit and for $\textit{Wi} \sim O(300)$ (within the Oldroyd-B model; see Khalid et al. Reference Khalid, Shankar and Subramanian2021b ). It is pertinent to note here that the Weissenberg number used in the present study and the one used in the earlier effort of Khalid et al. (Reference Khalid, Shankar and Subramanian2021b ) are related by a factor of three, as we elaborate in § 4.1. The role of decreasing $L$ on the CM is discussed subsequently in § 4.2.

Figure 10.

Eigenspectrum showing the presence of the centre mode in Dean flow of an Oldroyd-B fluid (with $n = 16$ , $\epsilon = 0.1$ , $\alpha = 0$ , $\beta = 0.994$ and $\textit{Wi} = 350$ ). The inset shows zoomed-in view highlighting the region near the unstable centre mode; the discrete eigenvalues obtained for $N = 600$ and $N = 700$ show good convergence.

4.1. Centre-mode instability within the Oldroyd-B framework

In the narrow-gap limit, $\epsilon \ll 1$ , the base state Dean profile approaches plane Poiseuille flow and we therefore begin with a validation of our numerical procedure by recovering the expected results (Khalid et al. Reference Khalid, Shankar and Subramanian2021b ) in this limit. To establish the connection between the two flows, note that the angular displacement $\theta$ in the Dean (cylindrical) geometry is the ratio of the arc length $x^*$ to the radial distance $r^*$ , that is, $\theta = {x^{*}}/{r^{*}}$ . In the narrow-gap limit, $r^* = R_1(1 + \epsilon y) \approx R_1$ , and the previous relation reduces to $\theta = x^*/R_1$ , leading to $n\theta = n x^*/R_1$ . Thus, $n/R_1$ in cylindrical coordinates, for Dean flow, corresponds to the dimensional streamwise wavenumber $k^*$ for plane Poiseuille flow; using $k^* = 2k/d$ gives $n = 2k/\epsilon$ . Now, Khalid et al. (Reference Khalid, Shankar and Subramanian2021b ) report the minimum critical Weissenberg number ( $\textit{Wi}_c = 973.8$ ) and critical wavenumber ( $k_c = 0.783$ ) for $\beta = 0.994$ , for the elastic CM in plane Poiseuille flow, with these values arising from use of the maximum speed ( $U_{\textit{max}}$ ) and half-gap width ( $d/2$ ) as the characteristic velocity and length scales; the Weissenberg number and non-dimensional wavenumber being defined as $\textit{Wi}_K = {\lambda U_{\textit{max}}}/{(d/2)}$ and $k = k^{*}d/2$ , respectively. In the present study, the Weissenberg number is defined as $\textit{Wi} = \lambda U_m/d$ . Noting that the ratio of maximum to mean speed is $3/2$ for plane Poiseuille flow, yields the following relation between the critical parameters: $\textit{Wi}/\textit{Wi}_K = 1/3$ , $n_c = 2k_c/\epsilon$ , in turn leading to $\textit{Wi}_c = 324.6$ and $n_c = 1566$ in the narrow-gap limit; such high $n$ do not necessitate an increased number of collocation points, for the spectral results shown in the following, since both the radial ( $\epsilon d$ ) and azimuthal ( $1/n \sim \epsilon$ ) scales decrease in direct proportion to $\epsilon$ .

Figure 11.

Centre-mode instability in Dean flow of an Oldroyd-B fluid at $\textit{Re} = 0$ , with $\alpha = 0$ : Neutral curves in the $\textit{Wi}$ – $n$ plane at different $\epsilon$ for $\beta = 0.994$ . For $\epsilon = 0.001$ and $0.1$ , we also plot the neutral curve (continuous lines) obtained for pressure-driven flow in a rectilinear channel (Khalid et al. Reference Khalid, Shankar and Subramanian2021b ); incipient departure from the $n \sim \epsilon ^{-1}$ scaling is evident in the comparison for $\epsilon = 0.1$ .

Figure 11. Long description

A scatter plot represents the relationship between two variables, Wi and n, on a log-log scale. The horizontal axis represents n, ranging from 10^0 to 10^3. The vertical axis represents Wi, ranging from 10^2 to 10^4. The plot includes data points for different values of epsilon (ε), with ε = 0.001 represented by gray squares, ε = 0.1 by red circles, ε = 0.5 by blue triangles, and ε = 1 by green inverted triangles. Magenta lines indicate trends or fits for the data points. The data points show distinct clusters and patterns for each value of epsilon, with the magenta lines fitting closely to the clusters. The overall trend indicates a negative correlation between Wi and n for each value of epsilon.

Figure 12.

Structure of the CM eigenfunctions in the $y$ – $\theta$ plane for Dean flow of an Oldroyd-B fluid with $\beta = 0.994$ and $\alpha = 0$ ; coloured contours represent the magnitude of the velocity perturbation $|\tilde {v}_{\theta }|$ . Panels (a)–(c) show the progression of the modal structure with increasing curvature, at a constant Weissenberg number ( $\textit{Wi} = 350$ ); panel (d) shows the modal structure at $\epsilon = 1$ for $\textit{Wi} = 230$ .

Figure 11 shows the neutral curves in the $\textit{Wi}$ – $n$ plane for $\epsilon$ ranging from $0.001$ to $1$ . Also shown are the neutral curves obtained by Khalid et al. (Reference Khalid, Shankar and Subramanian2021b ), highlighting the excellent agreement (for $\epsilon = 0.001$ ) between the Dean and channel flow cases. Although it was found convenient to vary $n$ continuously, when tracking the neutral points in the $\textit{Wi}$ – $n$ plane using the shooting procedure, it is worth keeping in mind that the discreteness of $n$ implies that the critical $\textit{Wi}$ corresponds to the integer $n$ closest to the minimum of the tongues in figure 11. There is a leftward movement of the neutral curves with increasing $\epsilon$ , as is expected from the $n \sim \epsilon ^{-1}$ scaling for $\epsilon \ll 1$ ; although this movement slows down for finite $\epsilon$ . There is, in addition, an overall lowering of the unstable tongue, resulting in a decrease in the critical Weissenberg number $\textit{Wi}_c$ with an increase in $\epsilon$ . Figure 12 shows CM eigenfunction contours in the $y$ – $\theta$ plane, over one azimuthal wavelength and for different $\epsilon$ . The plots for $\epsilon = 0.01$ and $0.1$ clearly reveal the ‘arrowhead’ structure identified in earlier theoretical and computational efforts for the channel flow (Page et al. Reference Page, Dubief and Kerswell2020; Khalid et al. Reference Khalid, Chaudhary, Garg, Shankar and Subramanian2021a ), and observed experimentally in flow through straight pipes by Choueiri et al. (Reference Choueiri, Lopez, Varshney, Sankar and Hof2021), and through curved (helical) pipes by Lu & Hof (Reference Lu and Hof2026). The arrowhead is centred around the base-state maximum and tilted in a sense commensurate with the shear of the plane Poiseuille profile. For $\epsilon = 0.5$ and $1$ , curvature effects cause the structures near the inner wall to be significantly diminished in both scale and amplitude.

Figure 13.

Eigenspectra showing the destabilising role of axially varying perturbations on the centre mode in Dean flow of an Oldroyd-B fluid. Data for $n = 1$ , $\epsilon = 0.1$ , $\beta = 0.98$ and $\textit{Wi} = 1500$ . The insets show the enlarged region of the spectra near the centre mode. (a) $\alpha = 0$ ; (b) $\alpha = 0.5$ .

Figure 13. Long description

Two scatter plots depict the eigenspectra of Dean flow in an Oldroyd-B fluid, showing the destabilizing role of axially varying perturbations on the center mode. Panel A: The scatter plot shows the eigenspectra for alpha equals 0. The x-axis represents the real part of the frequency (omega_r) ranging from 0 to 250, and the y-axis represents the imaginary part of the frequency (omega_i) ranging from -250 to 0. The plot includes data points for N equals 400 and N equals 450, with an inset zooming in on the region near the center mode. Panel B: The scatter plot shows the eigenspectra for alpha equals 0.5. The x-axis represents the real part of the frequency (omega_r) ranging from 0 to 250, and the y-axis represents the imaginary part of the frequency (omega_i) ranging from -500 to 0. The plot includes data points for N equals 400 and N equals 450, with an inset zooming in on the region near the center mode.

Figure 14.

(a) Neutral curves in the $\textit{Wi}$ – $\alpha$ plane for two ( $\epsilon$ , $n$ ) pairs for $\beta = 0.994$ . (b) Variation of $\omega _i$ with $\alpha$ showing the destabilisation of the centre mode in Dean flow for non-zero $\alpha$ . Data for the Oldroyd-B fluid with $n = 1$ , $\epsilon = 0.1$ and $\beta = 0.98$ .

Figure 14. Long description

The image contains two line graphs labeled as Panel A and Panel B. Panel A shows neutral curves in the plane for two pairs of parameters (ε, n) for Wi. The x-axis represents α on a logarithmic scale ranging from 10^-2 to 10^1, and the y-axis represents Wi on a logarithmic scale ranging from 10^-2 to 10^3. Two lines are plotted: one for ε = 0.1 & n = 16 in black and another for ε = 1 & n = 3 in red. Panel B illustrates the variation of ωi with α, showing the destabilization of the centre mode in Dean flow for non-zero α. The x-axis represents α ranging from 0 to 1, and the y-axis represents ωi ranging from -1.5 to 3. Three lines are plotted for different Wi values: 1520 in black, 1750 in red, and 2500 in blue. The legend indicates the corresponding Wi values for each line.

Next, figure 13 examines the role of axial variation ( $\alpha \neq 0$ ) of perturbations. In figure 13(a), for the chosen parameter values, the CM is stable for two-dimensional disturbances ( $n \neq 0$ and $\alpha = 0$ ). However, when three-dimensional disturbances are introduced (both $n \neq 0$ and $\alpha \neq 0$ ), surprisingly, the CM becomes unstable (figure 13 b) – the introduction of an axial variation destabilises the centre mode. However, a further increase in $\alpha$ was found to lead to eventual stabilisation. This implies that to obtain the critical $\textit{Wi}$ , it is in principle necessary to minimise over both $n$ and $\alpha$ . To this end, in figure 14(a), we show the neutral curves in the $\textit{Wi}$ – $\alpha$ plane for a fixed $n$ , with the latter taken to be the critical value for $\alpha \rightarrow 0$ . The monotonic increase in $\textit{Wi}_c$ with increasing $\alpha$ , for $\epsilon = 0.1$ , appears to be consistent with the existence of Squire’s theorem in the limit of $\epsilon \rightarrow 0$ (that is, for rectilinear shearing flows; see Bistagnino et al. Reference Bistagnino, Boffetta, Celani, Mazzino, Puliafito and Vergassola2007). That this is not the case is shown in figure 14(b), where we plot the variation of growth rate as a function of $\alpha$ for $\textit{Wi}$ significantly larger than the threshold for neutral stability, but for a $\beta$ ( $= 0.98$ ) that lies below the rectilinear channel-flow threshold for the chosen $\epsilon$ . Even for the relatively low value, $\epsilon = 0.1$ , an increase in $\alpha$ is seen to lead to an unstable CM, where no such instability existed for $\alpha = 0$ . The initial increase in growth rate is, however, followed by stabilisation for sufficiently large $\alpha$ . Returning to figure 14(a), the $\textit{Wi}$ – $\alpha$ neutral curve for $\epsilon = 1$ shows that, while there is an initial destabilisation (a decrease in $\textit{Wi}$ up until $\alpha \approx 2$ , in contrast to $\epsilon = 0.1$ ), this is followed by a stabilising effect at sufficiently large $\alpha$ . On the whole, both the growth rates and the critical $\textit{Wi}$ exhibit a non-monotonic dependence on $\alpha$ . Thus, there exist three-dimensional perturbations that are more unstable than two-dimensional ones, this nevertheless being consistent with the absence of Squire’s theorem for curvilinear geometries. The extent of decrease in $\textit{Wi}_c$ for $\epsilon = 1$ (from $227.76$ for $\alpha = 0$ to $202.03$ for $\alpha = 1.3$ ) is, however, not substantial and $\textit{Wi}_c$ for $\alpha = 0$ serves as a reasonably accurate estimate for the true threshold (obtained by minimising over both $n$ and $\alpha$ ). Thus, for the CM instability, in the ensuing discussion, we present the critical $\textit{Wi}_c$ obtained for two-dimensional disturbances.

Table 3.

Variation of critical Weissenberg number $(\textit{Wi}_c)$ with $\epsilon$ at $ \beta = 0.994$ and $\textit{Re} = 0$ for the centre mode in Dean flow of an Oldroyd-B fluid.

Table 3. Long description

A table with three rows and three columns. The columns are labeled with epsilon, n_c, and W_i_c. The rows are as follows: Row 1: epsilon, 0.001, n_c, 1566, W_i_c, 324.61. Row 2: epsilon, 0.1, n_c, 16, W_i_c, 323.38. Row 3: epsilon, 1, n_c, 3, W_i_c, 227.76.

Figure 15.

Centre-mode instability in Dean flow of an Oldroyd-B fluid at $\textit{Re} = 0$ , $\alpha = 0$ : (a) variation of critical Weissenberg number ( $\textit{Wi}_c$ ) with $\epsilon$ at different $\beta$ ; (b) variation of $\textit{Wi}_c$ with $(1 - \beta )$ at different $\epsilon$ . Here ‘CF’ represents the results for plane Poiseuille flow (Khalid et al. Reference Khalid, Shankar and Subramanian2021b ).

Figure 15(a) shows the reduction in $\textit{Wi}_c$ with increasing $\epsilon$ for different $\beta$ ; representative data, corresponding to this figure, is given in table 3 for completeness. Figure 15(b) shows the variation of the critical Weissenberg number $\textit{Wi}_c$ with $(1-\beta )$ for different $\epsilon$ . The CM instability is seen to be restricted to $\beta \geqslant 0.99$ for $\epsilon \leqslant 0.1$ , consistent with the $\beta$ -threshold of 0.9905 for plane Poiseuille flow (Khalid et al. Reference Khalid, Shankar and Subramanian2021b ). The instability must clearly be absent for $\beta = 1$ and we indeed find that $\textit{Wi}_c$ diverges as $(1-\beta )^{-1}$ for $\beta \rightarrow 1$ . Further, the unstable interval of $\beta$ increases with increasing $\epsilon$ , implying that the instability is present over a larger range of polymer concentrations as domain curvature becomes important. Interestingly, the crossing of the different curves in figure 15(a) shows that, in addition to the increase in the unstable $\beta$ -interval, it is the smallest $\beta$ that become the most unstable with increasing $\epsilon$ , in terms of having the smallest $\textit{Wi}$ -thresholds – this is consistent with both the widening and lowering of the unstable region in figure 15(b).

Figure 16.

Eigenspectrum showing the presence of unstable centre mode in Dean flow of a FENE-P fluid (with $n = 40$ , $\epsilon = 0.1$ , $\alpha = 0$ , $\beta = 0.98$ , $L = 100$ and $\textit{Wi} = 60$ ). The inset shows the zoomed-in view highlighting the region near the unstable centre mode; the discrete eigenvalues obtained for $N = 400$ and $N = 450$ show good convergence.

Figure 16. Long description

Panel: A scatter plot showing the eigenspectrum of a FENE-P fluid in Dean flow. The horizontal axis represents the real part of the eigenvalue (ωr) with units ranging from 0 to 350. The vertical axis represents the imaginary part of the eigenvalue (ωi) with units ranging from -200 to 0. The plot includes data points for two different values of N, specifically N = 400 and N = 450, represented by black circles and red asterisks, respectively. A red line connects the data points, indicating the presence of unstable center modes. The inset zooms in on the region near the unstable center mode, highlighting the discrete eigenvalues obtained for N = 400 and N = 450, which show good convergence. The inset’s horizontal axis ranges from 334 to 346, and the vertical axis ranges from -2 to 0.

4.2. Role of finite extensibility on CM in Dean flow

To begin with, in figure 16, we present the elastic spectrum for the FENE-P fluid with $L = 100$ and $\beta = 0.98$ , focusing on two-dimensional disturbances. Interestingly, the CM is unstable here, even at the relatively modest $\textit{Wi} = 60$ . The effect of decreasing $L$ on the CM has been shown to be quite subtle in the case of channel flow (Khalid et al. Reference Khalid, Badoni, Dutta, Naidu, Subramanian and Shankar2025), in that it could be uniformly stabilising or destabilising, or even non-monotonic; we find a similar trend for Dean flow in the following. The effect of introducing axial variation was also investigated for $L = 100$ and we again found (data not shown) a slight destabilisation as $\alpha$ is increased up to approximately 1; a further increase in $\alpha$ leads to stabilisation, similar to the Oldroyd-B case discussed previously (figure 14 a).

Khalid et al. (Reference Khalid, Badoni, Dutta, Naidu, Subramanian and Shankar2025) first showed, for plane Poiseuille flow, that the neutral curves, for the FENE-P model, form closed loops in the $\textit{Wi}$ – $k$ plane, as opposed to the open unstable tongues seen earlier in figure 11 for the Oldroyd-B case. The closing out of the tongues, leading to loops, is a manifestation of the expected stabilising role of shear thinning, which results in a weakening of elastic effects at sufficiently high $\textit{Wi}$ . As shown in figure 17(a), for Dean flow too, the neutral curves are closed loops in the $\textit{Wi}$ – $n$ plane. Importantly, the destabilising effect of curvature (i.e. increasing $\epsilon$ ) persists for the FENE-P model as well, as evident from the neutral loops in figure 17(a) shifting downward when $\epsilon$ is increased to $O(1)$ values. We also find very good agreement between the loop for $\epsilon$ $= 0.01$ with the one obtained by Khalid et al. (Reference Khalid, Badoni, Dutta, Naidu, Subramanian and Shankar2025) for channel flow. Furthermore, as shown in figure 17(b), the destabilising effect of $L$ on the CM, originally predicted for channel flow (Page et al. Reference Page, Dubief and Kerswell2020; Khalid et al. Reference Khalid, Badoni, Dutta, Naidu, Subramanian and Shankar2025), persists even for Dean flow at finite $\epsilon$ . Indeed, an obvious signature of the finite-extensibility-induced destabilisation is the fact that we continue to get unstable islands even for $\beta = 0.98$ in figure 17(a); the flow is stable for the Oldroyd-B model for this $\beta$ .

In Appendix C, we demonstrate that the CM instability is not governed solely by the ratio $\textit{Wi}/L$ (for different $(\textit{Wi}, L)$ pairs), although the base state profiles are, as was demonstrated in § 2.3. The eigenfunctions for various field variables, associated with the CM in Dean flow (for $L = 100$ and $\epsilon = 0.1$ ), are shown in Appendix D by way of comparison with those for the HSM. While key features of the CM eigenfunctions are rather robust to changes in $\epsilon$ and $L$ , they remain qualitatively different from the corresponding HSM eigenfunctions.

Figure 17.

Centre-mode instability in Dean flow of a FENE-P fluid: (a) neutral curves at different $\epsilon$ for $L = 100$ . For $\epsilon = 0.01$ and $0.1$ , we also plot the neutral curve (continuous lines) obtained for pressure-driven flow in a rectilinear channel (Khalid et al. Reference Khalid, Badoni, Dutta, Naidu, Subramanian and Shankar2025); (b) variation of $\textit{Wi}_c$ with $\epsilon$ at different $L$ . Data for both panels correspond to $\textit{Re} = 0$ , $\alpha = 0$ and $\beta = 0.98$ .

5.1. Role of finite extensibility on axisymmetric modes

Figure 20.

Neutral curves for HSM1 in the $\textit{Wi}$ – $\alpha$ plane for Dean flow of a FENE-P fluid at $\textit{Re} = 0$ , $\beta = 0.98$ , $n = 0$ : (a) $\epsilon = 0.01$ ; (b) $\epsilon = 0.1$ ; (c) $\epsilon = 0.5$ and (d) $\epsilon = 1$ .

Figure 20. Long description

Panel A: A line graph shows the relationship between Wi and alpha for epsilon equals 0.01. The x-axis represents alpha on a logarithmic scale ranging from 10^-2 to 10^2, and the y-axis represents Wi on a logarithmic scale ranging from 10^1 to 10^5. The graph includes lines for L equals 536, 600, and 700, and a line labeled Oldroyd-B. The lines for L values form closed loops, while the Oldroyd-B line decreases and then slightly increases. Panel B: A line graph shows the relationship between Wi and alpha for epsilon equals 0.1. The x-axis represents alpha on a logarithmic scale ranging from 10^-1 to 10^1, and the y-axis represents Wi on a logarithmic scale ranging from 10^1 to 10^4. The graph includes lines for L equals 175, 180, and 200, and a line labeled Oldroyd-B. The lines for L values form closed loops, while the Oldroyd-B line decreases and then slightly increases. Panel C: A line graph shows the relationship between Wi and alpha for epsilon equals 0.5. The x-axis represents alpha on a logarithmic scale ranging from 10^0 to 10^1, and the y-axis represents Wi on a logarithmic scale ranging from 10^1 to 10^4. The graph includes lines for L equals 90 and 100, and a line labeled Oldroyd-B. The lines for L values form closed loops, while the Oldroyd-B line decreases and then slightly increases. Panel D: A line graph shows the relationship between Wi and alpha for epsilon equals 1. The x-axis represents alpha on a logarithmic scale ranging from 10^-1 to 10^2, and the y-axis represents Wi on a logarithmic scale ranging from 10^1 to 10^4. The graph includes lines for L equals 73, 75, and 100, and a line labeled Oldroyd-B. The lines for L values form closed loops, while the Oldroyd-B line decreases and then slightly increases.

In figure 20, we show the neutral curves for HSM1 (the most critical mode in the Oldroyd-B model) in the $\textit{Wi}$ – $\alpha$ plane for $L$ ranging from the Oldroyd-B limit down to $100$ . The eigenspectra corresponding to these neutral curves were presented in figure 6. Owing to the stabilising role of finite extensiblity at large $\textit{Wi}$ , the neutral stability curves for finite $L$ are again closed loops in the $\textit{Wi}$ – $\alpha$ plane for the entire range of gap width ratios examined, analogous to the neutral loops for the CM in figure 17(a). Similar closing out of the neutral curves due to shear-thinning-induced weakening of the first normal stress difference has been demonstrated using the Giesekus (Oztekin et al. Reference Oztekin, Brown and McKinley1994) and FENE-CR (McKinley et al. Reference McKinley, Oztekin, Byars and Brown1995) models for the HSM in the cone-and-plate geometry. Unlike the Oldroyd-B model, which exhibits a kink due to a transition from stationary to propagating modes, for the $L$ shown in figure 20, the transition to instability occurs via propagating modes for all $\epsilon$ shown.

For the Dean flow examined here, shear-thinning-induced stabilisation, as expected, is strongest for the smallest $\epsilon$ ( $ = 0.01$ ), owing to the high $\textit{Wi}$ – the unstable loop in figure 20(a) already is vanishingly small at $L = 536$ . In contrast, for $\epsilon = 0.1$ , the loop remains finite in extent, even for $L = 175$ . This sensitivity to $\epsilon$ is further highlighted in figure 21, where we show neutral stability loops in the $\textit{Wi}$ – $\alpha$ plane for $L = 100$ , over a range of $\epsilon$ . The size of the unstable region decreases with decreasing $\epsilon$ for the said $L$ , with the HSM being absent for $\epsilon \lesssim 0.36$ . For smaller $\epsilon$ , it is necessary to further examine modes with $n \neq 0$ to ascertain the presence of an HSM instability, which we do in the next subsection.

Figure 21.

Neutral curves for HSM1 at different $\epsilon$ showing the absence of HSM1 for $\epsilon \lesssim 0.36$ . Data for $\textit{Re} = 0$ , $n = 0$ , $\beta = 0.98$ and $L$ = 100. The inset shows an enlarged view of the region $40 \lt \textit{Wi} \lt 200$ .

Figure 22.

Eigenspectra showing HSM1 and 2 in Dean flow of a FENE-P fluid at $\textit{Re} = 0$ , $\alpha = 7$ , $\beta = 0.98$ , $ \epsilon = 0.1, n = 1 , \textit{Wi} = 400$ at different $L$ : (a) $ L = 10^6$ ; (b) $ L = 10^3$ and (c) $ L = 100$ . The continuous spectra are indicated as ‘CS’.

Figure 22. Long description

Three scatter plots depict eigenspectra showing HSM1 and HSM2 in Dean flow of a FENE-P fluid. Panel A: A scatter plot with the horizontal axis labeled ωr and the vertical axis labeled ωi. The plot shows data points for N = 180 and N = 200, with HSM1, HSM2, and continuous spectra (CS) indicated. Panel B: A scatter plot with the horizontal axis labeled ωr and the vertical axis labeled ωi. The plot shows data points for N = 180 and N = 200, with HSM1, HSM2, and continuous spectra (CS) indicated. Panel C: A scatter plot with the horizontal axis labeled ωr and the vertical axis labeled ωi. The plot shows data points for N = 180 and N = 200, with HSM1, HSM2, and continuous spectra (CS) indicated. A zoom view is included to highlight specific details.

5.2. Role of finite extensibility on non-axisymmetric modes

To begin with, we recall the discussion in § 3.2 which had identified HSM2, a second discrete mode in the elastic spectrum, other than the dominant one (HSM1) that governed the transition for the Oldroyd-B case. Herein, we examine HSM1 and HSM2 for the FENE-P case in more detail. Thus, while figure 7(b) was for $L = 200$ , figure 22 examines the non-axisymmetric spectra (again for $n = 1$ ) for $L$ decreasing from $10^6$ (the Oldroyd-B limit for all practical purposes) down to 100. Decreasing $L$ has a stabilising effect on HSM1, as one might anticipate from the role of shear thinning on the hoop-stress instability in our earlier discussion – HSM1 is, in fact, a stable mode in figure 22(c) for $L= 100$ . Rather intriguingly, decreasing $L$ has an initial destabilising effect on HSM2, although an ultimately stabilising one. This is evident from figure 23, which shows the growth rate ( $\omega _i$ ) variation with $L$ for both HSM1 and HSM2 (for $n = 2$ ), with figure 23(b) presenting a magnified view of the interval $50 \lt L \lt 10^3$ . The non-monotonic variation for HSM2 leads to HSM1 being the critical mode for $L \gt 200$ and HSM2 being critical for $L \lt 200$ ; for $L = 100$ , HSM2 is the only unstable mode in the spectrum. An identical scenario prevails for $n = 1$ (not shown).

Figure 23.

Variation of $\omega _i$ with $L$ for the hoop-stress modes (HSM1 and HSM2) in Dean flow of a FENE-P fluid for $\epsilon = 0.1$ , $\alpha = 7$ , $n = 2$ , $\textit{Re} = 0$ , $\textit{Wi} = 130$ and $\beta = 0.98$ : (a) for $10 \lt L \lt 10^6$ and (b) zoomed-in view of panel (a) highlighting the cross-over from HSM1 to HSM2.

Figure 24.

Structure of the Dean flow hoop-stress mode (HSM2) for a FENE-P fluid, characterised via contours of the azimuthal velocity perturbation, $|\tilde {v}_\theta |$ , in the $r{-}z$ plane; $\beta = 0.98$ , $L = 100$ , $\textit{Wi} = 400$ , $n = 1$ , $\alpha = 7$ . In contrast to the centre-mode eigenfunction, which is localised near the velocity maximum, the HSM2 mode is localised near the inner wall.

The aforementioned destabilising influence of finite extensibility on HSM2 is reminiscent of its effect on the centre mode in channel flow (Buza et al. Reference Buza, Page and Kerswell2022b ; Khalid et al. Reference Khalid, Badoni, Dutta, Naidu, Subramanian and Shankar2025), and in Dean flow (as discussed in § 4.2). However, HSM2 is not a centre mode and has more in common with HSM1, as evident from the $r{-}z$ plane eigenfunction contours depicted in figures 24 and 25, and that are based on the amplitude of the azimuthal velocity perturbation; the choice of plane is on account of the primarily axisymmetric character of the HSMs. The only difference between the two HSMs is the region of localisation – the HSM2 (HSM1) eigenfunctions are localised near the inner (outer) radius. It is worth contrasting the structures in figures 24 and 25 with the CM arrowhead structures shown earlier in figure 12. This contrast is further demonstrated in Appendix D, which compares eigenfunctions of the different field variables for HSM1, HSM2 and CM.

Figure 25.

Structure of the Dean flow hoop-stress mode (HSM1) for a FENE-P fluid, characterised via contours of the azimuthal velocity perturbation, $|\tilde {v}_\theta |$ , in the $r{-}z$ plane. In contrast to the centre-mode eigenfunction, which is localised near the velocity maximum, the HSM1 mode is localised near the outer wall.

Figure 26.

Neutral curves for CM and non-axisymmetric HSM2 for the (a) Oldroyd-B and (b) FENE-P ( $L = 100$ ) models. Data for $n = 2$ , $\textit{Re} = 0$ , $\epsilon = 1$ , $\beta = 0.98$ . The results for axisymmetric HSM1 are shown for reference.

Figure 26. Long description

Panel A: A line graph titled Oldroyd-B shows the relationship between Wi and alpha. The x-axis is labeled alpha and ranges from 10^-1 to 10^1. The y-axis is labeled Wi and ranges from 10^0 to 10^3. The graph includes four lines representing different models: HSM1 (n = 0) in black, HSM1 (n = 2) in red, CM (n = 2) in blue, and HSM2 (n = 2) in green. The lines show various trends and curves. Panel B: A line graph titled L = 100 shows the relationship between Wi and alpha. The x-axis is labeled alpha and ranges from 10^-2 to 10^1. The y-axis is labeled Wi and ranges from 10^0 to 10^3. The graph includes three lines representing different models: HSM1 (n = 0) in black, HSM2 (n = 2) in green, and CM (n = 2) in blue. The lines show different trends and curves, with an oval highlighting a specific area of interest.

For the FENE-P case, the destabilising influence of finite extensibility on both CM and HSM2 might lead one to anticipate that these two modes might be connected in some regions of the parameter space. Figure 26 shows that this is not the case. This figure shows the neutral curves for HSM2 and CM, in the $\textit{Wi}$ – $\alpha$ plane, for both Oldroyd-B and FENE-P ( $L = 100$ ) fluids, for $\epsilon = 1$ . For the Oldroyd-B fluid, while the CM is stable for $\alpha = 0$ (when $\beta = 0.98$ ), it is destabilised by a non-zero $\alpha$ . In addition, there are neutral curves related to HSM1 and HSM2 (the HSM1 neutral curve for $n = 2$ is not shown). The critical $\textit{Wi}$ for the $n = 2$ mode in the FENE-P ( $L = 100$ ) model is significantly higher than that for the Oldroyd-B case, due to the stabilising influence of finite extensibility. The unstable region for HSM2 (with $n = 2$ ) is completely disconnected from the CM, both for Oldroyd-B and FENE-P fluids. As $\alpha$ is increased for a fixed $\textit{Wi}$ (in figure 26), we find that the CM becomes stable and eventually merges into the CS, a trend that is reminiscent of the CM in rectilinear channel and pipe flows (Chaudhary et al. Reference Chaudhary, Garg, Subramanian and Shankar2021; Khalid et al. Reference Khalid, Chaudhary, Garg, Shankar and Subramanian2021a ). Similarly, when we tracked HSM2 for a $\textit{Wi}$ from the unstable region, the $\omega _r$ for this mode is such that it is beyond the CS range and, hence, this mode circumvents the CS and eventually ends up below it, without crossing the same. This shows that the two modes (HSM2 and CM) do not continue to each other, as $\alpha$ is varied for a fixed $\textit{Wi}$ . While the discussion here pertains to HSM2 and CM, it is relevant to emphasise that for the larger gap-width ratio of $\epsilon = 1$ , HSM1 (with $n = 0$ ) is the most critical for both Oldroyd-B and FENE-P models.

Figure 27.

Neutral curves for HSM (1 and 2) at $\textit{Re} = 0$ , $\beta = 0.98$ and $L = 100$ : (a) $\epsilon = 0.1$ ; (b) $\epsilon = 0.2$ ; (c) $\epsilon = 0.5$ and (d) $\epsilon = 1$ . While the $n \neq 0$ modes (HSM2) are the most unstable for $\epsilon = 0.1$ and $0.2$ , and remain unstable even for $\epsilon = 0.5$ and $1$ , the $n = 0$ mode (HSM1) becomes critical at $\epsilon = 0.5$ and $1$ .

As seen in figure18, for the Oldroyd-B case, the axisymmetric HSM1 is the most unstable for any $\epsilon$ . However, the trend is very different for the FENE-P case, being sensitive to the choice of $\epsilon$ and $L$ . For $L \sim 100$ , the nature of the most unstable mode (viz., axisymmetric versus non-axisymmetric HSM2) depends on $\epsilon$ , as shown in figure 27. In the narrow-gap regime ( $\epsilon = 0.01 {-} 0.2$ ), the axisymmetric modes (both HSM1 and HSM2) remain stable for $L = 100$ , while a non-axisymmetric HSM2 becomes critical. Specifically, at $\epsilon = 0.1$ , $n = 2$ is the critical mode (HSM2), whereas at $\epsilon = 0.2$ , $n = 1$ emerges as the critical mode (HSM2). In contrast, for the larger values $\epsilon = 0.5$ and $1$ , the axisymmetric mode (HSM1) becomes critical. Notably, the range of $\epsilon$ between $0.1$ and $0.3$ corresponds to gap widths used in the Taylor–Dean experiments of Joo & Shaqfeh (Reference Joo and Shaqfeh1994), although the serpentine channel experiments of Groisman & Steinberg (Reference Groisman and Steinberg2001, Reference Groisman and Steinberg2004) used a higher $\epsilon \,(= 1)$ ; we return to the issue of the dominant mode in experiments in § 6.

For the Oldroyd-B model, $\textit{Wi}_c$ was seen to scale as $\epsilon ^{-1/2}$ for $\epsilon \ll 1$ . In contrast, and as already demonstrated in figure 20 and also in figure 27, one expects the HSM instability to be fully suppressed in the limit $\epsilon \rightarrow 0$ for any finite $L$ . This is illustrated in figure 28(a), where $\textit{Wi}_c$ is seen to diverge faster than $\epsilon ^{-1/2}$ for all finite $L$ . The data points in figure 28(a) for $\textit{Wi}_c$ (at different $\epsilon$ or $L$ ) could correspond to either HSM1 ( $n = 0$ ) or HSM1 ( $n \neq 0$ ) or HSM2 ( $n \neq 0$ ). Importantly, as $L$ is decreased to realistic values of $O(100)$ , HSM1 is fully suppressed, with a suggestion of the threshold diverging at a finite $\epsilon$ . This is in qualitative agreement with what one obtains at a scaling level by incorporating shear thinning in the Pakdel–McKinley criterion (see figure 15 of Castillo Sánchez et al. Reference Castillo Sánchez, Jovanović, Kumar, Morozov, Shankar, Subramanian and Wilson2022). Figure 28(b) shows the critical azimuthal wavenumber $n_c$ as a function of $\epsilon$ for different $L$ ; the notation used is similar to figure 28(a), in that open symbols pertain to HSM1, while filled symbols to HSM2. While in the Oldroyd-B limit, $n = 0$ remains critical for all $\epsilon$ , the critical azimuthal wavenumber for finite $L$ changes with $\epsilon$ and there does not appear to be a simple trend. Figure 28(c) shows the corresponding critical axial wavenumber $\alpha _c$ . The abrupt variation of $\alpha _c$ with $\epsilon$ in this figure is on account of the transition from HSM1 to HSM2, and also due to changes in $n_c$ .

Figure 28.

(a) $\textit{Wi}_c$ versus $\epsilon$ for HSM (1 and 2) in Dean flow at different $L$ illustrating the stabilising influence of finite extensibility. Here, filled symbols of a given shape denote HSM2 and the corresponding open ones denote HSM1; different symbol shapes corresponds to different $L$ . Corresponding (b) $n_c$ versus $\epsilon$ and (c) $n_c$ versus $\alpha _c$ . Data for $\textit{Re} = 0$ and $\beta = 0.98$ .

Figure 28. Long description

Panel A: A scatter plot showing Wi_c versus ε for different channel lengths (L) and models (HSM1 and HSM2). The x-axis represents ε on a logarithmic scale ranging from 10^-3 to 10^0, and the y-axis represents Wi_c on a logarithmic scale ranging from 10^0 to 10^2. Different symbols and colors denote different channel lengths and models, with filled symbols for HSM2 and open symbols for HSM1. Panel B: Multiple scatter plots showing n_c versus ε for different channel lengths and models. The x-axis represents ε ranging from 0 to 1, and the y-axis represents n_c ranging from 0 to 5. Different symbols and colors denote different channel lengths and models. Panel C: Multiple scatter plots showing α_c versus ε for different channel lengths and models. The x-axis represents ε ranging from 0 to 1, and the y-axis represents α_c ranging from 4 to 10. Different symbols and colors denote different channel lengths and models.