Inertial enhancement of the polymer diffusive instability

Abstract

Beneitez et al. (Phys. Rev. Fluids, vol. 8, 2023, L101901) have recently discovered a new linear ‘polymer diffusive instability’ (PDI) in inertialess rectilinear viscoelastic shear flow using the finitely extensible nonlinear elastic constitutive model of Peterlin (FENE-P) when polymer stress diffusion is present. Here, we examine the impact of inertia on the PDI for both plane Couette and plane Poiseuille flows under varying Weissenberg number ${W}$, polymer stress diffusivity $\varepsilon$, solvent-to-total viscosity ratio $\beta$ and Reynolds number ${Re}$, considering the FENE-P and simpler Oldroyd-B constitutive relations. Both the prevalence of the instability in parameter space and the associated growth rates are found to significantly increase with ${Re}$. For instance, as $Re$ increases with $\beta$ fixed, the instability emerges at progressively lower values of $W$ and $\varepsilon$ than in the inertialess limit, and the associated growth rates increase linearly with $Re$ when all other parameters are fixed. For finite $Re$, it is also demonstrated that the Schmidt number $Sc=1/(\varepsilon Re)$ collapses curves of neutral stability obtained across various $Re$ and $\varepsilon$. The observed strengthening of PDI with inertia and the fact that stress diffusion is always present in time-stepping algorithms, either implicitly as part of the scheme or explicitly as a stabilizer, implies that the instability is likely operative in computational work using the popular Oldroyd-B and FENE-P constitutive models. The fundamental question now is whether PDI is physical and observable in experiments, or is instead an artifact of the constitutive models that must be suppressed.

1. Introduction

The addition of polymers to a Newtonian solvent can induce dramatically different flow behaviours compared with those observed in the Newtonian fluid alone (Datta et al. 2022; Sánchez et al. 2022). In industrial processes, for instance, viscous polymer melts are susceptible to instabilities which constrain the maximum extrusion rate (Petrie & Denn 1976), while polymer additives are used in oil pipelines to reduce turbulent wall drag (Virk 1975). Two particularly important viscoelastic phenomena are the existence of ‘elastic turbulence’ (ET), a chaotic flow state sustained in the absence of inertia (Groisman & Steinberg 2000; Steinberg 2021), and ‘elasto-inertial turbulence’ (EIT), an inherently two-dimensional state arising when both inertia and elasticity are present (Samanta et al. 2013; Sid, Terrapon & Dubief 2018; Choueiri et al. 2021). While the initial pathway to ET in curvilinear geometries is understood (Larson, Shaqfeh & Muller 1990; Pakdel & McKinley 1996; Shaqfeh 1996; Datta et al. 2022), relatively little is known about what happens in rectilinear viscoelastic flows.

Initial progress in characterizing ET in rectilinear situations arose through consideration of Kolmogorov flow over a two-torus, where Boffetta et al. (2005) found a linear instability leading to ET (Berti & Boffetta 2010). Garg et al. (2018) subsequently discovered a centre-mode instability in viscoelastic pipe flow at finite Reynolds number $Re$, which was later also identified in plane Poiseuille flow (PPF, hereafter referred to as channel flow) (Khalid et al. 2021a) but notably not in plane Couette flow (PCF). Interestingly, this instability could only be traced down to $Re=0$ in channel flow (Khalid, Shankar & Subramanian 2021b; Buza, Page & Kerswell 2022b). The finite-amplitude state resulting from this instability is an ‘arrowhead’ travelling wave (Page, Dubief & Kerswell 2020; Buza et al. 2022a; Morozov 2022) which has been observed in channel flow EIT (Dubief et al. 2022) and, in retrospect, ET in two-dimensional Kolmogorov flow (Berti & Boffetta 2010). In channel flow, efforts have begun to establish a dynamical link between the arrowhead solution and both ET (Lellep, Linkmann & Morozov 2023) and EIT (Beneitez et al. 2023a). In the latter case in two dimensions, there does not appear to be a simple dynamical pathway between these arrowhead solutions, where the dynamics is concentrated near the midplane, and EIT (Beneitez et al. 2023a), which seems more dependent on a near-wall mechanism (Shekar et al. 2019, 2021; Dubief et al. 2022).

