2.1. Governing equations and parameters

We investigate the hydrodynamic stability of incompressible pipe Poiseuille flows based on the Oldroyd-B fluid model (Bird et al. Reference Bird, Curtiss, Armstrong and Hassager1987). The cylindrical coordinate system is used, with $r$, $\theta$ and $z$ denoting the radial, azimuthal and axial directions, respectively. The pipe radius $R$ and the centreline velocity $U_c$ are chosen to be the characteristic length and velocity scales to normalise the system. The non-dimensional perturbation equations are (Wan et al. Reference Wan, Sun and Zhang2021)

(2.1a)\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \quad \partial_t \boldsymbol{u} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{U} + \boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} + \boldsymbol{N}_{\boldsymbol{u}} ={-} \boldsymbol{\nabla} p + \frac{\beta}{Re} \nabla^{2} \boldsymbol{u} + \frac{1-\beta}{ReWi} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{c}, \end{gather}

(2.1b)\begin{gather}{\partial_t} \boldsymbol{c} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{C} - \boldsymbol{c} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{U} - (\boldsymbol{\nabla} \boldsymbol{u})^{T} \boldsymbol{\cdot} \boldsymbol{C} + \boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{c} - \boldsymbol{C} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} - (\boldsymbol{\nabla} \boldsymbol{U})^{T} \boldsymbol{\cdot} \boldsymbol{c} + \boldsymbol{N}_{\boldsymbol{c}} ={-} \frac{\boldsymbol{c}}{Wi}, \end{gather}

where $\boldsymbol {u} = (u_{r}, u_{\theta }, u_{z})^T$ is the perturbation velocity vector, $\boldsymbol {c}=(c_{rr}, c_{r\theta }, c_{rz}, c_{\theta \theta }, c_{\theta z}, c_{zz})^T$ the conformation tensor and $p$ the pressure. $\boldsymbol {U}=(0,0,U_z)$ and $\boldsymbol {C}=(1,0,WiU_z',1,0,1+2Wi^2U_z'^2)$ are the corresponding laminar base states, where $U_z=1-r^2$ and prime $'$ denotes differentiation with respect to $r$. The nonlinear terms are $\boldsymbol {N}_{\boldsymbol {u}} = \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {u}$ and $\boldsymbol {N}_{\boldsymbol {c}} = \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {c} - \boldsymbol {c} \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {u} - (\boldsymbol {\nabla }\boldsymbol {u})^{T} \boldsymbol {\cdot } \boldsymbol {c}$. The controlling parameters include the viscosity ratio $\beta ={\nu _s}/({\nu _s+\nu _p})$ (where $\nu _s$ is solvent viscosity and $\nu _p$ polymer viscosity), the Reynolds number $Re={U_c R}/({\nu _s+\nu _p})$ and the Weissenberg number $Wi={\lambda U_c}/{R}$ ($\lambda$: polymer relaxation time). We define the elasticity number as $E\equiv Wi/Re={\lambda (\nu _s+ \nu _p)}/{R^2}$, characterising the elastic effects of polymers. The no-slip boundary condition $\boldsymbol {u}(r=1)=\boldsymbol {0}$ is applied at the pipe wall. Note that, for the conformation tensor, it is not necessary to specify its boundary condition, because the linear operator in (2.1b) does not contain $r$-derivative terms (like $\partial /\partial _r$) of $\boldsymbol {c}$ (although such terms do exist in the nonlinear term $\boldsymbol {N}_{\boldsymbol {c}}$, which, in a weakly nonlinear framework, is constructed from the linear eigenvectors, as shown below in (A3)). Following previous works (Garg et al. Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018; Chaudhary et al. Reference Chaudhary, Garg, Subramanian and Shankar2021; Zhang Reference Zhang2021), we are interested in the axisymmetric disturbance because only this mode is found to be linearly unstable. Therefore, the symmetric conditions are imposed at the pipe centre, $u_r(r=0)= u'_z(r=0)=0$.

The equation system (2.1) will hereafter be referred to as the original equation system and we use subscript ‘$F$’ to mark it. Introducing $\boldsymbol {\gamma }_{{F}}\!=\!(u_{r{F}},u_{z{F}},p_{{F}},c_{rr{F}},c_{rz{F}},c_{\theta \theta {F}},c_{zz{F}})^T$, the equation system can be recast to a compact form of

(2.2)\begin{equation} (\boldsymbol{M}_{{F}} {\partial_t} -\boldsymbol{L}_{{F}} ) \boldsymbol{\gamma}_{{F}}= \boldsymbol{N}_{{F}}, \end{equation}

where the weight matrix $\boldsymbol {M}_{{F}}$, the linear operator $\boldsymbol {L}_{{F}}$ and the nonlinear operator $\boldsymbol {N}_{{F}}$ can readily be derived from (2.1). The linear mode is governed by the homogeneous system

(2.3)\begin{equation} (\boldsymbol{M}_{{F}} {\partial_t} -\boldsymbol{L}_{{F}} ) \boldsymbol{\gamma}_{{F}}= \boldsymbol{0}, \end{equation}

which admits the normal mode wave-like solution

(2.4)\begin{equation} {\boldsymbol{\gamma}}_{{F}}= \tilde{\boldsymbol{\gamma}}_{{F}}(r)\, {\rm e}^{{\rm i}\alpha (z-c t)} + {\rm c.c.}. \end{equation}

Here, $\tilde {\boldsymbol {\gamma }}_{{F}}$ is the linear eigenfunction, ${\rm i}$ is the imaginary unit, $\alpha$ is the axial wavenumber, $c=c_r + {\rm i} c_i$ is the complex phase speed ($\omega =\alpha c$ is the complex frequency) with $\alpha c_i$ (or $\omega _i$) representing the linear growth rate and ${\rm c.c}.$ represents the complex conjugate of its preceding term. By substituting (2.4) into (2.3), one obtains the following linear eigenvalue problem for the flow:

(2.5)\begin{equation} (c \tilde{\boldsymbol{M}}_{{F}}-\tilde{\boldsymbol{L}}_{{F}}) \tilde{\boldsymbol{\gamma}}_{{F}}=\boldsymbol{0}, \end{equation}