The very recent discovery of a new wall-mode ‘polymer diffusive instability’  (PDI) in plane Couette flow at $Re=0$ (Beneitez, Page & Kerswell 2023b), however, has added another intriguing possibility for the origin of ET. This instability is dependent on the existence of small but non-vanishing polymer stress diffusion which is invariably present in any time-stepping algorithm, whether added explicitly to stabilize a numerical scheme like a spectral method (see e.g. Dubief, Terrapon & Hof 2023) or arising implicitly such as through a finite difference formulation (see e.g. Zhang et al. 2015; Pimenta & Alves 2017). The PDI wall mode is primarily confined to a boundary layer of thickness $\sqrt {\varepsilon }$, where $\varepsilon$ is the (small) diffusion coefficient, travelling at the wall speed with a streamwise wavelength of the order of the boundary layer thickness. The instability is robust to the choice of boundary conditions applied to the polymer conformation equation, and has growth rates which remain $O(1)$ as $\varepsilon \rightarrow 0$. Direct numerical simulations (DNS) have demonstrated that PDI can lead to a sustained three-dimensional turbulent state, thus providing a potential mechanism for the origin of an ET-like state in the popular finitely extensible nonlinear elastic constitutive model of Peterlin (FENE-P) (Beneitez et al. 2023b).

While PDI has the potential to be a viscoelastic instability of significant importance, there is an important caveat: the instability emerges at small length scales which can, depending on the parameters, approach the order of the polymer gyration radius (Beneitez et al. 2023b), violating the continuum approximation. There is thus a question of whether the instability is physical or actually an artifact of the Oldroyd-B and FENE-P models. Either possibility has important implications: if the instability is a physical phenomenon then it provides a pathway to ET and EIT, albeit one which will likely be challenging to establish experimentally due to the small length scales involved; or, it is an artificial feature of the popular FENE-P model, which can compromise its predictions. It thus appears important to now establish the prevalence of PDI across a much wider region of parameter space than was considered in the initial study of Beneitez et al. (2023b). Therefore, we here map out the regions where PDI is operative at finite $Re$, considering both plane Couette flow and the more experimentally relevant channel flow scenarios. Both the prevalence of PDI and the associated growth rates are found to significantly increase at finite $Re$ and are relatively insensitive to the bulk flow geometry. Therefore, PDI is a candidate to trigger both ET and EIT in simulations using the FENE-P model.

2. Formulation

We consider the following dimensionless equations governing the flow of an incompressible viscoelastic fluid:

(2.1a,b)\begin{gather} Re\left(\frac{\partial\boldsymbol{u}}{\partial t}+\left(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\right)\boldsymbol{u}\right)={-}\boldsymbol{\nabla} p+\beta\nabla^{2}\boldsymbol{u}+\left(1-\beta\right)\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\tau}, \quad \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{gather}

(2.2)\begin{gather} \frac{\partial\boldsymbol{\boldsymbol{c}}}{\partial t}+\left(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\right) \boldsymbol{c}+\boldsymbol{\tau}=\boldsymbol{c}\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{u} +\left(\boldsymbol{\nabla}\boldsymbol{u}\right)^{T} \boldsymbol{\cdot}\boldsymbol{\boldsymbol{c}} +\varepsilon\nabla^{2}\boldsymbol{\boldsymbol{c}}, \end{gather}

where $\boldsymbol {u}=(u,v,w)$ and $p$ denote the velocity and pressure fields, respectively, and $\boldsymbol {\tau }$ denotes the polymeric contribution to the stress tensor. Following Beneitez et al. (2023b), the equations have been non-dimensionalized using the channel half-width $H$ and a characteristic flow speed $U_{0}$, taken to be the wall speed in plane Couette flow or the centreline velocity of the base channel flow. In (2.1a), the Reynolds number $Re:=U_{0} H/\nu _T$ describes the ratio of inertial to viscous forces (with $\nu _T = \nu _S + \nu _P$ denoting the total kinematic viscosity comprised of solvent and polymer components), and $\beta := \nu _S/\nu _T$ denotes the ratio of solvent to total viscosity. The polymeric stress tensor $\boldsymbol {\tau }$ may be described in terms of the polymer orientation through the conformation tensor $\boldsymbol {c}$ as in (2.2). We emphasize that the inclusion of a polymer stress diffusion term $\varepsilon \nabla ^{2}\boldsymbol {\boldsymbol {c}}$ (associated with diffusivity $\varepsilon := (Re\,Sc)^{-1}$, where $Sc$ denotes the Schmidt number) is the crucial ingredient for the PDI identified by Beneitez et al. (2023b). In the inertialess limit, $\varepsilon =Sc^{-1}$ when the governing equations are non-dimensionalized using viscous scales.

To close equations (2.1)–(2.2), the FENE-P constitutive model is used to relate $\boldsymbol {\tau }$ and $\boldsymbol {c}$:

(2.3a,b)\begin{equation} \boldsymbol{\tau}:=\frac{f\left(\mathrm{tr}\,\boldsymbol{c}\right)\boldsymbol{c}-\boldsymbol{\mathsf{I}}}{W},\quad f\left(s\right):=\left(1-\frac{s-3}{L^{2}}\right)^{{-}1}, \end{equation}

where $\boldsymbol{\mathsf{I}}$ is the identity matrix, $L$ denotes the maximum extensibility of the polymer chains and $W:=U_{0} \lambda /H$, the Weissenberg number, describes the ratio of time scales for polymer relaxation ($\lambda$) to the flow. In the limit $L\rightarrow \infty$, the simpler Oldroyd-B model is obtained. Inspection of (2.1)–(2.3) reveals five parameters of interest governing the flow dynamics: $Re$, $W$, $\beta$, $\varepsilon$, $L$.

We analyse the linear stability of (2.1)–(2.3) by perturbing them about their base state: $\boldsymbol {u}=\boldsymbol {U}+\boldsymbol {u}^{*}$, $p=P+p^{*}$, $\boldsymbol {\tau }=\boldsymbol {T}+\boldsymbol {\tau }^{*}$, $\boldsymbol {c}=\boldsymbol {C}+\boldsymbol {c}^*$. Our coordinate system ($x$, $y$, $z$) is aligned with the streamwise, wall-normal and spanwise directions of the channel, respectively. The base state $(U(\kern0.7pt y),C_{xx}(\kern0.7pt y),C_{yy}(\kern0.7pt y),C_{zz}(\kern0.7pt y),C_{xy}(\kern0.7pt y))$ satisfies

(2.4a)$$\begin{gather} -\partial_{x}P+\beta U''+\left(1-\beta\right)\partial_{y}T_{xy}=0, \end{gather}$$

(2.4b)$$\begin{gather}f\left(\mathrm{tr}\boldsymbol{C}\right)\boldsymbol{C}-\varepsilon W\boldsymbol{C}''-\begin{pmatrix}2WU'C_{xy}+1 & WU'C_{yy} & 0\\ WU'C_{yy} & 1 & 0\\ 0 & 0 & 1 \end{pmatrix}=0, \end{gather}$$

where primes indicate derivatives in the wall-normal ($y$) direction. We use pressure gradients $\partial _{x} P = 0$ and $\partial _{x} P = -2$ for plane Couette and channel flow, respectively. While we here compute the base flow accounting for the presence of a finite diffusivity $\varepsilon$ in (2.4b), it is worth noting that the presence of diffusion in the base flow is not required to induce PDI; the instability still arises if $\varepsilon =0$ in (2.4b) as the inclusion of $\varepsilon \neq 0$ does not significantly change the base state except in the limit of large $\varepsilon = O(1)$.

Normal mode solutions of the perturbed flow are sought using the ansatz $\phi ^{*}(x,y,t)=\tilde {\phi }(\kern0.7pt y){\rm e}^{{\rm i}k(x-ct)}$, where real-valued $k$ denotes the streamwise wavenumber and $c=c_{r}+{\rm i}c_{i}$ is a complex wave speed, with instability arising if $c_i> 0$. The perturbed state is governed by the following system of seven equations for $(\tilde {u},\tilde {v},\tilde {p},\tilde {c}_{xx},\tilde {c}_{yy},\tilde {c}_{zz},\tilde {c}_{xy})$:

(2.5a)$$\begin{gather} {\rm i}k\tilde{u}+\tilde{v}' = 0, \end{gather}$$

(2.5b)$$\begin{gather}Re\left(-{\rm i}kc\tilde{u}+\tilde{v}U'+{\rm i}kU\tilde{u}\right)+{\rm i}k\tilde{p}-\beta\left({-}k^{2}\tilde{u}+\tilde{u}''\right)-\left(1-\beta\right)\left({\rm i}k\tilde{\tau}_{xx}+\tilde{\tau}'_{xy}\right)=0, \end{gather}$$

(2.5c)$$\begin{gather}Re\left(-{\rm i}kc\tilde{v}+{\rm i}kU\tilde{v}\right)+ \tilde{p}'-\beta\left({-}k^{2}\tilde{v}+\tilde{v}''\right)-\left(1-\beta\right)\left({\rm i}k\tilde{\tau}_{xy}+\tilde{\tau}'_{yy}\right)=0, \end{gather}$$

(2.5d)$$\begin{gather}\left[\varepsilon k^{2}+{\rm i}k\left(U-c\right)\right]\tilde{c}_{xx}+\tilde{v}C'_{xx}+\tilde{\tau}_{xx}-\varepsilon\tilde{c}''_{xx}-2\left({\rm i}kC_{xx}\tilde{u}+C_{xy}\tilde{u}'+\tilde{c}_{xy}U'\right)=0, \end{gather}$$

(2.5e)$$\begin{gather}\left[\varepsilon k^{2}+{\rm i}k\left(U-c\right)\right]\tilde{c}_{yy}+\tilde{v}C'_{yy}+\tilde{\tau}_{yy}-\varepsilon\tilde{c}''_{yy} -2\left({\rm i}kC_{xy}\tilde{v}+C_{yy}\tilde{v}'\right)=0, \end{gather}$$