where $\tilde {\boldsymbol {M}}_{{F}}$ and $\tilde {\boldsymbol {L}}_{{F}}$ can be derived from $\boldsymbol {M}_{{F}}$ and $\boldsymbol {L}_{{F}}$, respectively, by replacing $\partial _z$ with ${\rm i} \alpha$.

2.2. Asymptotic analysis

Now, we conduct an asymptotic analysis of the centre-mode instability by assuming $Wi\gg 1$ and $E=O(1)$ ($Re \gg 1$). The polymer concentration is taken to be of an intermediate level, i.e. the viscosity ratio $(\beta,1-\beta )=O(1)$. Solving the linear eigenvalue problem (2.5), we can obtain the neutral curves in the $\alpha$–$Re$ plane. Figure 1 shows an example for $\beta =0.65$, which is adapted from figure 3(c) of Wan et al. (Reference Wan, Sun and Zhang2021). For a given $Wi$, the enclosed region by the solid line represents the unstable zone, and its onset at the lowest $Re$, marked by the red star, represents the linear critical condition. Increase of $Wi$ leads to a higher critical Reynolds number $Re_c$ and a lower critical wavenumber $\alpha _c$, as shown by the red arrow. The control parameters of our interest pertain to these linear critical conditions. From the figure, one can understand that the disturbances to be studied (especially for high $Wi$) are of small $\alpha$, indicating a long-wavelength nature. Moreover, under these linear critical conditions, the elasticity number $E$ is of $O(1)$, as illustrated in figure 11(e) of Wan et al. (Reference Wan, Sun and Zhang2021). In the figure, we also show the neutral curves for fixed $E$ where the linear critical conditions at low $E$, as traced by Chaudhary et al. (Reference Chaudhary, Garg, Subramanian and Shankar2021), are related to the short-wavelength regime and will not be considered here.

Figure 1.

Neutral curves in the $\alpha$–$Re$ plane at viscosity ratio $\beta =0.65$ for viscoelastic pipe flows based on the Oldroyd-B fluid model. The solid curves are computed for fixed $Wi$, with the linear critical conditions ($Re_c$ and $\alpha _c$) marked by the red stars; following the red arrow direction $\alpha _c$ decreases and is smaller than one, implying the long-wavelength nature of the centre mode. The dot-dashed curves are obtained for fixed $E$ (as shown in figure 18 of Chaudhary et al. (Reference Chaudhary, Garg, Subramanian and Shankar2021)). Adapted from figure 3(c) in Wan et al. (Reference Wan, Sun and Zhang2021).

In our asymptotic analysis, we assume $Wi^{-1} \ll \alpha$, which is consistent with the linear critical conditions in figure 1. Thus, a small parameter, $\sigma =(\alpha Wi)^{-1} \ll 1$, can be introduced. The complex phase speed is then expanded as

(2.6)\begin{equation} c=1+\sigma c_1+\cdots, \end{equation}

where $c_1$ is the phase speed correction to be solved for. From the balance of the leading-order terms in the linear system (2.5), we obtain that, in the bulk region where $r=O(1)$, $(\tilde u_{r{F}},\tilde p_{{F}},\tilde c_{rr{F}},\tilde c_{rz{F}},\tilde c_{zz{F}}) \sim (\alpha,1,\sigma ^{-1} ,\alpha ^{-1}\sigma ^{-2}, \alpha ^{-2}\sigma ^{-2})\tilde u_{z{F}}$. For example, from the continuity equation, we know that $\tilde {u}'_{r{F}} \sim {\rm i} \alpha \tilde {u}_{z{F}}$; thus, we have $\tilde {u}_{r{F}} \sim \alpha \tilde {u}_{z{F}}$. The leading-order perturbation field can then be rescaled as

(2.7)\begin{align} (u_{r{F}},u_{z{F}},p_{{F}},c_{rr{F}},c_{rz{F}},c_{zz{F}}) = (\alpha u_{r}, u_{z}, p,\sigma^{{-}1} c_{rr}, \alpha^{{-}1}\sigma^{{-}2} c_{rz}, \alpha^{{-}2}\sigma^{{-}2} c_{zz})+ \cdots. \end{align}

The long-wavelength nature of the instability determines that the radial momentum equation, to the leading order, reduces to $0=-p'_{F}$; so $c_{\theta \theta {F}}$ does not appear in the leading-order equation system. Because $c$ is close to unity, we introduce the relation

(2.8)\begin{equation} {\partial_t}={-}{\partial_z} + \sigma {\partial_\tau}, \end{equation}

where $\tau$ is a time scale related to $c_1$, recalling the expansion of $c$ in (2.6). Noting that in the long-wavelength limit, both $t$, $\tau$ and $z$ are of $O(\alpha ^{-1})$, we introduce

(2.9)\begin{equation} (\bar t, \bar \tau, \bar z)=\alpha (t, \tau, z)=O(1). \end{equation}

With the above assumptions, applying (2.7) in (2.2) and neglecting the $O(Wi^{-1})$ terms, we obtain the following asymptotic equations (referred to as the reduced system):

(2.10a,b)\begin{gather} 0= u_{r}'+ u_{r}/r + \partial_{\bar z} u_{z}, \quad 0={-} p', \end{gather}