(2.5f)$$\begin{gather}\left[\varepsilon k^{2}+{\rm i}k\left(U-c\right)\right]\tilde{c}_{zz}+\tilde{v}C'_{zz}+\tilde{\tau}_{zz}-\varepsilon\tilde{c}''_{zz} =0, \end{gather}$$

(2.5g)$$\begin{gather}\left[\varepsilon k^{2}+{\rm i}k\left(U-c\right)\right]\tilde{c}_{xy}+\tilde{v}C'_{xy}+\tilde{\tau}_{xy}-\varepsilon\tilde{c}''_{xy} -{\rm i}kC_{xx}\tilde{v}-U'\tilde{c}_{yy}-C_{yy}\tilde{u}'=0. \end{gather}$$

By expanding the variables in (2.5) in terms of a basis of Chebyshev polynomials, an eigenvalue problem is obtained which may be solved using standard Python libraries (Beneitez et al. 2023b). A basis of $N=300$ polynomials was used here to target the eigenvalues associated with PDI via an inverse iterations algorithm, which yielded sufficient convergence. Representative eigenvalue spectra are shown in figure 1 for plane Couette (PCF) and channel (PPF) geometries, illustrating the results using bases of both $N=300$ and $N=400$ Chebyshev polynomials to demonstrate convergence. The unstable PDI mode emerges with characteristic wave speeds in the vicinity of $|c_{r}|\approx 1$ and $c_r \approx 0$ for PCF and PPF geometries, respectively.

Figure 1. Eigenvalue spectra, plotted in terms of the real $c_r$ and imaginary $c_i$ components of the complex wave speed, obtained by solving the system (2.5a)–(2.5g) using the Oldroyd-B constitutive model at $Re=1000$, $\varepsilon =10^{-5}$ and $\beta =0.9$ for (a) plane Couette flow ($W=44.7$, $k=37$) and (b) channel flow ($W=22.1$, $k=52$). These parameters are chosen to roughly lie on the corresponding neutral curves plotted in figure 2(a). Results using bases of $N=300$ and $400$ Chebyshev polynomials are depicted by blue and red markers, respectively, to demonstrate convergence of the discrete unstable PDI modes highlighted in the zoomed inset above each panel.

Following Beneitez et al. (2023b), we choose boundary conditions that set $\varepsilon =0$ at the walls in the governing equations to match the typical configuration used in DNS. While $Sc\sim O(10^{6})$ is characteristic of polymer diffusion in a typical solvent like water (yielding $\varepsilon = 1/(Re\,Sc) \approx 10^{-9}-10^{-6}$ depending on $Re$), in DNS much larger values of $\varepsilon \approx O(10^{-3})$ are often employed to achieve numerical stability (see e.g. § 2.2 of Dubief et al. 2023). Thus, it is commonplace to set $\varepsilon =0$ at the walls in order to recover the idealized limit $\varepsilon \to 0$ (Sureshkumar, Beris & Handler 1997; Samanta et al. 2013; Dubief et al. 2022).

3. Results of linear stability analysis

The results of our linear stability analysis are presented for both the Oldroyd-B (§ 3.1) and FENE-P (§ 3.2) constitutive relations, delineating the trajectory of neutral curves (characterized by growth rate $\sigma := kc_i=0$) associated with PDI in the four-dimensional parameter space spanned by $(Re, W, \beta,\varepsilon )$. For a given point in parameter space, we perform a sweep over wavenumbers $k$ to ensure that we are tracking the most unstable mode. The results for plane-Couette (PCF) and channel (PPF) geometries are overlaid on the same axes using blue and red colouring, respectively. To promote a collapse of the curves for both geometries, the Weissenberg ($W$) axis is scaled by the respective wall shear rates: $U_{{wall}}'=\{ 1\ (\textrm {PCF}),\ 2\ (\textrm {PPF})\}$. As described in § 1, the PDI is a wall mode confined to a boundary layer of thickness $O(\sqrt {\varepsilon })$ and so its behaviour will primarily be influenced by the wall shear. However, we will demonstrate that this collapse of geometries fails in certain scenarios (such as for large $\varepsilon$), when the boundary layer grows sufficiently large to be influenced by the non-uniform shear profile away from the wall. In the case of finite $Re$, we also demonstrate that curves associated with variable $Re$ and $\varepsilon$ may be collapsed based on $Sc = 1/(\varepsilon Re)$.

3.1. The Oldroyd-B case

The trajectory of neutral curves in the $Re$–$W$–$\varepsilon$ volume are presented in figure 2 for various $\beta \in [0.7,0.98]$. Figure 2(a) considers the $Re$–$W$ plane at a fixed $\varepsilon = 10^{-5}$, while figure 2(b) considers the $\varepsilon$–$W$ plane at three fixed $Re\in \{ 0,1000,5000\}$. Regions of stability and instability are found to the left and right of each curve, respectively. For a given $\beta$, increasing $Re>0$ is found to promote instability at progressively lower values of both $W$ (figure 2a) and $\varepsilon$ (figure 2b). Notably, in the inertialess limit ($Re=0$), Beneitez et al. (2023b) found that the neutral curves tracked a roughly fixed $W \sim \beta / (1-\beta )$ for $\varepsilon \lesssim 10^{-2}$ (see their figure 2b) leading to the disappearance of PDI at finite $W$ in the limit $\beta \rightarrow 1$. Here, for $Re>0$, our figure 2(b) demonstrates that while for small $\varepsilon =O(10^{-7})$ the neutral curves do indeed still track a constant $W$, they then begin to significantly deviate to lower $W$ beyond $\varepsilon \gtrsim 10^{-6}$, with the $Re=5000$ case deviating at lower $\varepsilon$ than the $Re=1000$ case. Inertial effects thus play a significant role in promoting a greater prevalence of PDI in the parameter space. In figure 2(d), we demonstrate that the $Re=\{ 1000,5000\}$ curves plotted in figure 2(b) may be collapsed for both geometries based on the Schmidt number $Sc = 1/(\varepsilon Re)$.

Figure 2. Curves of neutral stability for plane Couette (‘PCF’, blue) and channel (‘PPF’, red) geometries, using the Oldroyd-B constitutive relation for five values of $\beta \in [0.7,0.98]$. (a) The $Re$–$W$ plane for fixed $\varepsilon = 10^{-5}$, noting that the PCF and PPF curves are virtually indistinguishable. (b) The $\varepsilon$–$W$ plane at three fixed $Re=\{ 0,1000,5000\}$. In panels (a,b), the $W$ axis is scaled by the wall shear rate: $U_{{wall}}'=\{ 1\ (\textrm {PCF}),\ 2\ (\textrm {PPF})\}$. (c) The streamwise wavenumber $k$ of the PDI, as a function of $\varepsilon$, along each of the neutral curves in panel (b). (d) A collapse of the $Re=\{ 1000,5000\}$ curves for both geometries from panel (b) based on the inverse Schmidt number $1/Sc = \varepsilon Re$. By plotting $1/Sc$ rather than $Sc$, small $\varepsilon$ occurs at the bottom of panel (d) facilitating an easier comparison with panel (b). (e,f) Colourmaps of the trace of the polymer conformation tensor $\mathrm {tr}(\boldsymbol {c})$ (red and blue denote positive and negative values, respectively), with contours of the stream function superimposed, for PPF eigenfunctions in the upper half-channel with $\beta =0.9$, at locations indicated by the square and triangular markers in panels (b,c). One wavelength $\lambda = 2{\rm \pi} / k$ of each eigenfunction is shown.

The streamwise wavenumbers $k$ associated with the neutral curves in figure 2(b) are plotted in figure 2(c) as a function of $\varepsilon$. At $Re=0$, $k$ follows the $1/\sqrt {\varepsilon }$ scaling reported by Beneitez et al. (2023b) for all $\varepsilon$, whereas for $Re>0$, $k$ deviates significantly from this scaling to plateau to a roughly constant value for $\varepsilon \gtrsim 10^{-5}$, with the unstable modes at higher $Re$ being more tightly confined to the wall (as indicated by a larger $k$). This deviation from the $1/\sqrt {\varepsilon }$ scaling for $Re>0$ corresponds to the previously noted deviation of the neutral curves away from a constant $W$ in figure 2(b), where the curves begin turning to lower $W$ in the vicinity of $\varepsilon \approx 10^{-6}-10^{-5}$. The behaviour of $k$ in figure 2(c) also explains the mismatch in the collapse of the PCF and PPF curves in figure 2(b) at higher $\varepsilon$ for $Re=0$ (see solid lines). Specifically, at higher $\varepsilon$, $k$ becomes sufficiently small at $Re=0$ such that the instability is no longer strongly confined to the wall (see square eigenfunction, figure 2e) and so the local wall shear, $U'_{wall}$, does not accurately describe the non-uniform shear profile influencing the instability. Conversely, the curves for $Re=\{ 1000,5000\}$ do remain collapsed at high $\varepsilon$, as the wavenumbers $k$ in figure 2(c) plateau to sufficiently large values such that the instability remains confined to the wall (see triangle eigenfunction, figure 2f). In figure 2(c), it is also worth noting that the prefactor of the PPF scaling (1/5) is roughly twice that of the PCF scaling (1/8), thus indicating that PDI is more tightly confined to the wall in the channel geometry.

3.2. The FENE-P case