(2.10c)\begin{gather}\sigma \partial_{\bar \tau} u_{z}=2r u_{r} + r^2 \partial_{\bar z} u_{z} + \sigma\beta E(u_{z}''+{u_{z}'}/{r}) - \partial_{\bar z} p + (1-\beta)E (c_{rz}'+{c_{rz}}/{r} + \partial_{\bar z}c_{zz}) - N_{u_{z}}, \end{gather}

(2.10d)\begin{gather}\sigma \partial_{\bar \tau} c_{rr} = 2\sigma u_{r}' - 4r\partial_{\bar z}u_{r} + r^2 \partial_{\bar z} c_{rr}-\sigma c_{rr} - N_{c_{rr}}, \end{gather}

(2.10e)\begin{gather}\sigma \partial_{\bar \tau} c_{rz} =2\sigma u_r - 2\sigma r u_{r}' + 8 r^2\partial_{\bar z} u_{r} +\sigma^2 u_{z}' - 2 \sigma r \partial_{\bar z} u_{z} - 2\sigma r c_{rr} + r^2 \partial_{\bar z} c_{rz}-\sigma c_{rz} - N_{c_{rz}}, \end{gather}

(2.10f)\begin{gather}\sigma \partial_{\bar \tau} c_{zz} ={-}16 r u_{r} - 4\sigma r u_{z}' + 16 r^2 \partial_{\bar z} u_{z} - 4 r c_{rz} + r^2 \partial_{\bar z} c_{zz}-\sigma c_{zz} - N_{c_{zz}}, \end{gather}

where the nonlinear terms are

(2.11a)$$\begin{gather} N_{u_{z}} = u_{r} u'_{z} + u_{z} {\partial_{\bar z}} u_{z}, N_{c_{rz}} = u_{r} c'_{rz} + u_{z} {\partial_{\bar z}} c_{rz} - \sigma c_{rr} u'_{z} - c_{rz} {\partial_{\bar z}} u_{z} - c_{rz} u'_{r} - c_{zz} {\partial_{\bar z}} u_{r}, \end{gather}$$

(2.11b)$$\begin{gather}N_{c_{rr}}=u_{r} c'_{rr} + u_{z}{\partial_{\bar z}} c_{rr} \!-\! 2 c_{rr} u'_{r} -2\sigma^{{-}1}c_{rz} {\partial_{\bar z}} u_{r},N_{c_{zz}}\!=\!u_{r} c'_{zz} +u_{z} {\partial_{\bar z}} c_{zz} -\! 2 c_{rz} u'_{z} - 2 c_{zz} {\partial_{\bar z}} u_{z}, \end{gather}$$

with no-slip boundary conditions applied at the pipe wall $r=1$ and symmetric conditions enforced at the pipe axis $r=0$, i.e. $u_r(1)=u_z(1)=0$ and $u_r(0)=u_z'(0)=0$. In (2.10) and (2.11), the $O(\sigma )$ and $O(\sigma ^2)$ terms are retained. If they are neglected in the bulk region, then we would obtain a solution which is singular at both $r=1$ and $r=0$. Thus, a wall layer and a central layer must be taken into account, being similar to the asymptotic structure in Dong & Zhang (Reference Dong and Zhang2022). Such an unstable mode is found to be possible only when $E\sim \sigma ^{-1}$, corresponding to the long-wavelength centre mode for $Re\gg Re_c$ in Dong & Zhang (Reference Dong and Zhang2022). However, when $Re$ is close to $Re_c$, the three layers merge together, and so the $O(\sigma )$ and $O(\sigma ^2)$ terms are kept. It should be noted that, now, we only have $\beta$, $E$ and $\sigma$ as control parameters in (2.10), compared with $\beta, Wi, Re, \alpha$ in the original system (2.5). Introducing $\boldsymbol {\gamma }=(u_{r},u_{z},p,c_{rr},c_{rz},c_{zz})^T$, the compact form of the asymptotic system (2.10) now reads

(2.12)\begin{equation} ( \boldsymbol{M} {\partial_{\bar \tau}} - \boldsymbol{L} ) \boldsymbol{\gamma} =\boldsymbol{N}. \end{equation}

Here, $\boldsymbol {M}$, $\boldsymbol {L}$ and $\boldsymbol {N}$ can be easily deduced by matching (2.12) with (2.10) and thus are not shown here.

2.3. Weakly nonlinear analysis of the asymptotic equation system

Following Wan et al. (Reference Wan, Sun and Zhang2021), we perform a standard multiple-scale expansion of the reduced equation system (2.12). To this end, the following expansions are applied as a series of a small quantity $\delta$:

(2.13a,b)\begin{gather} {\partial_{\bar \tau}} ={\partial_{\bar\tau_{0}}} + \delta {\partial_{\bar\tau_{1}}} + \delta^{2} {\partial_{\bar\tau_{2}}} + O(\delta^{3}), \quad {\partial_{\bar z}} = {\partial_{\bar z_{0}}} + \delta {\partial_{\bar z_{1}}} + O(\delta^{2}), \end{gather}

(2.13c,d)\begin{gather}\boldsymbol{\gamma} = \delta \boldsymbol{\gamma}_{1} + \delta^2 \boldsymbol{\gamma}_{2} + \delta^3 \boldsymbol{\gamma}_{3} + O(\delta^{4}), \quad E=E_c - E_c \delta^2 + O(\delta^4), \end{gather}

where $E_c= Wi/Re_c$ is the linear critical elasticity number. We use the notation $\delta$ to differentiate the present expansion from that in Wan et al. (Reference Wan, Sun and Zhang2021), where the small expansion parameter is $\epsilon$ for the original equation system. The corresponding variables in that paper will be added with a subscript ‘$origi$’ in the following discussion. Their (2.10) is copied below:

(2.14a,b)\begin{gather} {\partial_{t}} ={\partial_{t_{0,{origi}}}} + \epsilon {\partial_{t_{1,{origi}}}} + \epsilon^{2} {\partial_{t_{2,{origi}}}} + O(\epsilon^{3}), \quad {\partial_{z}} = {\partial_{z_{0,{origi}}}} + \epsilon {\partial_{z_{1,{origi}}}} + O(\epsilon^{2}), \end{gather}

(2.14c,d)\begin{gather}\boldsymbol{\gamma} = \epsilon \boldsymbol{\gamma}_{1,{origi}} + \epsilon^2 \boldsymbol{\gamma}_{2,{origi}} + \epsilon^3 \boldsymbol{\gamma}_{3,{origi}} + O(\epsilon^{4}), \quad Re=Re_c + \epsilon^2 + O(\epsilon^4). \end{gather}

The relation between these two expansion methods is briefly explained as follows. From (2.14d), we can get equivalently ${1}/{Re}={1}/{Re_c} - ({1}/{Re_c^2})\epsilon ^{2} + O(\epsilon ^{4})$. Multiplying this expansion with $Wi$ leads to

(2.15)\begin{equation} \frac{Wi}{Re}=\frac{Wi}{Re_c} - \frac{Wi}{Re_c^2}\epsilon^{2} + O(\epsilon^{4}), \quad \text{i.e.} \ E=E_c - E_c\frac{1}{Re_c}\epsilon^{2} + O(\epsilon^{4}). \end{equation}

Comparing the expansion (2.15) with the present expansion (2.13d) shows that these two parameter expansion methods are related by

(2.16)\begin{equation} \delta^2 = Re_c^{{-}1} \epsilon^{2}. \end{equation}

The operators in (2.12), which depend on $\bar \tau$ and $\bar z$, are expanded as $\boldsymbol {L} = \boldsymbol {L}_{0} + \delta \boldsymbol {L}_{1} + \delta ^{2} \boldsymbol {L}_{2} + O(\delta ^{3}), \boldsymbol {N} = \delta ^{2} \boldsymbol {N}_{2} + \delta ^{3} \boldsymbol {N}_{3} + O(\delta ^{4})$. Plugging these expansions along with those in (2.13) into (2.12), and collecting terms of the same order, a series of equations can be obtained:

(2.17a)\begin{gather} ( \boldsymbol{M} {\partial_{\bar \tau_{0}}} - \boldsymbol{L}_{0} ) \boldsymbol{\gamma}_{1} =\boldsymbol{0} \ \text{at order} \ \delta, \end{gather}

(2.17b)\begin{gather}( \boldsymbol{M} {\partial_{\bar \tau_{0}}} - \boldsymbol{L}_{0} ) \boldsymbol{\gamma}_{2} = ( \boldsymbol{L}_{1} - \boldsymbol{M} {\partial_{\bar \tau_{1}}} ) \boldsymbol{\gamma}_{1} + \boldsymbol{N}_{2} \ \text{at order} \ \delta^2, \end{gather}

(2.17c)\begin{gather}( \boldsymbol{M} {\partial_{\bar \tau_{0}}} - \boldsymbol{L}_{0} ) \boldsymbol{\gamma}_{3} = ( \boldsymbol{L}_{2} - \boldsymbol{M} {\partial_{\bar \tau_{2}}} ) \boldsymbol{\gamma}_{1} + ( \boldsymbol{L}_{1} - \boldsymbol{M} {\partial_{\bar \tau_{1}}} ) \boldsymbol{\gamma}_{2} + \boldsymbol{N}_{3} \ \text{at order} \ \delta^3. \end{gather}

Their solutions are assumed to take the wake-like form:

(2.18a)\begin{gather} \boldsymbol{\gamma}_{1} = A_1\tilde{\boldsymbol{\gamma}}_1 E_{{w}} + {\rm c.c.}, \quad \boldsymbol{\gamma}_{2} = |A_1|^2\tilde{\boldsymbol{\gamma}}_{20} + ({ \partial_{\bar{z}_1}}A_1 \tilde{\boldsymbol{\gamma}}_{21}E_{{w}} + {\rm c.c.}) + (A_1^2 \tilde{\boldsymbol{\gamma}}_{22} E_{{w}}^2 + {\rm c.c.}), \end{gather}

(2.18b)\begin{gather}\boldsymbol{\gamma}_{3} = (A_1\tilde{\boldsymbol{\gamma}}_{31,1} +{\partial_{\bar{z}_1 \bar{z}_1}}A_1 \tilde{\boldsymbol{\gamma}}_{31,2} + |A_1|^2 A_1\tilde{\boldsymbol{\gamma}}_{31,3})E_{{w}} + {\rm c.c.} + \cdots, \end{gather}

where $E_{{w}}=\exp ({{\rm i} \bar {z}_{0} - {\rm i} c_1 \bar {\tau }_{0}})$, $\tilde {\boldsymbol {\gamma }}$ with subscripts are the eigenfunctions, and $A_1=A_1(\bar {z}_{1},\bar {\tau }_{1},\bar {\tau }_{2})$ is the complex amplitude of the leading-order wave $\boldsymbol {\gamma }_{1}$. We use $A$ to denote the amplitude of the total disturbance $\boldsymbol {\gamma }$. Then, from the expansion (2.13c), we know that (to the leading order)

(2.19)\begin{equation} A\approx \delta A_1. \end{equation}

Similarly, from (2.14c) one obtains $A\approx \epsilon A_{1,{origi}}$.

Using (2.18) in (2.17) leads to a set of equations (and their complex conjugates which are omitted here) to be solved in the spectral space:

(2.20a)\begin{gather} ( {-{\rm i} c_1} \tilde{\boldsymbol{M}} - \tilde{\boldsymbol{L}}_{0}^{(1)} ) \tilde{\boldsymbol{\gamma}}_1 = {\bf 0}, \end{gather}

(2.20b)\begin{gather}{ \partial_{\bar{z}_1}}A_1 ({-{\rm i} c_1} \tilde{\boldsymbol{M}} - \tilde{\boldsymbol{L}}_{0}^{(1)} ) \tilde{\boldsymbol{\gamma}}_{21}= ( \tilde{\boldsymbol{L}}_{1}^{{\circ}} { \partial_{\bar{z}_1}} - \tilde{\boldsymbol{M}} { \partial_{\bar{\tau}_1}} ) A_1 \tilde{\boldsymbol{\gamma}}_1, \end{gather}

(2.20c)\begin{gather}{}[{(-{\rm i} c_1 + {\rm i} c_1^*)} \tilde{\boldsymbol{M}} - \tilde{\boldsymbol{L}}_{0}^{(0)} ] \tilde{\boldsymbol{\gamma}}_{20} = \tilde{\boldsymbol{N}}_{20}, ({-2{\rm i} c_1} \tilde{\boldsymbol{M}} - \tilde{\boldsymbol{L}}_{0}^{(2)} ) \tilde{\boldsymbol{\gamma}}_{22} = \tilde{\boldsymbol{N}}_{22}, \end{gather}

(2.20d)\begin{gather} \hspace{-5pc}({-{\rm i} c_1} \tilde{\boldsymbol{M}} - \tilde{\boldsymbol{L}}_{0}^{(1)} )(A_1\tilde{\boldsymbol{\gamma}}_{31,1} +{\partial_{\bar{z}_1 \bar{z}_1}}A_1 \tilde{\boldsymbol{\gamma}}_{31,2} + |A_1|^2 A_1\tilde{\boldsymbol{\gamma}}_{31,3}) \nonumber\\ \hspace{5pc} =( \tilde{\boldsymbol{L}}_{2E} - \tilde{\boldsymbol{M}} { \partial_{\bar{\tau}_2}} ) A_1 \tilde{\boldsymbol{\gamma}}_1 + ( \tilde{\boldsymbol{L}}_{1}^{{\circ}} { \partial_{\bar{z}_1}} - \tilde{\boldsymbol{M}} { \partial_{\bar{\tau}_1}} ) { \partial_{\bar{z}_1}}A_1 \tilde{\boldsymbol{\gamma}}_{21} + |A_1|^2 A_1 \tilde{\boldsymbol{N}}_{31}, \end{gather}

where the superscript ${}^*$ denotes the complex conjugate, and the explicit expressions of the various operators are given in Appendix A.

Equation (2.20a) with homogeneous boundary conditions forms an eigenvalue problem. Equation (2.20c) with non-singular linear operators is readily solvable. In order to ensure solutions for (2.20b) and (2.20d), the solvability conditions must be enforced to eliminate the secular terms on their right-hand sides. To do this, the adjoint of the linear problem (2.20a) is introduced as (based on the inner product defined in Wan et al. (Reference Wan, Sun and Zhang2021) via integration by parts Luchini & Bottaro Reference Luchini and Bottaro2014)

(2.21)\begin{equation} ({{\rm i} c_1^*} \tilde{\boldsymbol{M}}^{{{\dagger}}} - \tilde{\boldsymbol{L}}_{0}^{(1){{\dagger}}}) \tilde{\boldsymbol{\gamma}}_1^{{{\dagger}}} = \boldsymbol{0}, \end{equation}

with $\tilde {\boldsymbol {M}}^{{\dagger} }$ and $\tilde {\boldsymbol {L}}_{0}^{(1){{\dagger} }}$ described in Appendix A. Then, the solvability condition applied to (2.20d) leads to

(2.22)\begin{equation} {\partial_{\bar \tau_{2}}}A_1 = a_{1}A_1 + a_{2}{\partial_{\bar{z}_{1} \bar{z}_{1}}}A_1 + a_{3}{|A_1|}^{2}A_1, \end{equation}

where the coefficients are

(2.23a–d)\begin{equation} \left.\begin{array}{c@{}} a_{1} = \dfrac{{\langle \tilde{\boldsymbol{L}}_{2E}\tilde{\boldsymbol{\gamma}}_1, \tilde{\boldsymbol{\gamma}}_{1}^{{{\dagger}}} \rangle}_{s}}{{\langle \tilde{\boldsymbol{M}}\tilde{\boldsymbol{\gamma}}_1, \tilde{\boldsymbol{\gamma}}_{1}^{{{\dagger}}} \rangle}_{s}}, \quad a_{2} = \dfrac{{\langle (\tilde{\boldsymbol{L}}_{1}^{{\circ}} + c_{g}\tilde{\boldsymbol{M}})\tilde{\boldsymbol{\gamma}}_{21}, \tilde{\boldsymbol{\gamma}}_{1}^{{{\dagger}}} \rangle}_{s}}{{\langle \tilde{\boldsymbol{M}}\tilde{\boldsymbol{\gamma}}_1, \tilde{\boldsymbol{\gamma}}_{1}^{{{\dagger}}} \rangle}_{s}},\\ c_{g} ={-}\dfrac{ {\langle \tilde{\boldsymbol{L}}_{1}^{{\circ}} \tilde{\boldsymbol{\gamma}}_{1}, \tilde{\boldsymbol{\gamma}}_{1}^{{{\dagger}}} \rangle}_{s}}{{\langle \tilde{\boldsymbol{M}}\tilde{\boldsymbol{\gamma}}_1, \tilde{\boldsymbol{\gamma}}_{1}^{{{\dagger}}} \rangle}_{s}},\quad a_{3} = \dfrac{{\langle \tilde{\boldsymbol{N}}_{31}, \tilde{\boldsymbol{\gamma}}_{1}^{{{\dagger}}} \rangle}_{s}}{{\langle \tilde{\boldsymbol{M}}\tilde{\boldsymbol{\gamma}}_1, \tilde{\boldsymbol{\gamma}}_{1}^{{{\dagger}}} \rangle}_{s}}. \end{array}\right\} \end{equation}