We now consider how a finite polymer extensibility $L$ modifies the behaviour of the instability as compared with the Oldroyd-B case ($L\rightarrow \infty$) presented in § 3.1. Focusing first on $L=200$, figure 3 illustrates the much richer behaviour of the neutral curves for the FENE-P case, presented in an analogous manner to figures 2(a–d). In the $Re$–$W$ plane (figure 3a), a finite $L$ introduces two notable differences compared with Oldroyd-B. First, the neutral curves for a given $\beta$ now have a left-hand and right-hand branch and so the range of instability is bounded by an upper value of $W$. Second, there is now a critical value of $\beta$ ($\approx 0.865$ for both geometries) above which the neutral curves no longer intersect the $Re=0$ axis, at fixed $\varepsilon =10^{-5}$. Therefore, inertial effects are once again demonstrated to promote PDI, generating instability at finite $Re$ for ultradilute polymer solutions ($\beta \rightarrow 1$) that would otherwise remain stable at $Re=0$.

Figure 3. Curves of neutral stability using the FENE-P constitutive relation with a fixed extensibility $L=200$, presented analogously to the Oldroyd-B curves in figure 2 for various $\beta$. Curves are shown in (a) the $Re$–$W$ plane with a fixed $\varepsilon =10^{-5}$, (b) the $\varepsilon$–$W$ plane for fixed $Re=0$ and (c) the $1/Sc = \varepsilon Re$ vs $W$ plane, which collapses curves obtained at $Re=1000$ and $5000$ for both geometries. By plotting $1/Sc$ rather than $Sc$, small $\varepsilon$ occurs at the bottom of panel (c), facilitating an easier comparison with panel (b). Panel (d) illustrates the dependence of the optimal streamwise wavenumber $k$ on $\varepsilon$ for the PPF curves in panel (c) at $Re=1000$. Left-hand and right-hand branches of the curves in panel (c) are distinguished in panel (d) using solid and dashed lines, respectively. Comparison with the Oldroyd-B curves in figure 2(c) reveals that the left-hand branches behave similarly to the single Oldroyd-B branch, deviating from the $1/\sqrt {\varepsilon }$ scaling at large $\varepsilon$, while the right-hand branches retain this scaling to the highest $\varepsilon$ considered.

Neutral curves in the $\varepsilon$–$W$ plane for the inertialess case $Re=0$ are shown in figure 3(b), exhibiting notable differences to the Oldroyd-B $Re=0$ curves presented in figure 2(b). While for Oldroyd-B the plane Couette and channel curves adopt roughly the same shape, the FENE-P curves display dramatically different behaviours at large $\varepsilon$. Specifically, the plane Couette curves have an inverted ‘U’-shape, highlighting that the instability ceases to exist at large $\varepsilon$ for certain $\beta$ (e.g. see the $\beta =0.86$ curve which reaches a maximum value at $\varepsilon \approx 6\times 10^{-3}$). Conversely, the channel curves are roughly ‘U’-shaped, and the range of instability only increases with increasing $\varepsilon$.

At finite $Re$, as for the Oldroyd-B curves in figure 2(d), the FENE-P curves shown in figure 3(c) at $Re=\{ 1000,5000\}$ are again found to collapse for both geometries based on the Schmidt number $Sc=1/(\varepsilon Re)$. Comparing figures 2(d) and 3(c) reveals two main differences in the structure of the FENE-P and Oldroyd-B curves. First, for $\beta \approx 0.925$, in both geometries a pinch-off phenomenon occurs for the FENE-P case where the neutral curves form a ‘bubble’ of instability within the $\varepsilon -W$ plane in an otherwise stable region. Second, in the limit of $\varepsilon \rightarrow 0$ ($Sc \rightarrow \infty$), the neutral curves behave differently as compared with the $Re=0$ case reported by Beneitez et al. (2023b). In particular, there now appears to be a critical value of $\beta \approx 0.865$ (note the similarity to the critical $\beta$ in figure 3a corresponding to lift-off from the $Re=0$ axis) at which the two branches of the neutral curve asymptotically approach each other in the limit of small $\varepsilon$. Curves associated with $\beta$ greater than this critical value are thus observed to turn back at higher $\varepsilon$ (see e.g. the $\beta =0.9$ curves, figure 3c), suggesting that, at finite $Re$, PDI will not exist in the limit $\varepsilon \rightarrow 0$ for all $\beta$. We note another possibility, however, which is that an ‘hourglass’-like pinch-off behaviour occurs, in which the critical curves touch at some finite $\varepsilon$ (appearing here to be around $Sc \approx 10^{4}$), but then separate again for lower $\varepsilon$. The concave-up neutral curves (e.g. $\beta =0.9$) seen here, might then be reflected as concave-down branches at much lower $\varepsilon$ and instabilities for such $\beta$ could then still exist in the $\varepsilon \rightarrow 0$ limit. As $\varepsilon \rightarrow 0$ isolates the instability to an increasingly thin boundary layer at the wall, significantly increased computational power is required to resolve the neutral curves for $Sc > 10^{4}$ and so we do not further consider this behaviour here. It is also worth noting in figure 3(c) that while increasing $Re$ does not significantly influence the horizontal position of the neutral curves along the $W$ axis, the increase in inertial effects does induce a significant downward shift of the curves in $\varepsilon$ (by a factor proportional to $Re$), thus increasing the prevalence of PDI in parameter space.