The operation ${\langle \tilde {\boldsymbol {f}},\tilde {\boldsymbol {g}} \rangle }_{s} = \int _{0}^{1} \tilde {\boldsymbol {f}} \boldsymbol {\cdot } \tilde {\boldsymbol {g}}^* r \, \textrm {d} r$ is an inner product defined in the spectral space. The derived equation (2.22) is the third-order GLE for the asymptotic reduced equations in the large-$Wi$ limit, to be compared with the GLE for the original equations derived in Wan et al. (Reference Wan, Sun and Zhang2021). Of particular interest is the first Landau coefficient $a_3 = a_{3r} + \textrm {i} a_{3i}$ whose real part $a_{3r}$ indicates the primary bifurcation of the laminar flow around its linear critical conditions; its sign being positive (negative) denotes a subcritical (supercritical) bifurcation. To facilitate the comparison and discussion regarding the GLE, we re-write the GLE for the original equations from (2.24) of Wan et al. (Reference Wan, Sun and Zhang2021) (note that $A$ in Wan et al. (Reference Wan, Sun and Zhang2021) is denoted as $A_{1{,origi}}$ here):

(2.24)\begin{align} {\partial_{t_{2,{origi}}}}A_{1,{origi}} = a_{1,{origi}}A_{1,{origi}} + a_{2,{origi}}{\partial_{{z}_{1,{origi}} {z}_{1,{origi}}}}A_{1,{origi}} + a_{3,{origi}}{|A_{1,{origi}}|}^{2}A_{1,{origi}}, \end{align}

where the time scale ${t_{2,{origi}}}$ and space scale ${z}_{1,{origi}}$ have been given in (2.14a,b). This equation can also be expressed in time scale $t$, space scale $z$ with amplitude $A$ by inserting the transformation $\partial _{ t}\approx \epsilon ^2 \partial _{t_{2,{origi}}}$, ${\partial _{zz}} \approx \epsilon ^{2} {\partial _{z_{1,{origi}}z_{1,{origi}}}}$ and that for the disturbance amplitude $A\approx \epsilon A_{1,{origi}}$ into (2.24) as

(2.25)\begin{equation} {\partial_{t}}A = \epsilon^2 a_{1,{origi}}A + a_{2,{origi}}{\partial_{{z}{z}}}A + a_{3,{origi}}{|A|}^{2}A. \end{equation}

An important issue regarding the evaluation of $a_3$ in (2.22) is its uniqueness, i.e. $\tilde {\boldsymbol {\gamma }}_1$ should be normalised to make $a_3$ uniquely determined (Herbert Reference Herbert1980). We follow the normalisation method in Wan et al. (Reference Wan, Sun and Zhang2021), where the linear eigenfunction is normalised so that the square root of the total disturbance energy (kinetic energy plus elastic energy) equals one (see Appendix B). Therefore, the disturbance amplitude $|A|$ (or equivalently $\delta |A_1|$ as in (2.19)) in our expansion has the physical meaning of square root of the total disturbance energy. We hope that this consideration may facilitate comparisons of our results with experiments in the future.