Figure 4. (a) Dependence of the $\beta =0.87$ neutral curve presented in figure 3(a) on variable $\varepsilon$, with fixed $L=200$. (b) A collapse of the curves in panel (a) based on the inverse Schmidt number $1/Sc = \varepsilon Re$, plotted analogously to figures 2(d) and 3(c) such that small $\varepsilon$ appears at the bottom of the plot. The PCF curves exhibit a perfect collapse and are virtually indistinguishable, whereas the $\varepsilon = 10^{-3}$ and $10^{-4}$ channel curves intersect the $Re=0$ axis in panel (a) and thus diverge to infinite $Sc$ in panel (b). (c) Dependence of the $\beta =0.87$ neutral curve presented in figure 3(a) on variable $L$, with fixed $\varepsilon = 10^{-5}$. (d) Neutral curves in the $\beta$–$W$ plane for fixed $\varepsilon =10^{-5}$ and $Re=1000$.

The streamwise wavenumbers $k$ associated with the FENE-P channel (PPF) neutral curves at $Re=1000$ in figure 3(c) are presented in figure 3(d) as a function of $\varepsilon$, to compare with the scaling of the Oldroyd-B curves in figure 2(c). As for Oldroyd-B, the left-hand branches follow the $1/\sqrt {\varepsilon }$ scaling reported by Beneitez et al. (2023b) for $\varepsilon \lesssim 10^{-5}$, when the neutral curves are roughly independent of $W$ in figure 3(c). In contrast, the right-hand branches follow this $k$ scaling for the entire range of $\varepsilon$ considered here, thus explaining the slight mismatch in the collapse of the right-hand branches of the two geometries at low $\varepsilon$ in figure 3(c) (see e.g. the right-hand branches of the $\beta =0.7$ and 0.8 curves); in this regime, $k$ has become sufficiently small such that the instability is no longer confined to the wall and thus the wall shear $U_{{wall}}'$ is not entirely suitable for scaling the $W$ axis. In figure 3(d), we also note that for $\beta \gtrsim 0.864$, the pinchoff value seen in figure 3(c) at $Sc = O(10^{4})$, the neutral curves turn back to higher $\varepsilon$ before the $1/\sqrt {\varepsilon }$ scaling regime is reached.

In figure 4, we further consider the behaviour of the neutral curves presented in figure 3 by varying parameters that were previously held fixed. While the $Re$–$W$ plane was considered with a fixed $\varepsilon = 10^{-5}$ in figure 3(a), in figure 4(a) we focus on the $\beta =0.87$ neutral curve and demonstrate its dependence on variable $\varepsilon$, where it is found that decreasing $\varepsilon$ pushes the region of PDI to progressively higher $Re$. In figure 4(b), we collapse the curves from figure 4(a) based on the Schmidt number $Sc = 1/(\varepsilon Re)$. A near perfect collapse is observed for the PCF curves, while some channel (PPF) curves diverge to $Sc \rightarrow \infty$ due to their intersection with the $Re=0$ axis in figure 4(a). Additionally, having only considered a fixed polymer extensibility $L=200$ in figure 3, in figures 4(c,d) we consider the influence of variable $L$ on the neutral curves in the $Re$–$W$ and $\beta$–$W$ planes, respectively. In figure 4(c), decreasing $L$ induces a turn-back of the neutral curves at progressively higher values of $Re$, decreasing the extent of PDI in the $Re$–$W$ plane. Figure 4(d), at fixed $Re=1000$, displays the same qualitative features as those reported in the inertialess limit by Beneitez et al. (2023b) (see their figure 2c), again demonstrating that decreasing $L$ decreases the prevalence of PDI in the $\beta$–$W$ plane.

3.3. Growth rates

While we have thus far considered the behaviour of the neutral curves associated with PDI, it is also informative to consider how the growth rate of PDI evolves with $Re$ in regions of instability. Starting with each $\beta$ curve in figure 2(a) (Oldroyd-B), we first fix the value of $W$ at which the neutral curve intersects the $Re=0$ axis, and where the growth rate is thus zero by definition. At this fixed $W$, we then increase $Re$ incrementally and track the growth rate of the most unstable PDI mode, as shown in figure 5(a). In figure 5(b), we repeat this process for the FENE-P curves presented in figure 3(a). As some FENE-P curves do not intersect the $Re=0$ axis, in these cases we fix $W$ to correspond to the minimum $Re$ of the neutral curve, and then increase $Re$ from there. In figures 5(c,d), we again repeat this process for the neutral curves plotted in figures 4(b,c), respectively.