At the end of this section, we would like to discuss the limit of the Oldroyd-B model adopted in this work. The Oldroyd-B model allows for an infinite extension of polymers and does not account for the shear-thinning effects which can be significant at high $Wi$. The more realistic FENE-P model (finitely extensible nonlinear elastic model with Peterlin closure) overcomes these drawbacks. In the FENE-P model, a new parameter $L_{{max}}$ characterising the maximum statistical finite extensibility of polymers is introduced; the model reduces to the Oldroyd-B model when $L_{{max}} \to \infty$. In viscoelastic pipe flows, the centre-mode instability may disappear when $L_{{max}}$ is sufficiently small, as illustrated in table 1 of Zhang (Reference Zhang2021). Without a linear critical condition, the multiple-scale expansion method ceases to work. An alternative is the amplitude expansion method enabling a weakly nonlinear expansion of the disturbance around a weakly damped mode (instead of a neutral mode for the multiple-scale expansion) as in Meulenbroek et al. (Reference Meulenbroek, Storm, Bertola, Wagner, Bonn and van Saarloos2003) (such a scenario is also coined bifurcation from infinity). It would be interesting to see whether scaling laws to be presented exist in the results of the amplitude expansion method.

3.1. Numerical method and code validation

The reduced equations including (2.20a) and its adjoint, (2.20b) and (2.20c) are solved using a spectral collocation method. We avoid placing a grid point at $r=0$ following Mohseni & Colonius (Reference Mohseni and Colonius2000) and construct the differentiation matrices using the even–odd properties of the variables (Trefethen Reference Trefethen2000), i.e. $u_r$ and $c_{rz}$ are odd and $u_z$, $p$, $c_{rr}$ and $c_{zz}$ even for the axisymmetric mode. Validation in solving the original equations (in linear and weakly nonlinear phases) can be found in Wan et al. (Reference Wan, Sun and Zhang2021). Computations of the reduced asymptotic equations have been verified by comparing the results with those from the original equations as supplemented in Appendix C.

3.2. Determining the linear critical conditions

The multiple-scale expansion is commonly performed around linear critical conditions, which guarantees the convergence of the expansion (Fujimura Reference Fujimura1989). Next, we show how to determine these conditions. Figure 2(a) presents the linear growth rate of the most dangerous mode in the linear asymptotic system (2.20a) for a wide range of $\sigma$ at $\beta =0.5$. As we can see, a good agreement of the results from the two equation systems is achieved. Moreover, each curve shows a local peak and the peak corresponds to the linear critical condition when $E=0.16547$ (this critical condition is marked by the red star in panel (b)). Figure 2(b) shows the contours of $c_{1i}$ in the $\sigma$–$E$ plane for $\beta =0.5$, with the black curve being the neutral curve. The red star at the right end of the loop marks the linear critical condition (with $E_c=0.16547$ and $\sigma _c=0.04060$), which agrees well with $Wi/Re_c=0.16543$ and $1/(\alpha _c Wi)=0.04058$ for the original equation system at $\beta =0.5$, $Wi=650$, $Re_c=3929.166947$ and $\alpha _c=0.037909$. Because $E=Wi/Re$ and that in the original equations increasing $Re$ brings out instability, in the asymptotic equations decreasing $E$ renders the flow more unstable from the red star.

Figure 2.

(a) Comparison of $\sigma c_{1i}$ as a function of $\sigma$ at different $E$ and $\beta =0.5$. Results from the original equations are obtained at a large $Wi=650$ and then converted according to $\sigma c_1=\omega /\alpha -1$. (b) Contours of the imaginary part of $c_1$ (i.e. $c_{1i}$) in a $(\sigma,E)$ plane at $\beta =0.5$ with the black line being the neutral curve and the red star marking the linear critical condition. (c) Linear critical conditions $E_c$ and $\sigma _c$ as functions of $\beta$. Note that, in panel (c), $Re_c$ and $\alpha _c$ are first obtained from the original equations at a large $Wi=650$ and then converted to $E_c=Wi/Re_c$ and $\sigma _c=1/(\alpha _c Wi)$ for the comparison with the asymptotic results.

One of the advantages of the reduced system is that the number of governing parameters is reduced. Reflected in determining the critical conditions, we see that, at a given $\beta$, in the original equations the critical condition needs to be determined at each $Wi$, whereas there is only one critical condition in the $\sigma$–$E$ plane in the asymptotic equations. The variations of $E_c$ and $\sigma _c$ with $\beta$ are plotted in figure 2(c), where results obtained from the original equations at $Wi=650$ are converted accordingly and then superposed to illustrate the good agreement. We can also see that with increasing $\beta$, the critical $E_c$ increases whereas the critical $\sigma _c$ decreases.

3.3. Bifurcations and the scaling laws at large $Wi$

3.3.1. Explaining the scaling law of the first Landau coefficient

The GLE (2.22) governs the evolution of disturbance amplitudes around the linear critical conditions. Based on the original equations, Wan et al. (Reference Wan, Sun and Zhang2021) found that both subcritical and supercritical bifurcations exist in axisymmetric viscoelastic pipe flows of Oldroyd-B fluids in a large parameter space, and it is mainly the viscosity ratio $\beta$ (related to polymer concentration) that determines the bifurcation type. They reported a large-$Wi$ scaling law for $a_{3}$ (see their figure 13), which cannot be explained by a simple scaling analysis (as adopted by Chaudhary et al. (Reference Chaudhary, Garg, Subramanian and Shankar2021) to explain some linear scalings in their study). In the following, we use the derived reduced system to illustrate the scaling law.

Figure 3(a) shows the raw data of $a_{3r}$ of the original equations (symbols). After being multiplied by $Wi$ (symbols in panel (b)), $a_{3r}$ either exactly collapse on or gradually approach the $a_{3r}$ of the asymptotic equations (black curve in panel (b)) when $Wi$ increases, for all the $\beta$ investigated. The advantage of the reduced system is more manifest when $\beta$ is larger, which requires an even larger $Wi$ to present the scaling (see the red dot at $\beta =0.9$). In the large-$Wi$ limit, the flow bifurcation type is subcritical at large viscosity ratios $\beta$ (small polymer concentration) and changes to be supercritical when $\beta$ is small (large polymer concentration). The bifurcation boundary is $\beta _{{crit}}\approx 0.785$, which is approximately the same as that in Wan et al. (Reference Wan, Sun and Zhang2021) for large $Wi=650$. The link between $a_3$ obtained from the asymptotic equations and that from the original equations can be more clearly seen in the multiple-scale expansion of these two sets of equations, as follows.

Figure 3.

Scalings for GLE coefficients in the original equations and in the asymptotic equations. (a) Raw data of $a_{3r,{origi}}$; (b) scaled data $a_{3r,{origi}} Wi$ and results of $a_{3r,{asymp}}$ (larger $\beta$ needs higher $Wi$ to accurately present the scaling); (c) scaled data $a_{1r,{origi}}Wi^2$ and results of $a_{1r,{asymp}}E_c$. The imaginary parts of these coefficients also follow these scalings (not shown). In panels (a) and (b) data points above the horizontal dashed line correspond to subcritical bifurcations ($a_{3r}>0$) and those below indicate supercritical bifurcations ($a_{3r}<0$).

Noting that the GLE in (2.22) is in time scale $\bar {\tau }_2$ and space scale $\bar {z}_1$, our first step to build the link is to convert these scales to time scale $t_2$ and space scale $z_1$. To this end, we introduce the following expansions:

(3.1a,b)\begin{gather} {\partial_{\tau}} ={\partial_{\tau_{0}}} + \delta {\partial_{\tau_{1}}} + \delta^{2} {\partial_{\tau_{2}}} + O(\delta^{3}), \quad {\partial_{z}} = {\partial_{z_{0}}} + \delta {\partial_{z_{1}}} + O(\delta^{2}), \end{gather}

(3.1c,d)\begin{gather}{\partial_{\bar t}} ={\partial_{\bar t_{0}}} + \delta {\partial_{\bar t_{1}}} + \delta^{2} {\partial_{\bar t_{2}}} + O(\delta^{3}), \quad {\partial_{t}} ={\partial_{t_{0}}} + \delta {\partial_{t_{1}}} + \delta^{2} {\partial_{t_{2}}} + O(\delta^{3}). \end{gather}

Then, from $\bar t=\alpha t$ in (2.9), we obtain $\alpha \partial _{\bar t}={\partial _t}$. By comparing the terms of $O(\delta ^{2})$ in the expansions of $\partial _{\bar t}$ in (3.1c) and ${\partial _t}$ in (3.1d), we have $\alpha {\partial _{\bar t_{2}}}={\partial _{ t_{2}}}$. From the relations in (2.8) and (2.9), we know that ${\partial _{\bar t}}=-{\partial _{\bar z}} + \sigma {\partial _{\bar \tau }}$; considering the expansions of these partial derivatives in (3.1c) and (2.13a,b), we obtain ${\partial _{\bar t_{2}}}=\sigma {\partial _{\bar \tau _{2}}}$ at $O(\delta ^{2})$. Therefore, $\partial _{\bar \tau _{2}}=\sigma ^{-1}{\partial _{\bar t_{2}}} =\sigma ^{-1}\alpha ^{-1} \partial _{t_2}= (\alpha Wi)\alpha ^{-1} \partial _{t_2} =Wi \partial _{t_2}$. Using this relation in (2.22) results in (note that for the spatial derivative ${\partial _{\bar {z}_{1} \bar {z}_{1}}}=\alpha ^{-2} {\partial _{z_{1} z_{1}}}$ is used)

(3.2)\begin{equation} Wi \partial_{t_2} A_1 = a_{1}A_1 + \alpha^{{-}2} a_{2} {\partial_{z_{1} z_{1}}}A_1 + a_{3}{|A_1|}^{2}A_1. \end{equation}

The second step is to further transform (3.2) into a GLE in time scale $t$ and space scale $z$ with amplitude $A$. Inserting the transformation $\partial _{ t}\approx \delta ^2 \partial _{t_{2}}$, ${\partial _{zz}} \approx \delta ^{2} {\partial _{z_{1}z_{1}}}$ and that for the disturbance amplitude $A\approx \delta A_1$ into (3.2) leads to

(3.3)\begin{equation} \partial_{t} A =Wi^{{-}1}\delta^{2} a_{1}A + Wi^{{-}1}\alpha^{{-}2} a_{2} {\partial_{zz}}A + Wi^{{-}1}a_{3}{|A|}^{2}A. \end{equation}