Figure 5. Growth rates $\sigma := kc_{i}$ of the most unstable mode as a function of Reynolds number $Re$ obtained at a fixed $W$ corresponding to either the intersection of the neutral curve with the $Re=0$ axis or the minimum of the neutral curve in $Re$ for cases in which there is no intersection with $Re=0$. Growth rates are shown for neutral curves spanning various $\beta$ in (a) figure 2(a) (Oldroyd-B, fixed $\varepsilon =10^{-5}$) and (b) figure 3(a) (fixed $L=200$, $\varepsilon =10^{-5}$), and neutral curves at fixed $\beta =0.87$ spanning (c) various $\varepsilon$ in figure 4(b) (fixed $L=200$) and (d) various $L$ in figure 4(c) (fixed $\varepsilon =10^{-5}$).

In all cases, for neutral curves that intersect the $Re=0$ axis, the growth rate is observed to grow linearly with $Re$ as one moves away from the neutral curve, emphasizing the intensification of PDI due to the presence of inertia. Notably, the streamwise wavenumber $k$ remains virtually constant during this scaling, until $Re \gtrsim 10^3$ at which point the most unstable $k$ beings to vary and the linear scaling breaks down. For the subset of FENE-P curves that do not intersect the $Re=0$ axis (e.g. $\beta =0.87$, 0.88, 0.90 in figure 4b), it is also intriguing that the growth rates display a dramatic increase with $Re$ to quickly join the linear $Re$ scaling of the curves that do intersect the $Re=0$ axis. We note that the relative vertical translation of the various curves in figure 5 is due to differences in slope between the neutral curves at their intersection with the $Re=0$ axis; at a fixed $W$, these differences in slope will govern the rate at which one moves away from the neutral curve as $Re$ is increased, and hence the observed difference in the growth rate magnitudes.

4. Conclusions

In this study, we have demonstrated that the PDI is active in both plane Couette and channel flows with or without inertial effects, and that the instability intensifies with increasing Reynolds number $Re$. Through exploration of a variety of dimensionless parameters, we have found that PDI is operational across large regions of the parameter space including those relevant to many prior experiments (Choueiri, Lopez & Hof 2018; Qin et al. 2019; Choueiri et al. 2021; Jha & Steinberg 2021). In particular, increasing $Re$ enhances the prevalence of the instability, promoting instability at progressively smaller values of both $W$ and $\varepsilon$ than in the inertialess limit. Our results therefore significantly extend the conclusion of Beneitez et al. (2023b) that PDI could also present a possible transition mechanism to EIT as well as ET in FENE-P fluids.

The eigenfunction for PDI is a wall mode, confined to a boundary layer of thickness $\sqrt {\varepsilon }$. As a result, the neutral curves for plane Couette and channel flow are found to nearly overlap in most regions of parameter space when $W$ is scaled by the wall-shear rate. This collapse breaks down when the streamwise wavenumber $k$ approaches $O(1)$, as the instability is no longer confined to the wall and thus feels a non-monotonic shear profile in the channel's interior, as occurs for large $\varepsilon$ at $Re=0$ (but notably not at higher $Re$, see figure 2c), and small $\varepsilon$ at high $W$ in FENE-P fluids. The finite extensibility of the polymer chains ($L$) is also found to have a significant impact on the prevalence of PDI, as compared with that predicted by the Oldroyd-B model. Considering $L=200$, we found that for sufficiently high $\beta$, the instability is suppressed at $Re=0$ and only appears at progressively larger $Re$ (figure 3a). Similarly, beyond a critical value of $\beta$, the instability may also be suppressed at small values of $\varepsilon$ (figure 3c). Given that PDI emerges as a wall mode, we have also confirmed its presence in cylindrical pipe flow as well as Taylor–Couette flow.

The length scale associated with PDI raises some intriguing questions. As indicated by Beneitez et al. (2023b), PDI can emerge at a length scale roughly of the order of the polymer gyration radius, where the continuum assumptions of the FENE-P model may not hold. It is thus possible that PDI is an unphysical feature of the widely used FENE-P model, which would have significant implications for computational studies of ET, EIT and polymer drag reduction using this model. Until now, such studies may have been unknowingly influenced by PDI, due to the ubiquity of stress diffusion in numerical schemes, either introduced explicitly as a regularization term or arising implicitly through the discretization scheme. Assessing the relevance of PDI to real viscoelastic fluids is now a key challenge to be confronted.