This equation is to be compared with (2.25). Comparing (3.3) with (2.25) term by term results in (recall the relation in (2.16))

(3.4a)\begin{gather} a_{1,{origi}} =Wi^{{-}1} Re_c^{{-}1} a_{1,{asymp}} = Wi^{{-}2} E_c a_{1,{asymp}}, \end{gather}

(3.4b,c)\begin{gather}a_{2,{origi}} =Wi^{{-}1}\alpha^{{-}2} a_{2,{asymp}}, \quad a_{3,{origi}} = Wi^{{-}1} a_{3,{asymp}}, \end{gather}

where the subscript ‘$asymp$’ is additionally added to denote the coefficients in the asymptotic GLE (2.22), (3.2) and (3.3).

Since $a_{3,{asymp}}$ for a given $\beta$ is fixed at the critical state $(\sigma _c, E_c)$, the relation (3.4$c$) implies $a_{3,{origi}}$ is proportional to $Wi^{-1}$, in agreement with the scaling law observed numerically in Wan et al. (Reference Wan, Sun and Zhang2021) at large $Wi$. Thanks to the asymptotic equation system, we are able to easily identify more scaling results in the weakly nonlinear phase (see (3.4$a,b$)). For example, figure 3(c) shows the scaling law for the coefficient $a_1$ in GLE. We plot $a_{1r,{origi}}Wi^2$ and $a_{1r,{asymp}} E_c$ for comparison. They are found to agree well when $Wi$ is sufficiently large, confirming the relation (3.4$a$) since $Re_c=E_c^{-1}Wi$ by definition, i.e. $a_{1,{origi}}$ is proportional to $Wi^{-2}$.

3.3.2. Scaling law of the equilibrium amplitude $A_e$

According to the third-order GLE (Eq. (2.22) in the present asymptotic analysis and (2.24) in the original equations), finite-amplitude disturbances are required to trigger the subcritical transition, meaning that the flow is linearly stable but could be nonlinearly unstable. In this case, we can use the equilibrium solution of the GLE to quantify the amplitude threshold beyond which the transition occurs. For supercritical bifurcations, infinitesimal disturbances can cause the flow transition and the third-order nonlinearity stabilises the flow, leading to a saturated state whose amplitude could also be characterised by the equilibrium solution of the GLE. Based on the above scaling laws of the coefficients in the GLE, a scaling of the corresponding equilibrium amplitude of disturbance can be derived in the neighbourhood of linear critical conditions either from the original GLE (not shown here) or from the asymptotic GLE. We illustrate the derivation from the asymptotic GLE as follows.

Starting from the asymptotic third-order GLE (3.3), if the diffusion term of the disturbance amplitude is ignored, the equilibrium amplitude $A_e$ (i.e. the modulus of $A$) can be obtained by setting $\partial _t A = 0$ as (note that to calculate $A_e$ only the real parts of the coefficients are used; recall the parameter expansion $E=E_c - E_c \delta ^2 + O(\delta ^4)$ in (2.13d))

(3.5)\begin{align} 0 = \delta^2 a_{1r} |A| + a_{3r}{|A|}^{3} \rightarrow A_e=|A|=\sqrt{-\delta^2 \frac{a_{1r}}{a_{3r}}}=\sqrt{ \frac{a_{1r}}{a_{3r}} \frac{E-E_c}{E_c}}=\sqrt{ \frac{a_{1r}}{a_{3r}} \frac{Re_c-Re}{Re}}. \end{align}

Here, for supercritical bifurcations $Re>Re_c$, $a_{1r}>0$ and $a_{3r}<0$; for subcritical bifurcations $Re< Re_c$, $a_{1r}>0$ and $a_{3r}>0$. Therefore, $A_e$ is always a real number. This corresponds to the finite-amplitude equilibrium solution (near linear criticality) in the axisymmetric viscoelastic pipe flows calculated using GLE. Under the assumptions that $Wi \gg 1$ and $E=O(1)$ (thus $Re \gg 1$) and that the above analysis is restricted to the vicinity of the linear critical conditions ($|Re_c - Re|\ll Re_c$), there is the following scaling law of the equilibrium amplitude $A_e$ in terms of $Wi$:

(3.6)\begin{equation} A_e=\sqrt{ \frac{a_{1r}}{a_{3r}} \frac{Re_c-Re}{Re}} \approx \sqrt{ \frac{a_{1r}}{a_{3r}} \frac{Re_c-Re}{Re_c}} = \sqrt{ \frac{ a_{1r} E_c (Re_c-Re)}{a_{3r} Wi} } \propto Wi^{-{1}/{2}}. \end{equation}

As mentioned in the introduction, in both Newtonian channel and pipe flows, scaling laws of the amplitude threshold have been reported in subcritical regimes and favourable agreements between theoretical predictions (Chapman Reference Chapman2002; Waleffe & Wang Reference Waleffe and Wang2005) and experimental observations (Hof et al. Reference Hof, Juel and Mullin2003; Philip et al. Reference Philip, Svizher and Cohen2007; Lemoult et al. Reference Lemoult, Aider and Wesfreid2012) have been achieved. Scaling laws in viscoelastic flows have also been reported in the literature such as Jovanović & Kumar (Reference Jovanović and Kumar2010) and Morozov & van Saarloos (Reference Morozov and van Saarloos2007), both concerning elastic instabilities at vanishing $Re$. The scaling law of the equilibrium amplitude $A_e \propto Wi^{-1/2}$ derived in this work pertains to the EIT transition (in the parameter range that $Wi \gg 1$, $E=O(1)$). This result can be extended to other flow quantities such as the mean-flow distortion ($\delta ^2|A_1|^2\tilde {\boldsymbol {\gamma }}_{20}\approx |A|^2\tilde {\boldsymbol {\gamma }}_{20}=A_e^2\tilde {\boldsymbol {\gamma }}_{20} \propto Wi^{-1}$ at equilibrium states), which may be of interest to experimentalists who can measure the disturbance amplitude in either the flow field or the conformation tensor field.