Hostname: page-component-6898cfcb55-kwdml Total loading time: 0 Render date: 2026-11-11T09:32:32.450Z Has data issue: false hasContentIssue false
  • English
  • Français

Subcritical and supercritical bifurcations in axisymmetric viscoelastic pipe flows

Published online by Cambridge University Press:  21 October 2021

Dongdong Wan,
Guangrui Sun and
Mengqi ZhangOpen the ORCID record for Mengqi Zhang[Opens in a new window]
Dongdong Wan
Affiliation:
Department of Mechanical Engineering, National University of Singapore, 9 Engineering Drive 1, 117575 Singapore
Guangrui Sun
Affiliation:
Department of Mechanical Engineering, National University of Singapore, 9 Engineering Drive 1, 117575 Singapore
Mengqi Zhang*
Affiliation:
Department of Mechanical Engineering, National University of Singapore, 9 Engineering Drive 1, 117575 Singapore
*
†Email address for correspondence: mpezmq@nus.edu.sg
Rights & Permissions [Opens in a new window]

Abstract

Axisymmetric viscoelastic pipe flow of Oldroyd-B fluids has been recently found to be linearly unstable by Garg et al. (Phys. Rev. Lett., vol. 121, 2018, 024502). From a nonlinear point of view, this means that the flow can transition to turbulence supercritically, in contrast to the subcritical Newtonian pipe flows. Experimental evidence of subcritical and supercritical bifurcations of viscoelastic pipe flows have been reported, but these nonlinear phenomena have not been examined theoretically. In this work, we study the weakly nonlinear stability of this flow by performing a multiple-scale expansion of the disturbance around linear critical conditions. The perturbed parameter is the Reynolds number with the others being unperturbed. A third-order Ginzburg–Landau equation is derived with its coefficient indicating the bifurcation type of the flow. After exploring a large parameter space, we found that polymer concentration plays an important role: at high polymer concentrations (or small solvent-to-solution viscosity ratio $\beta \lessapprox 0.785$), the nonlinearity stabilizes the flow, indicating that the flow will bifurcate supercritically, while at low polymer concentrations ($\beta \gtrapprox 0.785$), the flow bifurcation is subcritical. The results agree qualitatively with experimental observations where critical $\beta \approx 0.855$. The pipe flow of upper convected Maxwell fluids can be linearly unstable and its bifurcation type is also supercritical. At a fixed value of $\beta$, the Landau coefficient scales with the inverse of the Weissenberg number ($Wi$) when $Wi$ is sufficiently large. The present analysis provides a theoretical understanding of the recent studies on the supercritical and subcritical routes to the elasto-inertial turbulence in viscoelastic pipe flows.

JFM classification

Non-Newtonian Flows: ViscoelasticityNonlinear Dynamical Systems: Bifurcation

Information

Type
JFM Papers
Information
Journal of Fluid Mechanics , Volume 929 , 25 December 2021, A16
Check for updates
Copyright
© The Author(s), 2021. Published by Cambridge University Press

1. Introduction

Upon a minute addition of polymers, turbulent drag reduction can occur in viscoelastic polymeric flows. This feature has attracted much attention from researchers for a long time. As one of the challenging problems, transition to turbulence in viscoelastic pipe flow has not been fully understood so far. Newtonian pipe flows are linearly stable at very high Reynolds numbers (Davey & Drazin Reference Davey and Drazin1969; Meseguer & Trefethen Reference Meseguer and Trefethen2003) (in fact, the current consensus holds that this type of flow is unconditionally linearly stable although there is no rigorous mathematical proof of this); thus, at Reynolds numbers explored so far, the laminar Newtonian pipe flow can only transition to turbulence via a subcritical route (meaning that nonlinearity destabilizes the flow because of the lack of a linear instability mechanism). It has been recently found that the balance of inertial, viscous and elastic effects can lead to a linear instability of the centre mode in viscoelastic pipe flows (of Oldroyd-B fluids), providing a possible supercritical transition route to turbulence in this flow (Garg et al. Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018). Much attention has been diverted to this novel instability; nevertheless, the subsequent nonlinear development of the centre-mode instability has not been well explored so far and the bifurcation type of this flow has not been investigated extensively in a large parameter space, both of which will be the topics of the present study. In the following subsections, we will first review the relevant works in the literature on viscoelastic shear flows and then discuss the position of the current work.

1.1. Linear stability/instability in viscoelastic shear flows

In the linear stability analysis of viscoelastic parallel shear flows, the plane Couette flow of upper convected Maxwell fluids (UCM, which is a limiting case of the Oldroyd-B model without solvent viscosity) was first studied, pioneered by the milestone work of Gorodtsov & Leonov (Reference Gorodtsov and Leonov1967). Many subsequent works confirmed the linear stability of this flow in the whole parameter space of Reynolds number $Re$ and Weissenberg number $Wi$ (Lee & Finlayson Reference Lee and Finlayson1986; Renardy Reference Renardy1992; Chokshi & Kumaran Reference Chokshi and Kumaran2009; Chaudhary et al. Reference Chaudhary, Garg, Shankar and Subramanian2019). In contrast, the plane Poiseuille flow of UCM fluids was found to present a more complex linear stability/instability diagram. The first investigation into the linear stability of this flow was carried out by Porteous & Denn (Reference Porteous and Denn1972). Their analysis at high $Re>1400$ showed that polymer elasticity (characterized by the elasticity number $E=Wi/Re$) destabilizes the flow, as the critical Reynolds number $Re_c$ decreases monotonically with increasing $E$ up to 0.0025. This effect is more profound for sufficiently high $E$ at which they found two new unstable modes, qualitatively different from the elastically modified Tollmien–Schlichting (TS) mode. After a more extensive exploration up to $E=0.0030$, the minimum $Re_c$ was found to exist at about $E=0.0025$ (Sureshkumar & Beris Reference Sureshkumar and Beris1995b). At low $Re$, however, the plane Poiseuille flow of highly elastic UCM fluids was found to be linearly stable (Ho & Denn Reference Ho and Denn1977; Lee & Finlayson Reference Lee and Finlayson1986), implying the necessity of a nonlinear instability to account for the ‘melt fracture’ phenomenon observed in a polymer extrusion process. The importance of such nonlinearity was indeed confirmed in a weakly nonlinear stability analysis by Meulenbroek et al. (Reference Meulenbroek, Storm, Morozov and van Saarloos2004), to be discussed below. In terms of the symmetry properties of eigenfunctions, all these studies found merely antisymmetric unstable modes, while a more recent comprehensive linear stability analysis performed by Chaudhary et al. (Reference Chaudhary, Garg, Shankar and Subramanian2019) revealed up to seven unstable modes with some of them being symmetric after an extensive search in the parameter space of the plane Poiseuille flow of UCM fluids. They interpreted these unstable modes as part of an infinite hierarchy of elasto-inertial instabilities resulting from the competing effects between polymer elasticity and fluid inertia.

The Oldroyd-B model, compared with the UCM model, more realistically models the effect of solvent viscosity, quantified by the ratio $\beta$ of solvent viscosity to total viscosity. By definition, the Oldroyd-B model reduces to the UCM model at $\beta =0$ and Newtonian fluid at $\beta =1$. In the same study mentioned above, Sureshkumar & Beris (Reference Sureshkumar and Beris1995b) also investigated the linear stability of plane Poiseuille flow of Oldroyd-B fluids. They reported a pronounced stabilizing effect of non-zero solvent viscosity (i.e.  $\beta >0$) at $E=0.001$, with the $Re_c$ increasing monotonically from the UCM limit to the Newtonian limit. At a fixed viscosity ratio ($\beta =0.5$ in their study), their $Re_{c}$-$E$ plot also presents a minimum similar to what they found for UCM fluids. The mechanism underlying such a minimum was revealed years later and it involves two competing contributions to the transport of perturbation vorticity at highly elastic regimes: one from the shear stress, stabilizing the flow, and the other from the normal stress, destabilizing the flow (Sadanandan & Sureshkumar Reference Sadanandan and Sureshkumar2002). Zhang et al. (Reference Zhang, Lashgari, Zaki and Brandt2013) proposed a dimensionless parameter characterizing the non-monotonic effect of elasticity on the linear stability, that is, the ratio of the polymer relaxation time to the characteristic instability time scale, applicable to both short-time and long-time horizons. In addition to these wall-mode instabilities reviewed above, a new centre-mode instability in plane Poiseuille flow of Oldroyd-B fluids has been reported recently by Khalid et al. (Reference Khalid, Chaudhary, Garg, Shankar and Subramanian2021). This novel mode only exists for dilute polymer concentrations (with the viscosity ratio $\beta >0.5$) and highly elastic fluids (with the elasticity number $E>0.01$). Moreover, the $Re_c$ they calculated is in a qualitative agreement with that estimated in the experiment conducted by Srinivas & Kumaran (Reference Srinivas and Kumaran2017).

Unlike the linear stability analyses of viscoelastic plane Poiseuille flows, very few studies have focused on viscoelastic pipe flows. Hansen (Reference Hansen1973) addressed the linear stability problem of the pipe flow with simplified UCM and Oldroyd-B fluid models. The author found that a polymer additive has a stabilizing effect on disturbances when the polymer relaxation time is small, while a destabilizing effect appears if this relaxation time exceeds a critical value. However, this analysis was restricted to high-phase-velocity axisymmetric disturbances and the conclusions were built upon an oversimplification of UCM and Oldroyd-B models by ignoring the convected derivatives in the polymer constitutive equations. Although no linear instability was explicitly found in this study, the author held that the ‘destabilizing effect’ identified might support the viewpoint that early turbulence observed in the viscoelastic pipe flow may possibly result from hydrodynamic instability (Forame, Hansen & Little Reference Forame, Hansen and Little1972). Such a speculation has proved to be reasonable in the work by Garg et al. (Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018) based on the complete UCM and Oldroyd-B models, as summarised below.

Using both the spectral collocation method and shooting method to solve the linear eigenvalue problem of the viscoelastic pipe flow, Garg et al. (Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018) identified an unstable mode travelling at a phase speed close to the maximum base flow, in contrast to the wall mode. The authors interpreted this unstable centre mode as an instability in velocity field reinforced by the polymeric force localized near the pipe centreline. The Oldroyd-B model with $\beta \to 0$ reaches the UCM limit, at which they found that the linear instability is absent, different from the plane Poiseuille flow of UCM fluids where linear instability exists (Porteous & Denn Reference Porteous and Denn1972; Sureshkumar & Beris Reference Sureshkumar and Beris1995b; Chaudhary et al. Reference Chaudhary, Garg, Shankar and Subramanian2019). The authors concluded that the novel linear instability is a subtle balance among fluid inertia, viscosity and polymer elasticity in viscoelastic pipe flows. Recently, the same research group extended their work to the exploration of instability in a larger parameter space and compared their predictions with direct numerical simulation (DNS) results and experimental measurements (Chaudhary et al. Reference Chaudhary, Garg, Subramanian and Shankar2021). Some of their new observations include that the eigenfunction of unstable modes is not localized near the pipe centreline as $\beta$ approaches one with $E$ of the order of 0.1; and that the laminar flow remains stable at quite low $Re$ regardless of the values of $E$ and $\beta$. However, one of their comparisons showed that there is one order of magnitude difference between the elasticity number $E$ at which they predicted instability and the $E$ measured in experiments near the onset of turbulence (Samanta et al. Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013). They attributed this discrepancy to the characterization procedure in measuring the polymer relaxation time in experiments. Nevertheless, the comparison of the predicted $Re_c$ with the transitional Reynolds number $Re_t$ observed in an experiment carried out by Chandra, Shankar & Das (Reference Chandra, Shankar and Das2018) showed general consistency.

Upon considering the shear-thinning effect by using the FENE-P model (finitely extensible nonlinear elastic model with Peterlin closure) in their scaling analysis, Chaudhary et al. (Reference Chaudhary, Garg, Subramanian and Shankar2021) were able to derive $Re_c \sim (E(1-\beta ))^{-5/8}$, which is consistent with the scaling $Re_t \sim (E(1-\beta ))^{-1/2}$ observed in the experiment. However, in the FENE-P model if the polymer maximum extensibility is small enough (when this parameter approaches infinity, the FENE-P model reduces to the Oldroyd-B model), the linear instability ceases to exist (Zhang Reference Zhang2021). The existence of such linear instability in viscoelastic pipe flow indicates the possibility of a supercritical bifurcation route. The bifurcation type can be more systematically studied in a weakly nonlinear analysis framework (as will be reviewed next) to further understand the nonlinear development.

1.2. Weakly nonlinear stability/instability in Newtonian and viscoelastic flows

When the amplitude of the disturbance increases to a certain degree, the linear framework becomes inapplicable and the effect of nonlinearity starts to manifest itself. In principle, it is difficult to analyse the nonlinearity, especially when it is strong, but if the weakly nonlinear phase is of interest, the flow can be studied in an analytical manner by applying a multiple-scale expansion method.

The weakly nonlinear stability theory has initially been developed in the context of Newtonian plane shear flows (Landau Reference Landau1944; Stuart Reference Stuart1960; Reynolds & Potter Reference Reynolds and Potter1967; Herbert Reference Herbert1983; Fujimura Reference Fujimura1989), reaching the conclusion that the transition in plane Poiseuille flow is subcritical, i.e. the nonlinearity will destabilize the flow just beyond the linear criticality. The weakly nonlinear stability theory was conventionally applied around the linear criticality in order to guarantee the convergence of the expansion scheme. Thus, historically, some controversies were caused when this theory was applied to the Newtonian pipe flow, which is linearly stable at all $Re$ investigated (Davey & Drazin Reference Davey and Drazin1969; Meseguer & Trefethen Reference Meseguer and Trefethen2003). Finite-amplitude equilibrium solutions were reported by Davey & Nguyen (Reference Davey and Nguyen1971) who adopted an equilibrium amplitude method (based on the false problem method developed by Reynolds & Potter Reference Reynolds and Potter1967) to analyse the axisymmetric Newtonian pipe Poiseuille flow. On the contrary, Itoh (Reference Itoh1977) found no such equilibrium solutions using Stuart's method. The issue was then partially solved by Davey (Reference Davey1978) who pointed out the major weakness of both expansion methods for the problems without neutral curves. Later, Patera & Orszag (Reference Patera and Orszag1981) conducted DNS and confirmed that there is no finite-amplitude equilibrium in the axisymmetric Newtonian pipe flow. From this perspective, the newly found linear instability in Garg et al. (Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018) can facilitate the application of the expansion method in viscoelastic pipe flows in general, because linear critical conditions now exist in this flow and can guarantee the convergence of the expansion scheme.

In addition to the Newtonian fluids, weakly nonlinear stability analyses have also proven to be useful in revealing the bifurcation for more complex fluids. For example, Bouteraa et al. (Reference Bouteraa, Nouar, Plaut, Métivier and Kalck2015) and Bouteraa & Nouar (Reference Bouteraa and Nouar2015) adopted an amplitude expansion method to investigate the bifurcation nature of Rayleigh–Bénard convention in shear-thinning fluids between two horizontal plates. In the viscoelastic flows, the corresponding weakly nonlinear stability analysis has been mainly performed by Morozov & van Saarloos and their co-workers. Their efforts aimed at understanding the melt fracture instability occurring in polymer extrusion through a die. In this case, the flow is dominated by polymer elasticity and also linearly stable. In the zero-$Re$ limit, Meulenbroek et al. (Reference Meulenbroek, Storm, Bertola, Wagner, Bonn and van Saarloos2003) studied the weakly nonlinear stability of pipe Poiseuille flow of UCM fluids and reported a subcritical instability for $Wi>5$. This subcritical mechanism has been experimentally demonstrated by Bonn et al. (Reference Bonn, Ingremeau, Amarouchene and Kellay2011) who observed large velocity fluctuations in this flow, supporting the existence of elastic turbulence at very low Reynolds numbers. An extension to plane Poiseuille flow showed that the subcritical instability exists for $Wi>4$ (Meulenbroek et al. Reference Meulenbroek, Storm, Morozov and van Saarloos2004). Both these subcritical instabilities are believed to be intrinsic routes to the melt fracture phenomenon (see Bertola et al. (Reference Bertola, Meulenbroek, Wagner, Storm, Morozov, van Saarloos and Bonn2003) for the corresponding experimental evidences). More recently, Morozov & van Saarloos (Reference Morozov and van Saarloos2019) further extended their analysis in the zero-$Re$ limit to the plane Poiseuille flow of Oldroyd-B fluids and again the instability was found to be subcritical. In addition to these weakly nonlinear analyses of Poiseuille flows, subcriticality was also found to exist in plane Couette flow of UCM fluids at very small Reynolds numbers (Morozov & van Saarloos Reference Morozov and van Saarloos2005). A more comprehensive discussion on these subcritical mechanisms in viscoelastic parallel shear flows can be found in their introductory essay (Morozov & van Saarloos Reference Morozov and van Saarloos2007). Overall, the works summarised above focused on the inertialess limit, relevant to the flow problems of their interest. The present work, instead, aims at performing weakly nonlinear stability analyses of the viscoelastic pipe flows in the elasto-inertial regime, where drag reduction happens.

1.3. Recent discussions on the bifurcation type in viscoelastic flows

Recently, in the research community of viscoelastic flows, elasto-inertial turbulence (EIT) has received much attention owing to its unique flow features relevant to the MDR (maximum drag reduction, an asymptotic statistical state first discussed by Virk Reference Virk1975), and has been investigated experimentally and numerically by many researchers (Dubief, Terrapon & Soria Reference Dubief, Terrapon and Soria2013; Samanta et al. Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013; Sid, Terrapon & Dubief Reference Sid, Terrapon and Dubief2018; Lopez, Choueiri & Hof Reference Lopez, Choueiri and Hof2019; Shekar et al. Reference Shekar, McMullen, Wang, McKeon and Graham2019,  Reference Shekar, McMullen, McKeon and Graham2020; Page, Dubief & Kerswell Reference Page, Dubief and Kerswell2020; Choueiri et al. Reference Choueiri, Lopez, Varshney, Sankar and Hof2021; Shekar et al. Reference Shekar, McMullen, McKeon and Graham2021), among many others. In the DNS study by Page et al. (Reference Page, Dubief and Kerswell2020) the exact coherent structures were calculated in subcritical two-dimensional viscoelastic channel flows of FENE-P fluids and the authors demonstrated that the flow possesses a subcritical transition mechanism in terms of both Reynolds number and Weissenberg number. As for the bifurcation type of pipe flows, since the transition to turbulence in Newtonian pipe flow is known to be subcritical, a naive speculation would be to argue that the viscoelastic pipe flow will also transition to turbulence only via subcritical routes. However, the finding of the centre-mode instability enriches the picture (Garg et al. Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018; Chaudhary et al. Reference Chaudhary, Garg, Subramanian and Shankar2021), implying and confirming the possibility of a supercritical transition route as envisioned by Graham (Reference Graham2014).

Samanta et al. (Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013) first studied experimentally the EIT phenomenon in a viscoelastic pipe flow at a high polymer concentration of 500 p.p.m. polyacrylamide solution, and found that the flow became turbulent at the same $Re$ (${\approx }800$) regardless of the perturbation level. Even though this non-hysteresis behaviour can be considered to support a supercritical transition (see also the discussion in Garg et al. Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018), the authors noted that a subcritical transition mechanism could not be ruled out, because it is difficult to reduce the disturbance in experiments to a quite low level for the possible subcriticality to manifest itself in this flow (indicating that the viscoelastic pipe flows may be very sensitive to disturbance). When the concentration is lower, Samanta et al. (Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013) found a clear hysteresis loop for 100 p.p.m. solutions and the flow is Newtonian-like, signifying a subcritical transition. The occurrence of such hysteresis seems to be dependent on the polymer concentration (related to viscosity ratio $\beta$), which implies that there may exist a boundary between supercritical and subcritical transitions according to this parameter. In a recent experiment of pipe flows of 600 p.p.m. polymer solutions (Choueiri et al. Reference Choueiri, Lopez, Varshney, Sankar and Hof2021), the measured pressure fluctuation amplitude $p_m$ grows continuously with increasing $Re$ near the instability onset, following a scaling of $p_m \propto \sqrt {Re-Re_c}$. Even though this continuous change is also suggestive of a supercritical transition, the authors warned the reader that other scaling relations may exist, considering experimental uncertainties. On the other hand, they revealed a disordered chevron pattern flow at $Re\sim 5$ while the lowest $Re_c$ predicted by the linear stability theory is of $O(80)$. This observation, along with the weakly chaotic fluid motion in the chevron pattern, is believed to be consistent with a subcritical scenario. However, due to experimental uncertainties, it is difficult to confirm the existence of this subcriticality. More coordinated comparisons between experiments and theoretical analyses need to be conducted in order to dispel the doubts on parameter choices and experimental uncertainties.

The current situation calls for a systematic investigation of the flow bifurcation in viscoelastic pipe flows from the perspective of governing equations. The possible existence of both subcriticality and supercriticality in a large parameter space in the reviewed experiments and numerical simulations implies a bifurcation boundary in viscoelastic pipe flows. This view seems to be supported also by the experimental observations in Chandra, Shankar & Das (Reference Chandra, Shankar and Das2020) where they suggested a possible crossover of the transition type from subcritical at low polymer concentrations to supercritical at high polymer concentrations in microtubes. We are motivated by these recent studies to conduct weakly nonlinear analyses of viscoelastic pipe flows to distinguish the two types of transition and identify the boundary between them.

1.4. The position of the current work

As reviewed above, there is currently no theoretical work studying the bifurcation mechanism near linear critical conditions in viscoelastic pipe flows; both supercritical and subcritical transitions have been observed in experiments where disturbances are likely of finite amplitude. It is difficult to distinguish between the genuine supercritical bifurcation and the subcritical bifurcation that is very sensitive to system-level disturbances (this kind of subcritical bifurcation may appear to be supercritical bifurcation in experiments as the laboratory background noise makes it difficult to differentiate the two; see Appendix A for an illustration). Such difficulty necessitates a systematic theoretical investigation. In the current work, we will perform a weakly nonlinear analysis of axisymmetric viscoelastic pipe flows (of Oldroyd-B fluids) based on multiple-scale expansion around the linear critical point (the damped mode was investigated in Meulenbroek et al. Reference Meulenbroek, Storm, Bertola, Wagner, Bonn and van Saarloos2003). The significance of the current work lies in (1) supplementing the works by Morozov & van Saarloos on the weakly nonlinear stability analysis of viscoelastic plane shear flows; (2) extending the works by Garg et al. (Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018) and Chaudhary et al. (Reference Chaudhary, Garg, Subramanian and Shankar2021) to the study of the weakly nonlinear phase of the linearly unstable centre mode; and (3) providing a theoretical investigation of the bifurcation types in viscoelastic pipe flows in a large parameter space to understand and reconcile the experimental observations and numerical results in the studies reviewed above.

The rest of the paper is organized as follows. Section 2 introduces the governing equation in two different formulations, the non-dimensional control parameters, the multiple-scale expansion method in the weakly nonlinear framework and the resulting Ginzburg–Landau equation. In § 3 the numerical method used to evaluate the Landau coefficients is briefly presented. We show the results in § 4, including a validation step, neutral curves in linear stability analysis, effects of nonlinearity in weakly nonlinear stability analysis, the bifurcation type of the flow and a scaling law of the Landau coefficient. We conclude the paper in § 5 with some discussions on the results. In the five appendices, we provide more information on the bifurcation types, the linear/nonlinear operators in the weakly nonlinear stability theory, validation by a DNS method and more results on the UCM pipe flows.

2. Problem formulation

2.1. Governing equations and parameters

We consider viscoelastic fluids in a circular pipe with density $\rho ^{*}$, dynamic viscosity of the solvent $\mu _{s}^{*}$ and the additional dynamic viscosity $\mu _{p}^{*}$ due to polymers, where ${}^*$ is used to indicate dimensional quantities. The governing equations are the incompressible Navier–Stokes equation and the constitutive equation modelling the polymer dynamics. The characteristic scales used to non-dimensionalise the governing equations include: the pipe radius $R^{*}$ as the length scale, the centreline velocity of the laminar flow $U_{c}^{*}$ as the velocity scale, $\rho ^* U_{c}^{*2}$ as the reference pressure and $\mu _{p}^{*} U_{c}^{*}/R^{*}$ for the polymeric shear stress. Then using a hat to denote a non-dimensionalised variable, we can write the non-dimensional continuity and Navier–Stokes equations as

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

where $\hat {\boldsymbol {u}} = (\hat {u}_{r}, \hat {u}_{\theta }, \hat {u}_{z})$ is the velocity vector with the subscripts $r$, $\theta$ and $z$ representing the radial, azimuthal and axial directions, respectively; $\hat {p}$ is the pressure; $\hat {\boldsymbol {\tau }}_{p}$ is the polymeric stress tensor, to be discussed below. The Reynolds number is defined as $Re = \rho ^{*} U_{c}^{*} R^{*} / \mu ^{*}$, and the viscosity ratio as $\beta = \mu _{s}^{*} / \mu ^{*}$, with $\mu ^{*} = \mu _{s}^{*} + \mu _{p}^{*}$ being the total dynamic viscosity.

To model the polymeric stress $\hat {\boldsymbol {\tau }}_{p}$, the Oldroyd-B model (Bird et al. Reference Bird, Curtiss, Armstrong and Hassager1987) is adopted in the present study. In this model, polymer chains are treated as non-interacting dumbbells with two beads connected with Hookean springs, described by an end-to-end vector $\boldsymbol {q}^{*}$. Coarse-graining modelling introduces the conformation tensor $\boldsymbol {c}^{*}=\langle \boldsymbol {q}^{*}\boldsymbol {q}^{*}\rangle$ to characterize the configuration of the polymers, where $\langle \cdot \rangle$ denotes ensemble average. The polymeric stress can then be expressed as $\boldsymbol {\tau }_{p}^{*} = ( \boldsymbol {c}^{*} - \boldsymbol {c}^{*eq} )\mu _{p}^{*} H^{*}/(\lambda ^{*} k_{B}^{*} T^{*})$, where $\lambda ^{*}$ is the polymer relaxation time, $k_{B}^{*}$ the Boltzmann constant, $T^{*}$ the absolute temperature and $H^{*}$ the spring constant. Normalizing $\boldsymbol {c}^{*}$ by $k_{B}^{*}T^{*}/H^{*}$ results in the non-dimensional expression $\hat {\boldsymbol {\tau }}_{p} = (\hat {\boldsymbol {c}}-\boldsymbol {I})/{Wi}$, where $\boldsymbol {I}$ denotes the identity matrix, corresponding to the equilibrium state $\boldsymbol {c}^{*eq}$, and the Weissenberg number is defined as $Wi = \lambda ^{*} U_{c}^{*} /R^{*}$, quantifying the polymer relaxation time to the flow turn-over time. Consequently, the evolution equation of the conformation tensor $\hat {\boldsymbol {c}}$ which has six components $(\hat {c}_{rr}, \hat {c}_{r\theta }, \hat {c}_{rz}, \hat {c}_{\theta \theta }, \hat {c}_{\theta z}, \hat {c}_{zz})$ in non-dimensional form reads as

(2.2)\begin{equation} \frac{\partial \hat{\boldsymbol{c}}}{\partial t} + \hat{\boldsymbol{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} \hat{\boldsymbol{c}} - \hat{\boldsymbol{c}} \boldsymbol{\cdot} (\boldsymbol{\nabla} \hat{\boldsymbol{u}}) - (\boldsymbol{\nabla} \hat{\boldsymbol{u}})^{\textrm{T}} \boldsymbol{\cdot} \hat{\boldsymbol{c}} ={-} \hat{\boldsymbol{\tau}}_{p}. \end{equation}

Despite the simplicity of the Oldroyd-B model, it has been reported to be able to reliably reproduce the purely elastic instabilities observed in experiments of viscometric flows (Shaqfeh Reference Shaqfeh1996) and discover the centre-mode linear instability in the stability analysis of viscoelastic pipe flows (Garg et al. Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018). This model has been used in many previous works on the stability analyses of viscoelastic flows (Sureshkumar & Beris Reference Sureshkumar and Beris1995b; Morozov & van Saarloos Reference Morozov and van Saarloos2007; Zhang et al. Reference Zhang, Lashgari, Zaki and Brandt2013; Morozov & van Saarloos Reference Morozov and van Saarloos2019; Khalid et al. Reference Khalid, Chaudhary, Garg, Shankar and Subramanian2021).

The governing equations (2.1a,b) and (2.2) admit a steady solution and it will hereafter be referred to as the laminar base flow. To perform the stability analysis, the state variables are decomposed as $\hat {\boldsymbol {u}} = \boldsymbol {U} + \boldsymbol {u}$, $\hat {p} = P + p$ and $\hat {\boldsymbol {c}} = \boldsymbol {C} + \boldsymbol {c}$, where the uppercase variables correspond to profiles of the laminar base flow and those lowercase variables without hat symbols denote the perturbations (note that the Reynolds decomposition of $\hat {\boldsymbol {c}}$ is formally presented here, but a geometric decomposition will be introduced later for $\hat {\boldsymbol {c}}$ in a different formulation of the same problem). By substituting this decomposition into (2.1a,b) and (2.2) and then subtracting equations for the laminar base flow, we can obtain the following nonlinear evolution equations for perturbations:

(2.3a)\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \quad \frac{\partial \boldsymbol{u}}{\partial t} + \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.3b)\begin{gather}\frac{\partial \boldsymbol{c}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{C} - \boldsymbol{c} \boldsymbol{\cdot} (\boldsymbol{\nabla} \boldsymbol{U}) - (\boldsymbol{\nabla} \boldsymbol{u})^{\textrm{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})^{\textrm{T}} \boldsymbol{\cdot} \boldsymbol{c} \!+\! \boldsymbol{N}_{\boldsymbol{c}} ={-} \frac{\boldsymbol{c}}{Wi}. \end{gather}

Here the nonlinear terms are $\boldsymbol {N}_{\boldsymbol {u}} = \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {u}$, $\boldsymbol {N}_{\boldsymbol {c}} = \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla } \boldsymbol {c} - \boldsymbol {c} \boldsymbol {\cdot } (\boldsymbol {\nabla } \boldsymbol {u}) - (\boldsymbol {\nabla }\boldsymbol {u})^{\textrm {T}} \boldsymbol {\cdot } \boldsymbol {c}$, and the laminar base flow profiles are

(2.4a)\begin{gather} U_r=U_{\theta}=0, \quad U_z = 1-r^2, \end{gather}
(2.4b)\begin{gather}C_{rr}=C_{\theta\theta}=1, \quad C_{r\theta}=C_{\theta z} = 0, \quad C_{rz}=WiU_z', \quad C_{zz}=1+2Wi^2 U_z'^2, \end{gather}

where the symbol prime $'$ denotes the derivative with respect to $r$.

In this work we will focus on extending the centre-mode instability as discovered by Garg et al. (Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018). Till now, only the axisymmetric mode in viscoelastic pipe flows have been found to be linearly unstable in the literature (Garg et al. Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018; Chaudhary et al. Reference Chaudhary, Garg, Subramanian and Shankar2021; Zhang Reference Zhang2021); therefore, it is legitimate to focus on the weakly nonlinear development of the axisymmetric mode in this work. For the axisymmetric viscoelastic pipe flow, the component-wise governing equations can be obtained from (2.3) as

(2.5a)\begin{gather} \frac{\partial u_{r}}{\partial r} + \frac{u_{r}}{r} +\frac{\partial u_{z}}{\partial z} = 0, \end{gather}
(2.5b)\begin{gather}\frac{\partial u_{r}}{\partial t} + U_{z} \frac{\partial u_{r}}{\partial z} + N_{u_{r}} ={-} \frac{\partial p}{\partial r} + \frac{\beta}{Re} \left( \nabla^{2}u_{r} - \frac{u_{r}}{r^{2}} \right) + \frac{1-\beta}{ReWi} \left( \frac{\partial c_{rr}}{\partial r} + \frac{c_{rr}}{r} - \frac{c_{\theta \theta}}{r} + \frac{\partial c_{rz}}{\partial z} \right), \end{gather}
(2.5c)\begin{gather}\frac{\partial u_{z}}{\partial t} + U_{z} \frac{\partial u_{z}}{\partial z} + u_{r} U_{z}' + N_{u_{z}} ={-} \frac{\partial p}{\partial z} + \frac{\beta}{Re} \left( \nabla^{2}u_{z} \right) + \frac{1-\beta}{ReWi} \left( \frac{\partial c_{rz}}{\partial r} + \frac{c_{rz}}{r} + \frac{\partial c_{zz}}{\partial z} \right), \end{gather}
(2.5d)\begin{gather}\frac{\partial c_{rr}}{\partial t} + u_{r} C_{rr}' + U_{z} \frac{\partial c_{rr}}{\partial z} - 2 C_{rr}\frac{\partial u_{r}}{\partial r} - 2 C_{rz}\frac{\partial u_{r}}{\partial z} + N_{c_{rr}} ={-} \frac{c_{rr}}{Wi}, \end{gather}
(2.5e)\begin{gather}\frac{\partial c_{rz}}{\partial t} + u_{r} C_{rz}' + U_{z} \frac{\partial c_{rz}}{\partial z} - c_{rr} U_{z}' - C_{rr}\frac{\partial u_{z}}{\partial r} - C_{rz}\frac{\partial u_{z}}{\partial z} - C_{rz}\frac{\partial u_{r}}{\partial r} - C_{zz}\frac{\partial u_{r}}{\partial z} + N_{c_{rz}} ={-} \frac{c_{rz}}{Wi}, \end{gather}
(2.5f)\begin{gather}\frac{\partial c_{\theta \theta}}{\partial t} + u_{r} C_{\theta \theta}' + U_{z} \frac{\partial c_{\theta \theta}}{\partial z} - 2 C_{\theta \theta}\frac{u_{r}}{r} + N_{c_{\theta \theta}} ={-} \frac{c_{\theta \theta}}{Wi}, \end{gather}
(2.5g)\begin{gather}\frac{\partial c_{zz}}{\partial t} + u_{r} C_{zz}' + U_{z} \frac{\partial c_{zz}}{\partial z} - 2 c_{rz} U_{z}' - 2 C_{rz}\frac{\partial u_{z}}{\partial r} - 2 C_{zz}\frac{\partial u_{z}}{\partial z} + N_{c_{zz}} ={-} \frac{c_{zz}}{Wi}, \end{gather}

where $\nabla ^{2}={\partial ^2}/{\partial r^2} + ({1}/{r})({ \partial }/{\partial r}) + {\partial ^2}/{\partial z^2}$ and the nonlinear terms are given as

(2.6a)\begin{gather} N_{u_{r}} = u_{r} \frac{\partial u_{r}}{\partial r} + u_{z} \frac{\partial u_{r}}{\partial z}, \quad N_{u_{z}} = u_{r} \frac{\partial u_{z}}{\partial r} + u_{z} \frac{\partial u_{z}}{\partial z}, \end{gather}
(2.6b)\begin{gather}N_{c_{rr}} = u_{r} \frac{\partial c_{rr}}{\partial r} + u_{z} \frac{\partial c_{rr}}{\partial z} - 2 c_{rr} \frac{\partial u_{r}}{\partial r} - 2 c_{rz} \frac{\partial u_{r}}{\partial z}, \end{gather}
(2.6c)\begin{gather}N_{c_{rz}} = u_{r} \frac{\partial c_{rz}}{\partial r} + u_{z} \frac{\partial c_{rz}}{\partial z} - c_{rr} \frac{\partial u_{z}}{\partial r} - c_{rz} \frac{\partial u_{z}}{\partial z} - c_{rz} \frac{\partial u_{r}}{\partial r} - c_{zz} \frac{\partial u_{r}}{\partial z}, \end{gather}
(2.6d)\begin{gather}N_{c_{\theta \theta}} = u_{r} \frac{\partial c_{\theta \theta}}{\partial r} + u_{z} \frac{\partial c_{\theta \theta}}{\partial z} - \frac{2 u_{r}}{r} c_{\theta \theta}, \end{gather}
(2.6e)\begin{gather}N_{c_{zz}} = u_{r} \frac{\partial c_{zz}}{\partial r} + u_{z} \frac{\partial c_{zz}}{\partial z} - 2 c_{rz} \frac{\partial u_{z}}{\partial r} - 2 c_{zz} \frac{\partial u_{z}}{\partial z}. \end{gather}

No-slip boundary conditions at the pipe wall for the velocity components are $u_{r}(1)=0$, $u_{z}(1)=0$, and we do not need to specify boundary conditions for the conformation tensor components at the wall. At the pipe axis, the conditions for both $\boldsymbol {u}$ and $\boldsymbol {c}$ are specified by the parity conditions as will be presented in § 3.

The above formulation in (2.5) having explicit pressure terms is referred to as the $up$-$c$ formulation hereafter. The traditional Reynolds decomposition has been applied in this formulation to the conformation tensor. However, this decomposition may jeopardise the positive definiteness of the conformation tensor $\hat {\boldsymbol {c}}$, leading to non-physical results. In the current work, we will also work with another formulation, i.e. the geometric decomposition proposed by Hameduddin et al. (Reference Hameduddin, Meneveau, Zaki and Gayme2018), Hameduddin, Gayme & Zaki (Reference Hameduddin, Gayme and Zaki2019) and Hameduddin & Zaki (Reference Hameduddin and Zaki2019) to decompose $\hat {\boldsymbol {c}}$. Another advantage of this formulation is that the elastic energy can be defined with a more physical significance (see the discussions on the inner product below). The perturbation $\boldsymbol {c}$ is expressed as $\boldsymbol {c} = \boldsymbol {F} \boldsymbol {\cdot } \boldsymbol {g} \boldsymbol {\cdot } \boldsymbol {F}^{\rm T}$, where $\boldsymbol {F}$ is a deformation gradient tensor and is given by the Cholesky decomposition of the conformation tensor of the base flow as $\boldsymbol {F} \cdot \boldsymbol {F}^{\rm T}=\boldsymbol {C}$; the perturbation tensor $\boldsymbol {g}$ can be interpreted as the polymer deformation relative to the mean configuration. By substituting the Cholesky decomposition into $\boldsymbol {c} = \boldsymbol {F} \boldsymbol {\cdot } \boldsymbol {g} \boldsymbol {\cdot } \boldsymbol {F}^{\rm T}$, the explicit expression of $\boldsymbol {c}=\boldsymbol {Q}\boldsymbol {g}$ can be obtained as (see also Zhang Reference Zhang2021)

(2.7)\begin{equation} \begin{pmatrix} c_{rr} \\ c_{rz} \\ c_{\theta\theta} \\ c_{zz} \end{pmatrix} = \begin{pmatrix} C_{rr} & 0 & 0 & 0\\ C_{rz} & S & 0 & 0 \\ 0 & 0 & C_{\theta\theta} & 0 \\ C_{rz}^{2}/C_{rr} & 2SC_{rz}/C_{rr} & 0 & S^2/C_{rr} \end{pmatrix} \begin{pmatrix} g_{rr} \\ g_{rz} \\ g_{\theta\theta} \\ g_{zz} \end{pmatrix}, \end{equation}

where the notation $S=\sqrt {C_{rr}C_{zz}-C_{rz}^2}$ is used. Besides, in the Navier–Stokes equations, we eliminate the pressure term using the continuity condition (so that the inner product of the variable array will bear the significance of the disturbance energy directly). To this end, we use the streamfunction $\phi$, which is related to the velocity components as $u_{r} = - ({1}/{r})({\partial \phi }/{\partial z})$, $u_{z} = ({1}/{r})({\partial \phi }/{\partial r})$. Instead of using $\phi$ directly for pipe flows, we will work with $\psi = \phi /r$ for the convenience of implementing the boundary conditions at $r = 0$ (Orlandi Reference Orlandi2012). The governing equation of $\psi$ reads as

(2.8)\begin{align} &\frac{\partial}{\partial t} \left( -\nabla^{2} \psi + \frac{\psi}{r^{2}} \right) \nonumber\\ &\quad = U_{z} \frac{\partial \nabla^{2} \psi}{\partial z} - \left( U_{z}'' - \frac{1}{r} U_{z}' + \frac{1}{r^2} U_{z} \right) \frac{\partial \psi}{\partial z}\nonumber\\ &\qquad - \frac{\beta}{Re} \left( \frac{\partial^{4}}{\partial r^{4}} + \frac{2}{r} \frac{\partial^{3}}{\partial r^{3}} - \frac{3}{r^{2}} \frac{\partial^{2}}{\partial r^{2}} + 2\frac{\partial^{4}}{\partial r^{2} \partial z^{2}} + \frac{3}{r^{3}} \frac{\partial}{\partial r} + \frac{2}{r} \frac{\partial^{3}}{\partial r \partial z^{2}} - \frac{3}{r^{4}} - \frac{2}{r^{2}} \frac{\partial^{2}}{\partial z^{2}} + \frac{\partial^{4}}{\partial z^{4}} \right) \psi\nonumber\\ &\qquad - \frac{1 - \beta}{ReWi} \left( F(g_{rr}) + F(g_{rz}) + \frac{C_{\theta\theta}}{r} \frac{\partial g_{\theta\theta}}{\partial z} + \frac{S^2}{C_{rr}} \frac{\partial^2 g_{zz}}{\partial r \partial z} + \left(\frac{S^2}{C_{rr}}\right)' \frac{\partial g_{zz}}{\partial z} \right) + N_{\psi}, \end{align}

where the nonlinear term is $N_{\psi } = {\partial N_{u_{z}}}/{\partial r} - {\partial N_{u_{r}}}/{\partial z}$; explicit expressions of the short notations $F(g_{rr})$ and $F(g_{rz})$, along with the governing equations of $(g_{rr},g_{rz},g_{\theta \theta },g_{zz})$, can be found in Appendix B; the boundary conditions for $\psi$ at the pipe wall are $\psi (1)=0$, ${\partial \psi }/{\partial r}|(1)=0$. This new equation system involving $\psi$ and $\boldsymbol {g}$ is denoted as the $\psi$-$g$ formulation.

A general compact form of the equation systems in the two formulations is

(2.9)\begin{equation} \boldsymbol{M} \frac{\partial \boldsymbol{\gamma}}{\partial t} = \boldsymbol{L} \boldsymbol{\gamma} + \boldsymbol{N}, \quad \text{or} \quad \left( \boldsymbol{M} \frac{\partial }{\partial t} - \boldsymbol{L} \right) \boldsymbol{\gamma} = \boldsymbol{N}, \end{equation}

where $\boldsymbol {\gamma }=(u_{r},u_{z},p,c_{rr},c_{rz},c_{\theta \theta }, c_{zz})^\textrm {T}$ in the $up$-$c$ formulation, and $\boldsymbol {\gamma }=(\psi,g_{rr},g_{rz},g_{\theta \theta }, g_{zz})^\textrm {T}$ in the $\psi$-$g$ formulation; $\boldsymbol {M}$ is the weight matrix, $\boldsymbol {L}$ the linear operator, $\boldsymbol {N}$ the nonlinear operator. The explicit expressions of these operators for the $up$-$c$ formulation are given in Appendix C.1 and those for the $\psi$-$g$ formulation can be derived in a similar way.

2.2. Multiple-scale expansion for the weakly nonlinear stability analysis

There are different ways to conduct the multiple-scale expansion of disturbances and equations. The specific expansion scheme below follows our previous work Zhang (Reference Zhang2016) (which followed and adapted the methods of Stewartson & Stuart Reference Stewartson and Stuart1971; Fujimura Reference Fujimura1989) and expands the time $t$, the spatial coordinate in the streamwise direction $z$, the disturbance $\boldsymbol {\gamma }$ and the governing parameter $Re$ as a function of a small dimensionless parameter $\epsilon$,

(2.10a)\begin{gather} \frac{\partial}{\partial t} = \frac{\partial}{\partial t_{0}} + \epsilon \frac{\partial}{\partial t_{1}} + \epsilon^{2} \frac{\partial}{\partial t_{2}} + O(\epsilon^{3}), \quad \frac{\partial}{\partial z} = \frac{\partial}{\partial z_{0}} + \epsilon \frac{\partial}{\partial z_{1}} + O(\epsilon^{2}), \end{gather}
(2.10b)\begin{gather}\boldsymbol{\gamma} = \epsilon \boldsymbol{\gamma}_{1} + \epsilon^2 \boldsymbol{\gamma}_{2} + \epsilon^3 \boldsymbol{\gamma}_{3} + O(\epsilon^{4}), \quad Re = Re_{c} + \epsilon^{2} + O(\epsilon^{4}). \end{gather}

There are three parameters $Re$, $Wi$ and $\beta$, while only $Re$ is perturbed around its linear critical value $Re_{c}$ at a given combination of $\beta$ and $Wi$. The expansion form of $Re$ is due to the parabolic shape of neutral curves around $Re_{c}$, originating from Stuart (Reference Stuart1960) in a weakly nonlinear analysis of Newtonian plane Poiseuille flows. Specifically, the small quantity $\epsilon$ can be considered as a measure of the distance between $Re$ and the linear critical condition $Re_c$. It should be mentioned that the definition of the expansion parameter $\epsilon$ is not unique and rescaling of the present $\epsilon = \sqrt {Re-Re_c}$ does not affect the Landau coefficient $a_3$ (derived in (2.24) to follow) as long as the computation is restricted to the linear critical conditions. An example illustrating this point can be found in table 1 in Zhang (Reference Zhang2016), where the Landau coefficient obtained from the weakly nonlinear analysis of Newtonian plane Poiseuille flow agrees well with that in Fujimura (Reference Fujimura1989) where a different expansion scheme of $Re$ was used. Since there are temporal and spatial derivatives $\partial /\partial t$ and $\partial /\partial z$ in operators as introduced in (2.9), these operators also need to be expanded as series of $\epsilon$,

(2.11a)\begin{gather} \boldsymbol{M} = \boldsymbol{M}_{0} + \epsilon \boldsymbol{M}_{1} + \epsilon^{2} \boldsymbol{M}_{2} + O(\epsilon^{3}), \quad \boldsymbol{L} = \boldsymbol{L}_{0} + \epsilon \boldsymbol{L}_{1} + \epsilon^{2} \boldsymbol{L}_{2} + O(\epsilon^{3}), \end{gather}
(2.11b)\begin{gather}\boldsymbol{N} = \epsilon^{2} \boldsymbol{N}_{2} + \epsilon^{3} \boldsymbol{N}_{3} + O(\epsilon^{4}). \end{gather}

The explicit expression of these subscale operators are provided in Appendix C.2 for the $up$-$c$ formulation.

Next, we substitute (2.10) and (2.11) into (2.9), and collect the terms at the same order of $\epsilon$. This step gives the equation at order $\epsilon$,

(2.12)\begin{equation} \left( \boldsymbol{M}_{0} \frac{\partial}{\partial t_{0}} - \boldsymbol{L}_{0} \right) \boldsymbol{\gamma}_{1} = 0, \end{equation}

the equation at order $\epsilon ^{2}$,

(2.13)\begin{equation} \left( \boldsymbol{M}_{0} \frac{\partial}{\partial t_{0}} - \boldsymbol{L}_{0} \right) \boldsymbol{\gamma}_{2} = \left( \boldsymbol{L}_{1} - \boldsymbol{M}_{1} \frac{\partial}{\partial t_{0}} - \boldsymbol{M}_{0} \frac{\partial}{\partial t_{1}} \right) \boldsymbol{\gamma}_{1} + \boldsymbol{N}_{2}, \end{equation}

and the equation at order $\epsilon ^{3}$,

(2.14)\begin{align} \left( \boldsymbol{M}_{0} \frac{\partial}{\partial t_{0}} - \boldsymbol{L}_{0} \right) \boldsymbol{\gamma}_{3} &= \left( \boldsymbol{L}_{2} - \boldsymbol{M}_{2} \frac{\partial}{\partial t_{0}} - \boldsymbol{M}_{1} \frac{\partial}{\partial t_{1}} - \boldsymbol{M}_{0} \frac{\partial}{\partial t_{2}} \right) \boldsymbol{\gamma}_{1} \nonumber\\ &\quad + \left( \boldsymbol{L}_{1} - \boldsymbol{M}_{1} \frac{\partial}{\partial t_{0}} - \boldsymbol{M}_{0} \frac{\partial}{\partial t_{1}} \right) \boldsymbol{\gamma}_{2} + \boldsymbol{N}_{3}. \end{align}

At order $\epsilon$, (2.12) is exactly the linearised equation in the framework of linear stability analysis. Its solution can be assumed to take the wave-like form of

(2.15)\begin{equation} \boldsymbol{\gamma}_{1}(z_{0},z_{1},r,t_{0},t_{1},t_{2};Re) = A(z_{1},t_{1},t_{2})\tilde{\boldsymbol{\gamma}}_1(r)\exp({\textrm{i} \alpha z_{0} + \mu t_{0}}) + \text{c.c.}, \end{equation}

and by substituting (2.15) into (2.13), then the second-order solution can be expressed in the form of

(2.16)\begin{equation} \boldsymbol{\gamma}_{2} = A^*A\tilde{\boldsymbol{\gamma}}_{20}(r) \!+\! \frac{ \partial A}{ \partial z_1}\tilde{\boldsymbol{\gamma}}_{21}(r) \exp({\textrm{i} \alpha z_{0} \!+\! \mu t_{0}}) + \text{c.c.} + A^2 \tilde{\boldsymbol{\gamma}}_{22}(r) \exp({2\textrm{i} \alpha z_{0} \!+\! 2\mu t_{0}}) \!+\! \text{c.c.}, \end{equation}

where c.c. stands for the complex conjugate of the preceding term; $\textrm {i}=\sqrt {-1}$ is the imaginary unit; $\alpha$ is the wavenumber; $\mu =-\textrm {i}\omega$ and $\tilde {\boldsymbol {\gamma }}_1$ are the complex eigenvalue and eigenfunction, respectively, of the eigenvalue problem $\mu \tilde {\boldsymbol {M}}_{0} \tilde {\boldsymbol {\gamma }}_1= \tilde {\boldsymbol {L}}_{0} \tilde {\boldsymbol {\gamma }}_1$ resulting from (2.12) (the symbol tilde $\tilde{\phantom{0}}$ denotes matrices/variables in the spectral space); $A$ is the complex amplitude of the disturbance (the physical meaning of $A$ depends on the normalization of $\tilde{\boldsymbol{\gamma}}_1$ as will be discussed in the end of this section), and its value cannot be determined from the linear equation as the linear problem can be arbitrarily scaled. Instead, the evolution equation of $A$ (i.e. Ginzburg–Landau equation) can be obtained by applying solvability conditions on (2.13) and (2.14) at orders $\epsilon ^{2}$ and $\epsilon ^{3}$, respectively. The aim is to remove the secular terms from the inhomogeneous terms (Bender & Orszag Reference Bender and Orszag1999) and this can be done via adjoint variables.

To formulate the adjoint problem corresponding to the linear eigenvalue problem (2.12), we introduce an inner product defined as

(2.17)\begin{align} \langle \boldsymbol{a},\boldsymbol{b} \rangle &= \frac{1}{2 L_p} \int_0^{L_p} \int_{0}^{1} \boldsymbol{a} \boldsymbol{\cdot} \boldsymbol{W} \boldsymbol{b} \, r\,\textrm{d}r \, \textrm{d}z_{0} = \frac{A_{\boldsymbol{a}}A_{\boldsymbol{b}}^*}{2} \int_{0}^{1} \tilde{\boldsymbol{a}} \boldsymbol{\cdot} \boldsymbol{W} {\tilde{\boldsymbol{b}}}^{*} \, r\,\textrm{d}r + \text{c.c.}\nonumber\\ & = A_{\boldsymbol{a}}A_{\boldsymbol{b}}^*{\langle \tilde{\boldsymbol{a}}, \tilde{\boldsymbol{b}} \rangle}_{s} + \text{c.c.}, \end{align}

where $L_p=2{\rm \pi} /\alpha$ is the wavelength; $\boldsymbol {a}$ and $\boldsymbol {b}$ are real-valued vectors in the same form as the linear solution in (2.15) ($A_{\boldsymbol {a}}$ and $A_{\boldsymbol {b}}$ are the corresponding complex amplitudes); $\boldsymbol {W}$ is a real coefficient matrix; ${\langle \tilde {\boldsymbol {a}}, \tilde {\boldsymbol {b}} \rangle }_{s}$ can be seen as the inner product defined in the spectral space. Under this definition, the inner product $\langle \boldsymbol {\gamma }_{1},\boldsymbol {M}_{0} \boldsymbol {\gamma }_{1} \rangle$ in the $\psi$-$g$ formulation represents the total energy (kinetic energy plus elastic energy) of the first-order disturbance with the coefficient matrix being specified as $\boldsymbol {W} = \textrm {diag}(I,wI,2wI,wI,wI)$, where $w = (1-\beta )/Re/Wi$. We illustrate this point as follows.

For the kinetic energy part in $\langle \boldsymbol {\gamma }_{1},\boldsymbol {M}_{0} \boldsymbol {\gamma }_{1} \rangle$, we single out $\langle \psi _{1}, (-\nabla ^2+{1}/{r^2}) \psi _{1} \rangle$ (where the element in $\boldsymbol {M}_{0}$ for $\psi$ is $-\nabla ^2+{1}/{r^2}$, see (2.8)) and one is able to get

(2.18)\begin{equation} E_{k} = \frac{1}{2 L_p} \int_0^{L_p} \int_0^1 \left(u_{r,1}^2 + u_{z,1}^2 \right) r\,\textrm{d}r \, \textrm{d}z_{0} = \left\langle \psi_{1}, \ \left(-\nabla^2+\frac{1}{r^2}\right) \psi_{1} \right\rangle, \end{equation}

which calculates the disturbance kinetic energy. The elastic energy part in $\langle \boldsymbol {\gamma }_1,\boldsymbol {M}_{0}\boldsymbol {\gamma }_1 \rangle$ can be directly written as

(2.19)\begin{equation} E_e =\frac{1}{2 L_p} \int_0^{L_p} \int_{0}^{1} \left( \frac{1-\beta}{Re Wi}(g_{rr,1}^{2}+2g_{rz,1}^{2}+g_{\theta\theta, 1}^{2}+g_{zz,1}^2) \right) r\,\textrm{d}r \, \textrm{d}z_{0}, \end{equation}

where we have used the same concept of geodesic distance as in the calculation of elastic energy proposed by Hameduddin et al. (Reference Hameduddin, Meneveau, Zaki and Gayme2018,  Reference Hameduddin, Gayme and Zaki2019). Thus, the inner product defined with the coefficient matrix $\boldsymbol {W}$ as above bears the physical meaning of total energy. Nevertheless, we find that the results in § 4 (especially the value of the Landau coefficient $a_3$ in (2.25)) do not depend on the specific choice of the coefficient matrix. Therefore, in the $up$-$c$ formulation, for simplicity, we use an identity matrix as the coefficient matrix, resulting in

(2.20)\begin{equation} \langle \boldsymbol{\gamma}_1,\boldsymbol{M}_{0}\boldsymbol{\gamma}_1 \rangle =\frac{1}{2 L_p} \int_0^{L_p} \int_{0}^{1} ( u_{r,1}^{2}+u_{z,1}^{2}+c_{rr,1}^{2}+c_{rz,1}^{2}+c_{\theta\theta,1}^{2}+c_{zz,1}^{2} ) \, r\,\textrm{d}r \, \textrm{d}z_{0}. \end{equation}

With the inner product defined, through integration by parts, we can obtain the adjoint equation (Luchini & Bottaro Reference Luchini and Bottaro2014)

(2.21)\begin{equation} \left\langle \left( \boldsymbol{M}_{0} \frac{\partial}{\partial t_{0}} - \boldsymbol{L}_{0} \right) \boldsymbol{\gamma}_{1}, \boldsymbol{\gamma}_{1}^{{{\dagger}}} \right\rangle = \left\langle \boldsymbol{\gamma}_{1}, \left( \boldsymbol{M}_{0}^{{{\dagger}}} \frac{\partial}{\partial t_{0}} - \boldsymbol{L}_{0}^{{{\dagger}}} \right) \boldsymbol{\gamma}_{1}^{{{\dagger}}} \right\rangle, \quad \!\rightarrow\! \quad \left( \boldsymbol{M}_{0}^{{{\dagger}}} \frac{\partial}{\partial t_{0}} - \boldsymbol{L}_{0}^{{{\dagger}}} \right) \boldsymbol{\gamma}_{1}^{{{\dagger}}} \!=\! 0, \end{equation}

where the symbol dagger ${\dagger}$ denotes adjoint variables and operators; explicit expressions of the adjoint operators $\boldsymbol {M}_{0}^{{\dagger} }$ and $\boldsymbol {L}_{0}^{{\dagger} }$ can be found in Appendix C.3. The boundary conditions for the adjoint problem are the same as those for the direct problem: at the pipe wall, the adjoint velocity components $u_{r,1}^{{\dagger} }(1)=0$, $u_{z,1}^{{\dagger} }(1)=0$, and there is no need to specify boundary conditions for the adjoint conformation tensor components; at the pipe axis, parity conditions are applied to $\boldsymbol {u}_1^{{\dagger} }$ and $\boldsymbol {c}_1^{{\dagger} }$. As mentioned earlier, in the inhomogeneous equations (2.13) and (2.14), the secular terms – terms related to the wave component $\exp ({\textrm {i} \alpha z_{0} + \mu t_{0}})$ and its complex conjugate – on the right-hand side may cause the solution to blow up and, thus, the solvability conditions need to be enforced, dictating the projection of the inhomogeneous terms of (2.13) and (2.14) on the adjoint waves $\boldsymbol {\gamma }_{1}^{{\dagger} }$ to be zero, i.e.

(2.22a)\begin{gather} \left\langle \left( \boldsymbol{L}_{1} - \boldsymbol{M}_{1} \frac{\partial}{\partial t_{0}} - \boldsymbol{M}_{0} \frac{\partial}{\partial t_{1}} \right) \boldsymbol{\gamma}_{1} + \boldsymbol{N}_2, \boldsymbol{\gamma}_{1}^{{{\dagger}}} \right\rangle = 0, \end{gather}
(2.22b)\begin{gather}\left\langle \left( \boldsymbol{L}_{2} - \boldsymbol{M}_{2} \frac{\partial}{\partial t_{0}} - \boldsymbol{M}_{1} \frac{\partial}{\partial t_{1}} - \boldsymbol{M}_{0} \frac{\partial}{\partial t_{2}} \right) \boldsymbol{\gamma}_{1} + \left( \boldsymbol{L}_{1} \!-\! \boldsymbol{M}_{1} \frac{\partial}{\partial t_{0}} \!-\! \boldsymbol{M}_{0} \frac{\partial}{\partial t_{1}} \right) \boldsymbol{\gamma}_{2} \!+\! \boldsymbol{N}_3, \boldsymbol{\gamma}_{1}^{{{\dagger}}} \right\rangle = 0. \end{gather}

By singling out the terms related to the wave component $\exp ({\textrm {i} \alpha z_{0} + \mu t_{0}})$ and its complex conjugate, and expressing the above projection in the spectral space, the following equation with group velocity $c_{g}$ can be derived from (2.22a):

(2.23)\begin{equation} \frac{\partial A}{\partial t_{1}} + c_{g}\frac{\partial A}{\partial z_{1}} = 0 \quad \text{with} \ c_{g} ={-}\frac{ {\langle (\tilde{\boldsymbol{L}}_{1}^{{\circ}} - \mu_{c}\tilde{\boldsymbol{M}}_{1}^{{\circ}} )\tilde{\boldsymbol{\gamma}}_{1}, \tilde{\boldsymbol{\gamma}}_{1}^{{{\dagger}}} \rangle}_{s}}{{\langle \tilde{\boldsymbol{M}}_{0}\tilde{\boldsymbol{\gamma}}_1, \tilde{\boldsymbol{\gamma}}_{1}^{{{\dagger}}} \rangle}_{s}}. \end{equation}

Here $\mu _{c}=\mu _{c,r} + \textrm {i} \mu _{c,i}$ is the eigenvalue at the critical condition $(Re_c, \alpha _{c})$, i.e. $\mu _{c,r}$ (bearing the meaning of growth rate) should be zero theoretically or a very small value numerically. The superscript $^\circ$ means that the $\partial /\partial z_{1}$-related derivatives and the $A$-related coefficients have been moved out of the corresponding operators (see Appendix C.2 for details in the $up$-$c$ formulation). A Ginzburg–Landau equation (GLE) governing the evolution of the complex amplitude $A$ can be derived from (2.22b),

(2.24)\begin{equation} \frac{\partial A}{\partial t_{2}} = a_{1}A + a_{2}\frac{\partial^{2} A}{\partial z_{1}^{2}} + a_{3}{|A|}^{2}A, \end{equation}

where the three coefficients can be expressed as

(2.25a)\begin{gather} a_{1} = \frac{{\langle \tilde{\boldsymbol{L}}_{2Re}\tilde{\boldsymbol{\gamma}}_1, \tilde{\boldsymbol{\gamma}}_{1}^{{{\dagger}}} \rangle}_{s}}{{\langle \tilde{\boldsymbol{M}}_{0}\tilde{\boldsymbol{\gamma}}_1, \tilde{\boldsymbol{\gamma}}_{1}^{{{\dagger}}} \rangle}_{s}}, \end{gather}
(2.25b)\begin{gather}a_{2} = \frac{{\langle (\tilde{\boldsymbol{L}}_{2}^{{\circ}} - \mu_{c}\tilde{\boldsymbol{M}}_{2}^{{\circ}} - c_{g} \tilde{\boldsymbol{M}}_{1}^{{\circ}})\tilde{\boldsymbol{\gamma}}_{1} + (\tilde{\boldsymbol{L}}_{1}^{{\circ}} - \mu_{c}\tilde{\boldsymbol{M}}_{1}^{{\circ}} - c_{g}\tilde{\boldsymbol{M}}_{0}^{{\circ}})\tilde{\boldsymbol{\gamma}}_{21}, \tilde{\boldsymbol{\gamma}}_{1}^{{{\dagger}}} \rangle}_{s}}{{\langle \tilde{\boldsymbol{M}}_{0}\tilde{\boldsymbol{\gamma}}_1, \tilde{\boldsymbol{\gamma}}_{1}^{{{\dagger}}} \rangle}_{s}}, \end{gather}
(2.25c)\begin{gather}a_{3} = \frac{{\langle \tilde{\boldsymbol{N}}_{3}^{{\circ}}, \tilde{\boldsymbol{\gamma}}_{1}^{{{\dagger}}} \rangle}_{s}}{{\langle \tilde{\boldsymbol{M}}_{0}\tilde{\boldsymbol{\gamma}}_1, \tilde{\boldsymbol{\gamma}}_{1}^{{{\dagger}}} \rangle}_{s}}. \end{gather}

Among the above coefficients, the complex Landau coefficient $a_{3}=a_{3,r}+\textrm {i} a_{3,i}$ is of particular interest in this study as the sign of its real part can reveal the bifurcation type of the flow studied (i.e. the bifurcation is supercritical if $a_{3,r}<0$ and subcritical if $a_{3,r}>0$) as long as the amplitudes of the higher-order coefficients do not increase significantly fast with right signs. The present GLE is truncated to the third order (up to $a_3$). The expansion (even to higher orders) is known to be well convergent a priori as long as the calculation is performed near the linear critical condition (Herbert Reference Herbert1980; Sen & Venkateswarlu Reference Sen and Venkateswarlu1983), which will be followed here. Besides, in the present study, we have also used DNS to determine the value of $a_3$ and the DNS results validate the bifurcation type (see Appendix D), suggesting that the present truncation of GLE up to third order is reasonable. We mention in passing that even though we focus solely on the significance of $a_3$ in this paper, the application of the GLE to our flow can be general (see Cross & Hohenberg Reference Cross and Hohenberg1993). Another issue to note is that, to avoid ambiguity, $\tilde {\boldsymbol {\gamma }}_1$ should be normalized so that the Landau coefficient $a_3$ can be uniquely determined. In the current work, we normalize $\tilde {\boldsymbol {\gamma }}_1$ so that

(2.26)\begin{equation} \sqrt{\frac{1}{2} \int_{0}^{1} \left( \left(|\tilde{u}_{r,1}|^2 + |\tilde{u}_{z,1}|^2 \right) + \frac{1-\beta}{Re Wi} \left( |\tilde{g}_{rr,1}|^{2}+2|\tilde{g}_{rz,1}|^{2}+|\tilde{g}_{\theta\theta,1}|^{2}+|\tilde{g}_{zz,1}|^2\right) \right) r\,\textrm{d}r} = 1. \end{equation}

With this normalization, the physical meaning of $|A|$ in (2.15) is the square root of the total energy (kinetic energy plus elastic energy) in both formulations (conversion is needed where different variables are used).

3. Numerical method

A spectral collocation method has been adopted to solve the linear equation (2.12) at order $\epsilon$, (2.13) at order $\epsilon ^{2}$ for the second-order solution $\boldsymbol {\gamma }_{2}$ and the adjoint equation (2.21), as well as to evaluate the Landau coefficients in (2.25). To enforce the no-slip boundary condition for velocity at the pipe wall, we remove the corresponding rows and columns of the node $r=1$ in the matrices. We apply no boundary conditions for the conformation tensor $\boldsymbol {c}$ or $\boldsymbol {g}$ at the pipe wall. As for the singularity problem at $r=0$, we avoid placing a node at $r=0$ (Mohseni & Colonius Reference Mohseni and Colonius2000) and the derivative matrices are constructed with an even–odd property following Trefethen (Reference Trefethen2000). In our axisymmetric pipe flows, in the $up$-$c$ formulation, $u_{r}$, $c_{rz}$ are odd functions and $u_z$, $p$, $c_{rr}$, $c_{\theta \theta }$, $c_{zz}$ are even; in the $\psi$-$g$ formulation, $\psi$, $g_{rz}$ are odd and $g_{rr}$, $g_{\theta \theta }$, $g_{zz}$ are even. The same parity conditions are applied to the corresponding adjoint variables. The present weakly nonlinear code is adapted from the code used by Zhang (Reference Zhang2016), which has been validated by comparing with the results in the literature. The adaptation is mainly for implementing the constitutive governing equations for polymers, while the main framework of the code remains unchanged.

4. Results and discussions

4.1. Validation

The linear instability of viscoelastic pipe flow of Oldroyd-B fluids was firstly reported by Garg et al. (Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018). A typical parameter setting at which they found a single unstable mode is $Re=800$, $\alpha =1$, $Wi=65$ and $\beta =0.65$. We solved the linear eigenvalue problem with the same parameters and validated our code by comparing our results with theirs. Firstly, the convergence of the eigenvalue of the unstable mode has been examined by varying the total number of interior nodes $N$ within $0< r<1$ in the spectral collocation method. The results are listed in table 1 for both the $up$-$c$ formulation and the $\psi$-$g$ formulation. It is clear that, in the $up$-$c$ formulation, $N=125$ gives well-converged results (up to 9–11 decimal numbers), while in the $\psi$-$g$ formulation, it is more difficult to converge to the same accuracy as the $up$-$c$ formulation. Secondly, we compare the present eigenspectra of Oldroyd-B pipe flow with that obtained by Garg et al. (Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018) in figure 1(a). An overall good agreement can be observed regarding the locations of the continuous spectra and the discrete modes (the data labelled with a blue dot were manually digitized from figure 1 in Garg et al. Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018). Moreover, we validate our result of linear stability analysis of UCM pipe flow by comparing with that reported in Chaudhary et al. (Reference Chaudhary, Garg, Subramanian and Shankar2021). As shown in figure 1(b), the variations of the growth rate $c_i$ and phase speed $c_r$ of the least stable mode with elasticity number $E$ at $Re=6000, \alpha =2$ agree very well with theirs.

Figure 1.

(a) Comparison of the eigenspectra of viscoelastic pipe flow of Oldroyd-B fluids at $Re=800$, $\alpha =1$, $Wi=65$, $\beta =0.65$; the inset is an enlarged view near the unstable mode. (b) Effect of $E$ on the eigenvalue of the least stable centre mode ($c=\omega /\alpha$ is used to facilitate the comparison and $c_i$ is the imaginary part and $c_r$ the real part) for UCM fluids at $Re=6000$ and $\alpha =2$.

Table 1.

Convergence of the eigenvalue of the single unstable mode in two different formulations ($up$-$c$ and $\psi$-$g$) of the viscoelastic pipe flow at $Re=800$, $\alpha =1$, $Wi=65$ and $\beta =0.65$.

In addition to the convergence check of the unstable mode in the linear eigenvalue problem, we have also examined the convergence of the group velocity $c_g$ in (2.23) and the coefficients $a_1$, $a_2$, $a_3$ in (2.25). To calculate these coefficients, the critical parameter needs to be determined first. As will be presented in the next subsection, the neutral curve (see figure 3(a) for an example at $Wi=65$, $\beta =0.65$) is in the form of a closed loop and we only consider the critical condition at the left end where transition from laminar to turbulent happens upon increasing $Re$. At $Wi=65$, $\beta =0.65$, the left-end critical parameters are $Re_{c}=265.572335$, $\alpha _{c}=0.515525$. The convergence of the critical eigenvalue $\mu _{c}$, the group velocity $c_{g}$ and the three coefficients $a_{1}$, $a_{2}$ and $a_{3}$ in the GLE are examined using the $\psi$-$g$ formulation as shown in table 2. We can see that the resolution $N=400$ is high enough to get well-converged results while $N=200$ also gives results with little loss of accuracy. In addition, the results obtained from the $up$-$c$ formulation are shown in table 3. It is clear that $c_{g}$, $a_{1}$, $a_{2}$ and $a_{3}$ agree well with those obtained from the $\psi$-$g$ formulation where the first three to six non-zero digits are identical. In this case, the real part of the Landau coefficient $a_3$ is negative, indicating a supercritical bifurcation at these parameters. There is a residual of the order of $10^{-6} \sim 10^{-7}$ in the imaginary part of the group velocity $c_g$ which is supposed to be a real number, and the reason may be attributed to numerical discretization errors. Such kind of residual has also been reported in weakly nonlinear analyses of other fluid systems (for example, the residual $O(10^{-7})$ in Kolyshkin & Ghidaoui (Reference Kolyshkin and Ghidaoui2003) and the residual $O(10^{-5})$ in Gao et al. Reference Gao, Sergent, Podvin, Xin, Le Quéré and Tuckerman2013; Zhang Reference Zhang2016). By comparing the inner product $\langle \boldsymbol {\gamma }_{1},\boldsymbol {M}_{0} \boldsymbol {\gamma }_{1} \rangle$ in the $up$-$c$ (2.20) and $\psi$-$g$ (2.18), (2.19) formulations, one can see that the $\psi$-$g$ formulation has clear physical meanings whereas the $up$-$c$ formulation does not bear any physical significance. Nevertheless, the results of the Landau coefficient are consistent. We have additionally formulated the problem (details not shown in this paper) with variables $u_{r}$, $u_{z}$, $p$, $g_{rr}$, $g_{rz}$, $g_{\theta \theta }$ and $g_{zz}$. All these three formulations give the same results, indicating that the specific definition of the inner product does not affect the values of Landau coefficients in the present weakly nonlinear stability analysis. In the following result section, the $up$-$c$ formulation will be used to generate most of the results.

Table 2.

Convergence of the Landau coefficients of the viscoelastic pipe flow of Oldroyd-B fluids at $Wi=65$, $\beta =0.65$, ${Re}_{c}=265.572335$ and $\alpha _{c}=0.515525$ using the $\psi$-$g$ formulation. Note that $a_{1}^* = a_{1} \times {10}^{5}$.

Table 3.

Convergence of the Landau coefficients of the viscoelastic pipe flow of Oldroyd-B fluids at $Wi=65$, $\beta =0.65$, ${Re}_{c}=265.572335$ and $\alpha _{c}=0.515525$ using the $up$-$c$ formulation. Note that $a_{1}^* = a_{1} \times {10}^{5}$.

Next, we validate the group velocity $c_{g}$ and the coefficient $a_{2}$ (listed in tables 2 and 3) in another way by examining the dispersion relation between the eigenvalue $\mu = \mu _{r} + \textrm {i} \mu _{i}$ and the wavenumber $\alpha$ near the critical wavenumber $\alpha _{c}$. Firstly, the growth rate $\mu _{r}(\alpha )$ can be approximated by a parabola function of $\alpha$ around the critical parameter $(Re_{c},\alpha _{c})$ (Huerre & Rossi Reference Huerre and Rossi1998). In the present case, figure 2(a) indeed presents such a parabola and the fitted function is shown in (4.1a) below. Secondly, the relationship $c_{g}=-\partial \mu _{i}/ \partial \alpha$ holds at $\alpha _{c}$. To verify this, $\mu _{i}$ as a function of $\alpha$ is plotted in figure 2(b) and a second-order polynomial is again used to fit the data (even though the curve looks like a straightline); evaluating $-\partial \mu _{i}/ \partial \alpha$ at $\alpha _{c}$ from function (4.1b) gives $1.012771$ which is quite close to the result $c_{g} = 1.012781 - 0.000001\textrm {i}$ in table 2 and $c_{g} = 1.012787 - 0.000003\textrm {i}$ in table 3. The fitted functions of $\mu_r$ and $\mu_i$ are

(4.1a)\begin{gather} \mu_{r,fit} ={-}0.033015 \alpha^{2} + 0.034052 \alpha - 0.008780, \end{gather}
(4.1b)\begin{gather}\mu_{i,fit} = 0.015435 \alpha^{2} - 1.028685 \alpha - 0.014776. \end{gather}

Thirdly, following Stewartson & Stuart (Reference Stewartson and Stuart1971), the coefficient $a_{2}$ is related to the eigenvalue $\mu$ in the form of

(4.2)\begin{equation} a_{2} ={-} \frac{1}{2} \frac{\partial^2 \mu_{r}}{\partial \alpha^{2}} - \frac{1}{2} \frac{\partial^2 \mu_{i}}{\partial \alpha^{2}}\textrm{i}. \end{equation}

Substitution of the fitted functions (4.1) into (4.2) gives $a_{2,fit}=0.033015-0.015435\textrm {i}$, in good agreement with $a_{2}=0.033028-0.015449\textrm {i}$ listed in table 2 and $a_{2}=0.033033-0.015465\textrm {i}$ in table 3. In Appendix D we additionally use DNS to qualitatively verify our calculations (in terms of the linear coefficient $a_1$ and Landau coefficient $a_3$) to confirm that the bifurcation type is correct.

Figure 2.

Dispersion relation at $Re_{c}=265.572335$, $\alpha _{c}=0.515525$, $Wi=65$ and $\beta =0.65$: (a) $\mu _r$ and (b) $\mu _i$ as functions of $\alpha$.

4.2. Neutral stability curves and linear instability of UCM pipe flows

Since the weakly nonlinear stability analysis is performed near the critical condition, it is necessary to calculate first the critical conditions and be aware of where the critical conditions are in the parameter space. Thus, in this section we present the neutral curves. In figure 3(a) a typical neutral curve at $\beta =0.65$ and $Wi=65$ is plotted, which appears in the form of a closed loop. In the study of viscoelastic channel flows by Chaudhary et al. (Reference Chaudhary, Garg, Shankar and Subramanian2019), the neutral curves also appear in the form of loops with the interior area being linearly unstable. We look at the critical conditions at the two ends in terms of $Re$. The critical Reynolds number and wavenumber are $Re_{c}=265.572335$, $\alpha _{c}=0.515525$ at the left end, and $Re_{c}=3553.575$, $\alpha _{c}=4.7337$ at the right end. Although they are both centre modes, the corresponding eigenfunctions are different, as plotted in figure 3(b). It appears that the velocity eigenfunctions corresponding to the right critical point (see the dashed and dotted curves) are more localized near the pipe axis than those of the left one. In the remaining part of the paper, we will restrict our discussion to the left-end critical point, as it is related to the transition from a laminar state to turbulence when we increase $Re$. Figure 3(c) shows neutral curves in the $(Re,\alpha )$ plane for varying elasticity numbers $E$ at a fixed viscosity ratio $\beta =0.65$ in a log-log plot, following the same format in Chaudhary et al. (Reference Chaudhary, Garg, Subramanian and Shankar2021). We present four dot–dashed curves corresponding to $E=0.02,0.04,0.12,0.22$ and superpose the data points for $E=0.12$ (blue square) extracted from Chaudhary et al. (Reference Chaudhary, Garg, Subramanian and Shankar2021) (see their figure 18a). Our results are in good agreement with theirs. Other cases of $E$ also show good agreement, but are not presented in order to not overload the figure. Similar to figure 3(a), the neutral curves can also be plotted at a constant $Wi$, which again exhibit in the form of a closed loop. In figure 3(d) the same neutral curves are displayed in a linear-linear plot. The purpose of showing these plots is to figure out where the critical conditions are, so that we can perform weakly nonlinear stability analysis around them.

Figure 3.

(a) Neutral curve in a $(\alpha,Re)$ plane with the colourbar for levels of growth rate and the area inside the loop is linear unstable. (b) Amplitudes of velocity eigenfunctions, where solid and dash–dotted curves are for the left-end critical point $Re_{c}=265.572335$, $\alpha _{c}=0.515525$ and dashed and dotted curves for the right-end critical point $Re_{c}=3553.575$, $\alpha _{c}=4.7337$. (c) Neutral curves in a log-log plot at a fixed viscosity ratio $\beta =0.65$ for various elasticity numbers $E$ and Weissenberg numbers $Wi$, where the discrete data points labelled by blue squares are manually digitized from figure 18(a) in Chaudhary et al. (Reference Chaudhary, Garg, Subramanian and Shankar2021). (d) Neutral curves in a linear-linear plot at various Weissenberg numbers $Wi$ for $\beta =0.65$.

For the left-end critical point, the corresponding adjoint eigenfunctions of the linear problem are shown in figures 4(a)–4(c). As the adjoint problem is linear and can be arbitrarily scaled, we normalize the adjoint eigenfunction so that the maximum amplitude of the adjoint axial velocity $\max (|u_{z,1}^{{\dagger} }|)=1$. With this normalisation, the adjoint eigenfunction can be uniquely determined and it is noted that this normalization process does not affect the value of the Landau coefficient $a_3$ (see (2.25c) ). In figure 4 we show the adjoint modes in the first row and the corresponding direct modes in the second row for a comparison. As shown in panel $(a)$, the adjoint velocity components $|u_{r,1}^{{\dagger} }|, |u_{z,1}^{{\dagger} }|$ obtained from the $up$-$c$ and $\psi$-$g$ formulations agree well with each other. On the other hand, figures 4(b) and 4(c) plot respectively the adjoint conformation tensor components $g_{rz,1}^{{\dagger} }$ vs $c_{rz,1}^{{\dagger} }$ and $g_{zz,1}^{{\dagger} }$ vs $c_{zz,1}^{{\dagger} }$ in the two formulations. We find that even though for the linear variables, we have $\boldsymbol {c}=\boldsymbol {Q}\boldsymbol {g}$ (as introduced in (2.7) and plotted in figure 4d–f), there is no simple relation between $\boldsymbol {c}^{{\dagger} }$ and $\boldsymbol {g}^{{\dagger} }$, probably due to the non-normality of $\boldsymbol {Q}$ and the linear operators. Still, the Landau coefficients obtained from these two different formulations are in good agreement as discussed in the validation subsection, suggesting that the solvability conditions described in (2.22) have been correctly applied with the aid of the adjoint variables.

Figure 4.

Adjoint and direct eigenfunctions at $Re_{c}=265.572335$, $\alpha _{c}=0.515525$, $Wi=65$ and $\beta =0.65$, corresponding to the left-end critical point in figure 3(a): (a) adjoint velocity components $\tilde {u}_{z,1}^{{\dagger} }$ and $\tilde {u}_{r,1}^{{\dagger} }$; (b) adjoint conformation tensor components $\tilde {g}_{rz,1}^{{\dagger} }$ and $\tilde {c}_{rz,1}^{{\dagger} }$; (c) adjoint conformation tensor components $\tilde {g}_{zz,1}^{{\dagger} }$ and $\tilde {c}_{zz,1}^{{\dagger} }$; (d) velocity components $\tilde {u}_{z,1}$ and $\tilde {u}_{r,1}$; (e) conformation tensor components $\tilde {c}_{rr,1}$ and $\tilde {c}_{rz,1}$; (f) conformation tensor components $\tilde {c}_{\theta \theta,1}$ and $\tilde {c}_{zz,1}$. The conformation tensor components $\tilde {c}_{rr,1}$, $\tilde {c}_{rz,1}$, $\tilde {c}_{\theta \theta,1}$, $\tilde {c}_{zz,1}$ in the $\psi$-$g$ formulation are transformed from $\tilde {g}_{rr,1}$, $\tilde {g}_{rz,1}$, $\tilde {g}_{\theta \theta,1}$, $\tilde {g}_{zz,1}$ according to (2.7).

Before ending this section, we mention that pipe flow of UCM fluids can be linearly unstable. The UCM fluids can be obtained from the Oldroyd-B fluids with $\beta =0$ (zero solvent viscosity) and, thus, the flow is dominated by inertia, polymer viscosity and elasticity. To the best of our knowledge, no previous works have reported linear instability in pipe flows of UCM fluids. After exploring a large parameter space of axisymmetric UCM flows covering $0.5<\alpha <3$, $100< Re<10000$ and $0< E<1$, Chaudhary et al. (Reference Chaudhary, Garg, Subramanian and Shankar2021) found stable modes in this flow. They concluded that solvent viscosity is important in constituting the linear instability in viscoelastic pipe flows along with the elasticity and fluid inertia. However, we find that the viscoelastic pipe flow of UCM fluids can be linearly unstable. Figure 5 shows the eigenspectra at two sets of parameters: one is at $Re=1000$, $\alpha =0.2$, $Wi=60$ and the other is at $Re=260$, $\alpha =0.94$, $Wi=14.5$. For both of them, there is a single unstable mode and the convergence has been examined by increasing the total number of nodes $N$ (we have additionally used the Petrov–Galerkin method (Meseguer & Trefethen Reference Meseguer and Trefethen2003) to double-check the existence and convergence of this mode using the code in Zhang (Reference Zhang2021); the details are not shown). Because this finding is only a byproduct of our numerical investigation of weakly nonlinear stability of viscoelastic pipe flows, we place the discussions of the linear instability of UCM fluids in Appendix E.

Figure 5.

Eigenspectra obtained from a linear stability analysis of pipe Poiseuille flows of UCM fluids: (a) at $Re=1000$, $\alpha =0.2$ and $Wi=60$; (b) at $Re=260$, $\alpha =0.94$ and $Wi=14.5$.

4.3. Bifurcation studies in a weakly nonlinear framework

The present weakly nonlinear stability analysis is performed around the linear critical condition. For all the results presented here, the neutral mode with almost zero growth rate is well separated from the continuous spectrum (by using a sufficiently large $N$) and, thus, its eigenfunction can be safely used to construct the higher-order solutions for the evaluation of the Landau coefficient. In what follows, we will successively discuss supercritical/subcritical bifurcations and a scaling law for the Landau coefficient. We find that the main parameter that differentiates the bifurcation types appears to be $\beta$. Thus, we will first consider different ranges of $\beta$ in these two scenarios and then delimit the boundary between them.

4.3.1. Supercritical bifurcations at small $\beta$

As mentioned in the previous sections, the real part of the Landau coefficient $a_3$ (denoted as $a_{3,r}$) determines the bifurcation type: a negative (positive) value indicates a supercritical (subcritical) bifurcation; see (2.24). In this subsection, we consider relatively small values of $\beta$. Firstly, the critical Reynolds number $Re_c$ and wavenumber $\alpha _c$ at which the Landau coefficient $a_3$ is evaluated are shown in figure 6(a,b). On the one hand, at a fixed $\beta$, increasing the Weissenberg number $Wi$ leads to high $Re_c$ and low $\alpha _c$, implying a strong stabilizing effect of polymer elasticity on the laminar base flow. On the other hand, at a fixed $Wi$, increasing polymer concentrations (decreasing $\beta$ from approximately 0.7) also results in an increase of $Re_c$ and a decrease of $\alpha _c$. Secondly, figure 6(c) shows the corresponding real parts (growth rate) of the eigenvalue $\mu _r$ of the almost neutral modes. For all the five values of $\beta$, the growth rates are all indeed very small of the order of $O(10^{-11})$, around the linear instability onset. Thirdly, figure 6(d) illustrates the corresponding values of $a_{3,r}$ for different values of $Wi$ and $\beta$. Clearly, all the curves fall in the negative range of $a_{3,r}$, suggesting that supercritical bifurcations exist in the parameter space explored here. Thus, the lowest-order nonlinear term in the governing equation stabilizes the disturbances and drives them to saturate. Besides, increasing $Wi$ or decreasing $\beta$ generally corresponds to a smaller amplitude of $a_{3,r}$, indicating that the degree of supercriticality is reduced. Thus, the stabilizing effect of nonlinearity becomes weaker when increasing $Wi$ or decreasing $\beta$ in their respective ranges, as shown in figure 6(d). From the slopes of these curves, one may infer that the degree of supercriticality probably approaches an asymptote at sufficiently high $Wi$ for each $\beta$.

Figure 6.

Critical conditions and the Landau coefficient for pipe flow of Oldroyd-B fluids at various $Wi$ and $\beta$ (small and moderate values): (a) critical Reynolds number $Re_{c}$; (b) critical wavenumber $\alpha _{c}$; (c) linear growth rate $\mu _{r}$; (d) Landau coefficient $a_{3,r}$. All the cases here are supercritical as $a_{3,r}$ is negative.

In the UCM limit $\beta =0$, we have reported the linear instability in the last subsection. From the neutral curves (see figure 19 in Appendix E), the critical condition $Re_{c}$ and $\alpha _{c}$ at the leftmost end of the multiple loops can be determined, as shown in figure 7(a,b). Basically, decreasing $Wi$ leads to smaller $Re_c$ and larger $\alpha _c$. One can observe a discontinuity in the curve around $Wi=19$ whose origin results from the complexity of the neutral curves at low $Wi$, as we explain now. The neutral curve at $Wi=19$ (shown in figure 19b) consists of two separate loops. The left tips of them are denoted as point A (for the smaller loop) and point B (for the larger loop), respectively. In figure 7 the corresponding $Re$ and $\alpha$ at point A and point B are represented by red filled squares and blue filled triangles at $Wi=19$ (see the vertical lines in panels $(a)$ and $(b)$ and their insets). When $Wi>19$, we use the critical point A in the weakly nonlinear analysis because it is the lowest $Re$ which is the critical $Re_c$ in a global sense. With $Wi$ reducing to $18$ and smaller, the smaller loop disappears and we thus use the critical point B in the weakly nonlinear stability analysis, resulting in the apparent jump of $\alpha _{c}$ by about 0.1. Such discontinuity can also be seen from the group velocity $c_{g,r}$ in figure 7(c), where with decreasing $Wi$ from 150 the group velocity $c_{g,r}$ firstly goes down and then suddenly increases at $Wi=19$ and decreases again. As for the Landau coefficient, $a_{3,r}$ undergoes a minimum at about $Wi=25$ and also there is a clear discontinuity at $Wi=19$. In spite of such complex behaviour, all the values of $a_{3,r}$ for $Wi$ ranging from 14.5 to 150 are negative, again suggesting supercritical bifurcations for UCM pipe Poiseuille flow at the inertial regime. In the literature, Meulenbroek et al. (Reference Meulenbroek, Storm, Bertola, Wagner, Bonn and van Saarloos2003) reported subcritical bifurcations in the inertialess regime in their weakly nonlinear analysis of the viscoelastic pipe Poiseuille flow of UCM fluids. The weakly nonlinear analysis was also extended to plane Poiseuille flow (Meulenbroek et al. Reference Meulenbroek, Storm, Morozov and van Saarloos2004) and plane Couette flow (Morozov & van Saarloos Reference Morozov and van Saarloos2005) of UCM fluids, both resulting in subcritical bifurcations in the inertialess regime. It is thus interesting as a future work to conduct weakly nonlinear stability analysis of these flows with relatively strong inertia to assess its effect.

Figure 7.

Critical conditions and the Landau coefficient at various $Wi$ for pipe flow of UCM fluids ($\beta =0$): (a) critical Reynolds number $Re_{c}$; (b) critical wavenumber $\alpha _{c}$; (c) group velocity $c_{g,r}$; (d) Landau coefficient $a_{3,r}$. All the cases here are supercritical as $a_{3,r}$ is negative.

Next, we discuss the possible link between our theoretical prediction of flow bifurcations and the flow transitions observed in experiments. As briefly reviewed in the introduction, Samanta et al. (Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013) observed no hysteresis of pressure fluctuations in pipe Poiseuille flow at high polymer concentrations of 500 p.p.m. solution (approximately equivalent to $\beta =0.69$), suggesting the existence of supercritical bifurcations (see figure 2B in their paper). In this case, the transitional Reynolds number is $Re_t \approx 800$ at which our weakly nonlinear analysis also predicts a supercritical bifurcation but at a high Weissenberg number $Wi \approx 250$ (at which the linear growth rate is almost zero). Note that we confine ourselves to the linear critical condition in our calculation in order that weakly nonlinear results are guaranteed to converge. If we use $Wi=15$ (the maximum value of transition $Wi$ in Samanta et al. Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013) in our weakly nonlinear analysis anyway, we can formally also get a negative $a_{3,r}$ (again the convergence of this result may be questionable). In any case, in experiments, infinitesimal perturbations are difficult to realise and realistic disturbances are always to some extent of finite amplitude, so it is difficult to differentiate the genuine supercriticality from the subcriticality which is very sensitive to noise-level disturbances (see Appendix A for more discussions). But the present theoretical prediction seems to lend some supports to a supercritical bifurcation to a certain degree.

Another experimental study conducted very recently by Choueiri et al. (Reference Choueiri, Lopez, Varshney, Sankar and Hof2021) showed that in the pipe flow of a 600 p.p.m. polymer solution, the measured pressure fluctuation amplitude $p_m$ grows continuously with increasing $Re$ near the instability onset $Re_t\approx 18$, following $p_m \propto \sqrt {Re-Re_t}$ which signifies a supercritical transition. In this experiment, the viscosity ratio $\beta$ is about 0.56, at which our weakly nonlinear analysis predicts supercritical bifurcations in a wide range of $Wi$ but at a higher critical Reynolds number $Re_c \sim 150$ than the instability threshold $Re_t=18$ observed in experiments. Concerning this difference, it should be noted that, in view of the experimental fact that the onset of structural turbulence is characterized by a critical shear rate instead of a critical Reynolds number (Samanta et al. Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013), the critical Reynolds number actually varies with pipe radius as $Re_c \sim R^2$ (Ram & Tamir Reference Ram and Tamir1964). Therefore, the low Reynolds number $Re_t=18$ is expected to become higher (may exceed $Re \sim O(150)$ predicted by linear theory) should a pipe with a larger radius be used. Despite the dependency of $Re_c$ on the flow geometry, our prediction of the intrinsic bifurcation type seems to be consistent with their observations. In fact, when Choueiri et al. (Reference Choueiri, Lopez, Varshney, Sankar and Hof2021) investigated the flow pattern at various $Re$ (at a viscosity ratio of about $\beta =0.18$ at which our calculation also predicts supercritical bifurcations), they observed the scaling $p_m \propto \sqrt {Re-Re_t}$ near the onset, which seems to support a supercritical bifurcation at high polymer concentrations (even though the authors were cautious enough to not exclude possibilities of other scaling relations due to the possibly sensitive nature of the viscoelastic pipe flow to the disturbances). As for the flow pattern, they observed chevron shape streaks at low $Re$, which resembles the flow patterns predicted by linear stability theory. These flow structures correspond to our linear solution in (2.15) (and we have checked that they look similar). We now present the second-order solution in (2.16). The axial velocity $\tilde {u}_{z,21}$ at the same parameter $Wi=203, \beta =0.56$ as in Choueiri et al. (Reference Choueiri, Lopez, Varshney, Sankar and Hof2021) is shown in figure 8(a). As one can see, they also exhibit such kind of chevron shape patterns. The same figure also shows the contours of the trace of the conformation tensor $\textrm {tr}(\tilde {\boldsymbol {c}}_{21})$ figure 8(b). It is clear that these polymers are stretched along the streamwise direction.

Figure 8.

Flow pattern of the second-order solution at $Re_c = 1054.35$, $\alpha _{c}=0.1359$, $Wi=203$ and $\beta =0.56$. Results for (a) $\tilde {u}_{z,21}/U_z$; (b) $\textrm {tr}(\tilde {\boldsymbol {c}}_{21})$ along with the blue curves for streamlines.

4.3.2. Subcritical bifurcations at large $\beta$

In the last subsection the viscosity ratio $\beta$ is restricted to small and moderate values and we found supercritical bifurcations in viscoelastic pipe flows in the weakly nonlinear stability analysis. In this subsection our focus is on relatively larger values of $\beta$, corresponding to low polymer concentrations.

The effects of $\beta$ and $Wi$ on the critical Reynolds number and wavenumber, as shown in figure 9(a,b), are similar to those at low viscosity ratios: the flow is stabilized at higher values of $Wi$ and lower values of $\beta$ (both can be interpreted as strong viscoelastic effects due to polymers). The asymptotic trend shows that the critical $Re_c$ may be very large for sufficiently high $Wi$. As for the Landau coefficient $a_{3,r}$ shown in figure 9(c), it can be seen that its values are all positive, indicating subcritical bifurcations in the viscoelastic pipe flow with low polymer concentrations. The degree of subcriticality (characterized by the amplitude of $a_{3,r}$) becomes smaller for higher $Wi$ or lower $\beta$. At high $\beta$ (or low polymer concentrations), Newtonian-like regimes in the form of hysteresis and intermittency were observed in the experiments conducted by Samanta et al. (Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013). Hysteresis was observed in their experiment of viscoelastic pipe flows at a low polymer concentration of 100 p.p.m. (equivalent to $\beta \approx 0.92$, see their figure 2A). Extrapolation of the Landau coefficient $a_{3,r}$ based on figure 9(c) suggests also subcritical bifurcations at $\beta \approx 0.92$ for a wide range of Weissenberg numbers, qualitatively consistent with the experimental observations. More specifically, at the two ends of the $Wi$-range in figure 9(c), for $\beta =0.92$, the Landau coefficient at $Wi=600, Re_c=325.106379,\alpha _{c}=0.203808$ is $a_{3,r}=202$; the Landau coefficient at $Wi=70, Re_c=87.694576,\alpha _{c}=1.863525$ is $a_{3,r}=1655$, both showing subcritical bifurcations. Therefore, the current calculation supports the subcritical bifurcation in viscoelastic pipe flows at relatively high $\beta$.

Figure 9.

Critical conditions and the Landau coefficient for pipe flow of Oldroyd-B fluids at various $Wi$ and $\beta$ (large values): (a) critical Reynolds number $Re_{c}$; (b) critical wavenumber $\alpha _{c}$; (c) Landau coefficient $a_{3,r}$. All the cases here are subcritical as $a_{3,r}$ is positive.

Since both supercritical and subcritical bifurcations exist in the present viscoelastic pipe flows, there must be a boundary in the parameter space delimiting the bifurcation types, and this will be discussed in the following subsection.

4.3.3. Bifurcation boundary delimited by $\beta$

Figures 10(a) and 10(b) illustrate how the critical parameters $Re_c$ and $\alpha _c$ vary with the viscosity ratio $\beta$. This variation has been observed and discussed in the previous subsections and, thus, will not be repeated here. What is of significance here is the continuous but non-monotonic change of $a_{3,r}$ with the variation of $\beta$ from 0.1 to 0.9 as plotted in figure 10(c). Two key observations can be made as follows. Firstly, at a fixed $Wi$, with $\beta$ increased, the Landau coefficient $a_{3,r}$ decreases first and reaches a minimum and then increases; at about $\beta \approx 0.785$, the sign of $a_{3,r}$ changes from negative to positive, suggesting the change of bifurcation type from supercritical to subcritical. Figure 10(d) shows a close-up of the critical $\beta$ at which the bifurcation type changes within $0.75 \leqslant \beta \leqslant 0.80$. It can be seen that for $300\leqslant Wi\leqslant 650$, the bifurcation boundary is almost fixed at $\beta \approx 0.785$ while it slightly reduces for $Wi=150$. In general, the alteration of bifurcation type from supercritical to subcritical can be attributed to the interplay between the flow field and polymer dynamics (as $\beta$ is related to the ratio between solvent viscosity and polymer viscosity). At large $\beta$, viscous effects mainly result from the solvent and thus cause Newtonian-like subcriticality, while polymer dynamics becomes more significant at small $\beta$, resulting in the change to the supercritical bifurcation. Secondly, by comparing the curves at different $Wi$, one can find that the degrees of both supercriticality and subcriticality reduce at high Weissenberg numbers (e.g. the curve for $Wi=650$ is closer to the horizontal reference line $a_{3,r}=0$ than that for $Wi=150$). This indicates that the destabilizing and stabilizing of the nonlinear effects in a viscoelastic pipe flow are more prone to the changes of polymer concentration ($1-\beta$) when the elastic effect is relatively small, whereas their degrees are less influenced by the polymer concentration when the elastic effect is strong. The same data of $Re_c$ and $\alpha _c$ are further plotted in $Wi\text {--}\beta$ planes, as shown in figure 11(a–d). The panels (b,d) extend the range of $Wi$ to lower values where we cannot find linear instability for some $\beta \in [0.1,0.9]$. The meaning of the red dashed lines is the bifurcation boundary to be further discussed in figure 12 shortly. In panels (e,f) we present the data in terms of $E_c=Wi/Re_c$ (not to be confused with a general elastic number $E=Wi/Re$ and note that in this work we fix $Wi$ to perturb $Re$, instead of fixing $E$ to perturb $Re$). We can see that in the range of high $Wi$ (see figure 11e), these contours are almost vertical, showing that $E_c$ is basically independent of $Wi$ at the linear critical conditions. Moreover, the bifurcation boundary $\beta \approx 0.785$ almost coincides with the contour level of $E_c \approx 0.46$; high levels of $E_c$ are always related to high values of $\beta$.

Figure 10.

Continuous variation of the Landau coefficient with viscosity ratio $\beta$ at various Weissenberg numbers $Wi$: (a) critical Reynolds number $Re_{c}$; (b) critical wavenumber $\alpha _{c}$; (c) Landau coefficient $a_{3,r}$; (d) Landau coefficient $a_{3,r}$ near the bifurcation boundary. Data points below the horizontal dashed line in (c,d) denote supercritical bifurcations, while those above the line denote subcritical bifurcations.

Figure 11.

Critical conditions in a large $(Wi,\beta )$ plane. (a,b) Critical Reynolds number $Re_c$; (c,d) critical wavenumber $\alpha _c$; (e,f) elasticity number at the critical condition $E_c=Wi/Re_c$. The right column is presented because for some values of $Wi,\beta$, there is no linear instability. The dashed red line labels the bifurcation boundary, see the discussion in figure 12.

In figure 12(a) the values of $a_{3,r}$ are displayed and the bifurcation boundary within the parameter space of $100\leqslant Wi\leqslant 650$ and $0.1 \leqslant \beta \leqslant 0.9$ is illustrated in a $(Wi,\beta )$ plane. The dashed red line delimits the bifurcation types in this plane; it is almost a vertical line, with $\beta \approx 0.785$. The highest degree of supercriticality (white area) is located at the smallest $Wi=100$ and moderate $0.3 \lesssim \beta \lesssim 0.7$, while the highest degree of subcriticality (black area) can be found at the largest $\beta$ and relatively lower $Wi$. When $Wi$ is less than about 60, linear instability ceases to exist for some values of $\beta$ and, thus, we do not apply weakly nonlinear analysis to these cases. These results are presented in curves in panel ($b$). Within the parameter range of $20\leqslant Wi\leqslant 60$ and $0.1\leqslant \beta \leqslant 0.9$, data points not displayed correspond to the absence of linear instability. For example, figure 11(b) shows that with the decreasing of $Wi$ down to 20, the linear instability only exists at low $\beta$, implying that polymer effects are essential to the linear instability. This scenario agrees with that in Chaudhary et al. (Reference Chaudhary, Garg, Subramanian and Shankar2021) where a lower $\beta$ is needed to render the flow linearly unstable at smaller $Wi$ (see their figure 23b). As for the Landau coefficient $a_{3,r}$ shown in figure 12(b), its variation with $\beta$ at various $Wi$ is similar to that observed in figure 10(c) with a minimum at approximately $\beta \in (0.5,0.6)$ and a change of bifurcation type at approximately $\beta \approx 0.785$.

Figure 12.

Landau coefficient in a $(Wi,\beta )$ plane at linear critical conditions for high $Wi$ (a) and in a curve plot for low $Wi$ (b). The dashed red line representing $\beta \approx 0.785$ in (a) marks the bifurcation boundary, the same as those in figure 11; this boundary almost coincides with the contour of $E_c \approx 0.46$ in figure 11(e).

Based on the above results, one can see that the key parameter governing the bifurcation type is the viscosity ratio $\beta$ (sweeping in $Wi$ does not change the bifurcation type in general). The change of the bifurcation type due to the viscosity ratio $\beta$ seems to be supported by the experimental observations in Samanta et al. (Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013), where by increasing the polymer concentration (equivalent to decreasing $\beta$) from the Newtonian fluid, hysteresis behaviour can be firstly observed but then disappears when the polymer concentration is beyond 200 p.p.m. ($\beta \approx 0.855$). Although their critical viscosity ratio $\beta \approx 0.855$ (for the existence of hysteresis or not) is slightly larger than our prediction $\beta \approx 0.785$, both their experiments and our predictions imply the existence of a bifurcation boundary defined by $\beta$. The reason for the discrepancy may be that they did not dedicate themselves in determining the critical value of $\beta$ for the bifurcation type and that our results are based on Oldroyd-B axisymmetric flows (other more realistic models such as FENE-P may improve the results). A similar scenario is observed in a more recent experiment by Choueiri, Lopez & Hof (Reference Choueiri, Lopez and Hof2018) where Newtonian-like transition (known to be subcritical) persists at low polymer concentrations but is replaced by a qualitatively different (likely to be supercritical) instability at high polymer concentrations. In addition, Chandra et al. (Reference Chandra, Shankar and Das2020) has also suggested the possible crossover of transition type from subcritical to supercritical with the former at low polymer concentrations and the latter at high polymer concentrations in viscoelastic microtube flows. For viscoelastic channel flows, Graham (Reference Graham2014) depicted a phase diagram demonstrating that: subcritical bifurcations exist at small elasticity numbers $E$ where the turbulence is Newtonian-like; meanwhile, increasing $E$ to a sufficiently high value results in early turbulence. He speculated that there is a possible linear instability, which has been confirmed theoretically by Garg et al. (Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018) and Chaudhary et al. (Reference Chaudhary, Garg, Subramanian and Shankar2021). A similar diagram for viscoelastic pipe flow in Chaudhary et al. (Reference Chaudhary, Garg, Subramanian and Shankar2021) showed that subcritical axisymmetric elasto-inertial structures may exist at low elasticity numbers $E$, while they may change to be supercritical at high $E$. In this diagram, the bifurcation boundary seems to be defined by $E$, instead of $\beta$. But it should be noted that their diagram is mainly for $\beta \sim 0.60$ and the situation is likely to be different for other values of $\beta$, according to our results here.

As elucidated above, a single parameter, viscosity ratio $\beta$ in the present case, can change the bifurcation type of viscoelastic pipe flows. Simple flows, such as Newtonian plane Poiseuille flow and pipe Poiseuille flow, are governed by a single control parameter $Re$ and there is only one type of bifurcation. As the flow system becomes complex, where, e.g. two and more competing mechanisms are present, the bifurcation type may be different in different parts of the parameter space. For example, a weakly nonlinear analysis by Graham (Reference Graham1998) showed that, in circular Couette flows (with imposed axial Poiseuille flow) of UCM fluids, increasing the axial Weissenberg number can change the bifurcation type from subcritical to supercritical. Such change of bifurcation types also exists in Dean flows with an imposed axial Poiseuille flow of both UCM and Oldroyd-B fluids (Ramanan, Kumar & Graham Reference Ramanan, Kumar and Graham1999). Fujimura & Kelly (Reference Fujimura and Kelly1997) found that for stably stratified plane Poiseuille flows, the subcritical bifurcation exists in a wide range of Prandtl numbers while the bifurcation changes to be supercritical in a narrow range of $Pr<0.17$. In this case, the bifurcation boundary in the parameter space is defined by both the Prandtl number and the Richardson number (see their figure 5). For another example of a numerical study of electro-convection of viscoelastic fluids between two horizontal plates subjected to unipolar injection, Su et al. (Reference Su, Zhang, Luo and Yi2020) found that, at $\beta =0.8$, the intrinsic subcritical bifurcation at $Wi=0.2$ changes to be supercritical when $Wi$ increases to 0.5 (see their figures 5 and 6). In this flow, it is $Wi$ that alters the bifurcation type, instead of $\beta$. In addition, it is worth mentioning that the supercritical bifurcations in Su et al. (Reference Su, Zhang, Luo and Yi2020) are also caused by viscoelastic effects, consistent with our present finding about supercriticality at low and moderate $\beta$ where viscoelastic effects are strong. Another example is taken from the work by Wu et al. (Reference Wu, Traoré, Pérez and Zhang2015), where a multi-physics problem, i.e.  the electro-thermo-convection between two infinite plane plates subjected to unipolar charge injection, was investigated. They found that the bifurcation type can be changed by varying the electric Rayleigh number or Prandtl number or ion mobility number (see the definitions of these dimensionless parameters in their paper). Such complex characteristics of the flow result from the superposition of two different systems: one is the Rayleigh–Bénard convection, well known to be supercritical; the other is the unipolar charge injection induced convection, which is subcritical (Zhang Reference Zhang2016). From this point of view, our present results of the change of bifurcation type can also be interpreted as an intermediate state between the subcritical Newtonian pipe flow and the supercritical UCM pipe flow. Thus, at finite viscosity ratios, there exists a change of the bifurcation type.

4.3.4. A scaling law of the Landau coefficient $a_3$

In the end, we report a scaling law of the Landau coefficient $a_{3}$. In the results of linear stability analysis, Chaudhary et al. (Reference Chaudhary, Garg, Subramanian and Shankar2021) has reported scalings of critical Reynolds number $Re_c \propto [E(1-\beta )]^{-3/2}$ and critical wavenumber $\alpha _c \propto [E(1-\beta )]^{-1/2}$ for $E(1-\beta ) \ll 1$; they justified these scalings by means of a scaling analysis for the boundary layer near the pipe centreline $r=0$. For the present weakly nonlinear analysis, we have also observed a scaling of the Landau coefficient $a_3$, as shown in figure 13(a). The real part $a_{3,r}$ scales with $1/Wi$ (with $Wi\in [100,650]$) at fixed values of $\beta$ and the scaling also holds for the imaginary part $a_{3,i}$ (not shown). At a first glance, the scaling at a large value of $\beta$ looks inferior; this is because the $Wi$ is not large enough for these $\beta$. Figure 13(b) shows the variation of $a_{3,r}Wi$ as a function of $Wi$ at several viscosity ratios. When $\beta =0.1$, $0.3$, $0.5$ or $0.7$, the variation of $a_{3,r}Wi$ is large at low Weissenberg numbers $Wi<100$ and the value of $a_{3,r}Wi$ remains almost constant in the range of $100< Wi<650$. Within the same range of $Wi$ for $\beta =0.9$, however, the value of $a_{3,r}Wi$ changes with $Wi$ significantly and only asymptotically approaches a constant value for much higher $Wi$. These results suggest that the scaling result of $a_{3,r}Wi$ is a high-$Wi$ phenomenon (strong elasticity) and larger $\beta$ (lower polymer concentrations) entails higher $Wi$ to present the scaling law. Considering the scaling results in the linear regime found by Garg et al. (Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018) and Chaudhary et al. (Reference Chaudhary, Garg, Subramanian and Shankar2021) and the scaling in the weakly nonlinear regime discovered in the current work, it would be very interesting to reveal the connection of these two scalings; however, due to space limit, this direction of work, and more possible scalings, will be pursued in a future work.

Figure 13.

A scaling law of the Landau coefficient $a_3$: (a) collapse of the $a_{3,r}Wi$-$\beta$ curves for different $Wi$; (b) values of $a_{3,r}Wi$ as a function of $Wi$ for some $\beta$. The horizontal dashed line in (a) marks the bifurcation boundary.

5. Conclusions

In this paper we have performed a weakly nonlinear stability analysis of viscoelastic pipe flows of Oldroyd-B fluids based on a multiple-scale expansion method. Specifically, the disturbance is restricted to be axisymmetric, the same as that in Garg et al. (Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018) where the linear instability was reported for the first time. It is the existence of such linear instability that makes our weakly nonlinear analysis feasible and, thus, the convergence of multiple-scale expansion near the linear critical point can be guaranteed. This is a different situation compared with Meulenbroek et al. (Reference Meulenbroek, Storm, Bertola, Wagner, Bonn and van Saarloos2003,  Reference Meulenbroek, Storm, Morozov and van Saarloos2004) where neutral curves do not exist and they had to consider bifurcation from infinity. Resulting from the multiple-scale expansion of space, time, disturbance and control parameter, a complex GLE governing the evolution of disturbance amplitude has been derived. The corresponding Landau coefficient has been evaluated based on two different formulations of the same problem, yielding results in good agreement. In addition, DNS have been conducted to verify the bifurcation type predicted by the present weakly nonlinear analysis. The key finding is that both supercritical and subcritical bifurcations exist in viscoelastic pipe flows and the viscosity ratio (or the polymer concentration) is the control parameter defining the bifurcation type.

The present weakly nonlinear analysis has been performed around the linear instability point. In our revisit of the linear eigenvalue problem, the results we have obtained (in terms of the eigenspectra and neutral curves) agree quantitatively well with those in the literature (Garg et al. Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018; Chaudhary et al. Reference Chaudhary, Garg, Subramanian and Shankar2021). Moreover, we have found that the pipe Poiseuille flow of UCM fluids (with $\beta =0$) can be linearly unstable and such linear instability ceases to exist for $Wi<14.5$. The corresponding eigenfunction is not localized near the pipe axis, instead, it goes through the whole pipe cross-section. This result supplements those in Garg et al. (Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018) and Chaudhary et al. (Reference Chaudhary, Garg, Subramanian and Shankar2021) where the authors find stable modes in a wide range of governing parameters for the UCM fluids. In addition, the corresponding weakly nonlinear results show that the bifurcation in UCM pipe Poiseuille flow in the inertial regime is supercritical. It is noted that Meulenbroek et al. (Reference Meulenbroek, Storm, Bertola, Wagner, Bonn and van Saarloos2003) reported subcritical instability of pipe Poiseuille flow of UCM fluids in the inertialess limit; and its counterpart in channel is also subcritical (Meulenbroek et al. Reference Meulenbroek, Storm, Morozov and van Saarloos2004).

With increasing $\beta$ from the UCM limit, the above-mentioned supercritical bifurcation continues to exist at high polymer concentrations (small $\beta$) but the flow becomes subcritical for low polymer concentrations. In the experimental study of viscoelastic pipe Poiseuille flow by Samanta et al. (Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013), they also observed hysteresis in changing $Re$ at low polymer concentrations (signifying subcriticality), and the hysteresis disappeared at high polymer concentrations (implying supercriticality despite experimental uncertainties). These experimental results indicate that changing polymer concentration (or $\beta$) may affect the role the (lowest-order) nonlinearity plays in the flow bifurcation (being stabilizing as in supercritical bifurcations or destabilizing as in subcritical bifurcations). In terms of the bifurcation boundary defined by the viscosity ratio $\beta$, the present analysis predicts $\beta \approx 0.785$ (above which, the flow is subcritical and below which the flow is supercritical), close to the experimental value of $\beta \approx 0.855$ in Samanta et al. (Reference Samanta, Dubief, Holzner, Schäfer, Morozov, Wagner and Hof2013). The discrepancy may be attributed to several reasons, including that the experiments were not designed particularly for determining this bifurcation boundary in $\beta$ and that we have only considered the Oldroyd-B model. Another experimental study by Chandra et al. (Reference Chandra, Shankar and Das2020) on the flow transition in viscoelastic mircotube flows is also supportive of such a bifurcation boundary in terms of $\beta$. Our results can be viewed as a theoretical exploration of the supercritical and subcritical bifurcations in viscoelastic pipe flows from the perspective of the equations. Lastly, a scaling law of $a_3$ with $Wi^{-1}$ is observed for values of $\beta \in [0.1,0.9]$ when $Wi$ is sufficiently large. It is conjectured that this scaling in the weakly nonlinear stability analysis is related to the scaling laws in the linear regime as found by Garg et al. (Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018) and Chaudhary et al. (Reference Chaudhary, Garg, Subramanian and Shankar2021). This newly discovered scaling law in the weakly nonlinear phase will be explained in a future work.

As mentioned above, the current results are based on viscoelastic pipe flows of Oldroyd-B fluids. For more accurate quantitative comparisons to experiments, more realistic constitutive models deserve to be further analysed. For example, in Chaudhary et al. (Reference Chaudhary, Garg, Subramanian and Shankar2021) the authors had to resort to the FENE-P model to consider the shear-thinning effect in order to better match the scaling in experiments. The more realistic constitutive model may improve the prediction of the bifurcation boundary we found in this work. Besides, it may also be interesting to investigate the bifurcation boundary in other viscoelastic flows and possible scaling laws of the Landau coefficients there. In revising this paper, we found that Buza, Page & Kerswell (Reference Buza, Page and Kerswell2021) submitted a weakly nonlinear stability analysis of viscoelastic channel flows. Lastly, one can also consider to perturb $Re$ or $Wi$ with fixed $E$ and $\beta$, as the elastic number can facilitate the interpretation of the results.

Funding

D.W. is supported by a PhD scholarship (No. 201906220200) from the China Scholarship Council and the NUS research scholarship. M.Z. acknowledges the financial support of Ministry of Education, Singapore (a Tier 2 grant with the WBS no. R-265-000-661-112).

Declaration of interests

The authors report no conflict of interest.

Appendix A. An illustration of bifurcation types

We present here the bifurcation types as shown in figure 14. The solid curves represent stable solutions and the dashed curves unstable solutions. The grey dotted circle represents a neighbourhood of the linear critical condition, where our computation is carried out. Panel ($a$) shows a supercritical bifurcation (with $a_{3,r}<0$ in (2.24)) and panel ($b$) a subcritical one (with $a_{3,r} > 0$). The arrows in the figures indicate that unstable solutions will eventually become stable (linearly or nonlinearly), but via different bifurcation routes. Routes 1 to 4 represent stable spiral, supercritical bifurcation, subcritical bifurcation and stable spiral, respectively. Panel ($c$) shows a subcritical bifurcation but may appear as a supercritical one in experiments (black dots). The red dotted line represents the disturbance level of the system. This case is relevant to a flow which is very sensitive to the external disturbance (including the systematic one).

Figure 14.

Illustration of different bifurcation types ($y$-axis: disturbance amplitude $|A|$; $x$-axis: the control parameter $Re$) near the linear critical points (marked with grey dotted circles). Solid curves are stable solutions and dashed curves unstable solutions; arrows indicate the evolution direction of solutions. (a) Supercritical bifurcation ($a_{3,r}<0$); (b) subcritical bifurcation ($a_{3,r}>0$); (c) a subcritical bifurcation but may appear to be supercritical in experiments (black dots). The red dotted line indicates the systematic disturbance, which cannot be entirely removed. The black dots represent experiment measurements with arbitrary error bars.

Appendix B. Governing equations in the $\psi$-$g$ formulation

The short notations $F(g_{rr})$ and $F(g_{rz})$ in the governing equation of $\psi$ are

(B1a)\begin{align} F(g_{rr})&=C_{rz} \frac{\partial^2 g_{rr}}{\partial r^2} + \left( 2 C'_{rz} + \frac{C_{rz}}{r} \right) \frac{\partial g_{rr}}{\partial r} + \left( C''_{rz} + \frac{C'_{rz}}{r} - \frac{C_{rz}}{r^2} \right) g_{rr} \nonumber\\ &\quad + \left(\frac{C_{rz}^2}{C_{rr}} - C_{rr}\right) \frac{\partial^2 g_{rr}}{\partial r \partial z} + \left( \left(\frac{C_{rz}^2}{C_{rr}}\right)' - C'_{rr} - \frac{C_{rr}}{r}\right) \frac{\partial g_{rr}}{\partial z} - C_{rz}\frac{\partial^2 g_{rr}} {\partial z^2}, \end{align}
(B1b)\begin{align} F(g_{rz})&=S \frac{\partial^2 g_{rz}}{\partial r^2} + \left( 2 S' + \frac{S}{r} \right) \frac{\partial g_{rz}}{\partial r} + \left( S'' + \frac{S'}{r} - \frac{S}{r^2} \right) g_{rz} \nonumber\\ &\quad + \frac{2S C_{rz}}{C_{rr}} \frac{\partial^2 g_{rz}}{\partial r \partial z} + \left(\frac{2S C_{rz}}{C_{rr}}\right)' \frac{\partial g_{rz}}{\partial z} - S\frac{\partial^2 g_{rz}} {\partial z^2}, \end{align}

and the governing equations of $\boldsymbol {g}$ using the geometric decomposition are

(B2a)\begin{equation} \frac{\partial g_{rr}}{\partial t} ={-} U_z \frac{\partial g_{rr}}{\partial z} - \frac{ g_{rr}}{Wi} - 2 \frac{\partial^2 \psi}{\partial r \partial z} + \frac{C'_{rr}}{C_{rr}} \frac{\partial \psi}{\partial z} - \frac{2 C_{rz}}{C_{rr}} \frac{\partial^2 \psi}{\partial z^2} - \frac{N_{c_{rr}}}{C_{rr}}, \end{equation}
(B2b)\begin{align} \frac{\partial g_{rz}}{\partial t} &={-} U_z \frac{\partial g_{rz}}{\partial z} - \frac{ g_{rz}}{Wi} + \frac{C_{rr}}{S} \frac{\partial^2 \psi}{\partial r^2} + \frac{C_{rr}}{ r S} \frac{\partial \psi}{\partial r} + \frac{2C_{rz}}{S} \frac{\partial^2 \psi}{\partial r \partial z} -\frac{C_{rr}}{r^2 S } \psi \nonumber\\ &\quad + \left( \frac{C_{rz}}{ r S} + \frac{C'_{rz}}{S} - \frac{C'_{rr} C_{rz}}{C_{rr} S} \right) \frac{\partial \psi}{\partial z} + \left( \frac{2 C_{rz}^2}{C_{rr} S} - \frac{ C_{zz}}{S} \right) \frac{\partial^2 \psi}{\partial z^2}\nonumber\\ &\quad + \frac{C_{rr} U'_z}{S} g_{rr} + \frac{C_{rz} N_{c_{rr}}}{C_{rr}S} - \frac{N_{c_{rz}}}{S}, \end{align}
(B2c)\begin{equation} \frac{\partial g_{\theta\theta}}{\partial t} ={-} U_z \frac{\partial g_{\theta\theta}}{\partial z} - \frac{ g_{\theta\theta}}{Wi} + \left(\frac{C'_{\theta\theta}}{C_{\theta\theta}} - \frac{2}{r}\right) \frac{\partial \psi}{\partial z} - \frac{N_{c_{\theta\theta}}}{C_{\theta\theta}}, \end{equation}
(B2d)\begin{align} \frac{\partial g_{zz}}{\partial t} &={-} U_z \frac{\partial g_{zz}}{\partial z} - \frac{ g_{zz}}{Wi} + 2 \frac{\partial^2 \psi}{\partial r \partial z} + \left( \frac{2}{r} + \frac{2 S'}{S} - \frac{C'_{rr}}{C_{rr}} \right) \frac{\partial \psi}{\partial z}\nonumber\\ &\quad + \frac{2 C_{rz}}{C_{rr}} \frac{\partial^2 \psi}{\partial z^2} + \frac{2C_{rr}U'_{z}}{S}g_{rz} - \frac{C_{rz}^2 N_{c_{rr}}}{C_{rr} S^2} + \frac{2C_{rz} N_{c_{rz}}}{ S^2} - \frac{C_{rr} N_{c_{zz}}}{ S^2}. \end{align}

Appendix C. Operators in the $up$-$c$ formulation

C.1. Original operators

In the $up$-$c$ formulation the weight matrix $\boldsymbol {M}$, linear operator $\boldsymbol {L}$ and nonlinear operator $\boldsymbol {N}$ are given as

(C1)\begin{gather} \boldsymbol{M} = \textrm{diag}(I,I,0,I,I,I,I), \end{gather}
(C2)\begin{gather}\boldsymbol{L} = \begin{pmatrix} L_{r,r} & 0 & L_{r,p} & L_{r,rr} & L_{r,rz} & L_{r,\theta\theta} & 0 \\ L_{z,r} & L_{z,z} & L_{z,p} & 0 & L_{z,rz} & 0 & L_{z,zz} \\ L_{p,r} & L_{p,z} & 0 & 0 & 0 & 0 & 0 \\ L_{rr,r} & 0 & 0 & L_{rr,rr} & 0 & 0 & 0\\ L_{rz,r} & L_{rz,z} & 0 & L_{rz,rr} & L_{rz,rz} & 0 & 0\\ L_{\theta\theta,r} & 0 & 0 & 0 & 0 & L_{\theta\theta,\theta\theta} & 0 \\ L_{zz,r} & L_{zz,z} & 0 & 0 & L_{zz,rz} & 0 & L_{zz,zz} \\ \end{pmatrix}, \end{gather}

where the non-zero elements are

(C3)\begin{gather} \left.\begin{gathered} L_{r,r} ={-}U_{z} \frac{\partial }{\partial z} + \frac{\beta}{Re} \left( \nabla^{2} - \frac{1}{r^{2}} \right), \quad L_{r,p}={-} \frac{\partial}{\partial r}, \\ L_{r,rr}=\frac{1-\beta}{ReWi} \left( \frac{\partial }{\partial r} + \frac{1}{r} \right),\quad L_{r,rz}=\frac{1-\beta}{ReWi} \frac{\partial }{\partial z}, \\ L_{r,\theta\theta}=\frac{1-\beta}{ReWi} \left(- \frac{1}{r} \right),\quad L_{z,r}={-}U_{z}', \\ L_{z,z}={-} U_{z} \frac{\partial }{\partial z} + \frac{\beta}{Re} \nabla^{2},\quad L_{z,p}={-} \frac{\partial}{\partial z}, \\ L_{z,rz}=\frac{1-\beta}{ReWi} \left( \frac{\partial }{\partial r} + \frac{1}{r} \right),\quad L_{z,zz}=\frac{1-\beta}{ReWi} \frac{\partial }{\partial z}, \\ L_{p,r}=\frac{\partial }{\partial r} + \frac{1}{r},\quad L_{p,z}=\frac{\partial }{\partial z}, \\ L_{rr,r}={-} C_{rr}' + 2 C_{rr}\frac{\partial }{\partial r} + 2 C_{rz}\frac{\partial }{\partial z},\quad L_{rr,rr}={-} U_{z} \frac{\partial }{\partial z} - \frac{1}{Wi}, \\ L_{rz,r}={-} C_{rz}' + C_{rz}\frac{\partial }{\partial r} + C_{zz}\frac{\partial }{\partial z},\quad L_{rz,z}= C_{rr}\frac{\partial }{\partial r} + C_{rz}\frac{\partial }{\partial z}, \\ L_{rz,rr}= U_{z}',\quad L_{rz,rz}={-} U_{z} \frac{\partial }{\partial z} - \frac{1}{Wi}, \\ L_{\theta\theta,r}={-} C_{\theta \theta}' - 2 C_{\theta \theta}\frac{1}{r},\quad L_{\theta\theta,\theta\theta}={-} U_{z} \frac{\partial }{\partial z} - \frac{1}{Wi}, \\ L_{zz,r}={-} C_{zz}',\quad L_{zz,z}= 2 C_{rz}\frac{\partial }{\partial r}+ 2 C_{zz}\frac{\partial }{\partial z}, \\ L_{zz,rz}= 2 U_{z}',\quad L_{zz,zz}={-} U_{z} \frac{\partial}{\partial z} - \frac{1}{Wi}, \end{gathered}\right\} \end{gather}
(C4)\begin{gather} \boldsymbol{N} ={-} (N_{u_{r}},N_{u_{z}},0,N_{c_{rr}},N_{c_{rz}},N_{c_{\theta \theta}},N_{c_{zz}})^\textrm{T}. \end{gather}

C.2. Multiple-scale expansion of operators

In the $up$-$c$ formulation the weight matrix $\boldsymbol {M}$ is a constant matrix and, thus,

(C5)\begin{equation} \boldsymbol{M}_0 = \boldsymbol{M}, \quad \boldsymbol{M}_1 = \boldsymbol{M}_2 = \boldsymbol{0}, \quad \text{so} \ \boldsymbol{M}_1^\circ = \boldsymbol{M}_2^\circ = \boldsymbol{0}. \end{equation}

The subscale operators for $\boldsymbol {L} = \boldsymbol {L}_{0} + \epsilon \boldsymbol {L}_{1} + \epsilon ^{2} \boldsymbol {L}_{2} + O(\epsilon ^{3})$ are described as follows. Operator $\boldsymbol {L}_{0}$ can be obtained by simply replacing ${\partial }/{\partial z}$ with ${\partial }/{\partial z_0}$, and also replacing $Re$ with $Re_c$ in $\boldsymbol {L}$. Operator $\boldsymbol {L}_{1}$ is not so straightforward and is given as

(C6)\begin{equation} \boldsymbol{L}_1 = \begin{pmatrix} L_{r,r1} & 0 & 0 & 0 & L_{r,rz1} & 0 & 0 \\ 0 & L_{z,z1} & L_{z,p1} & 0 & 0 & 0 & L_{z,zz1} \\ 0 & L_{p,z1} & 0 & 0 & 0 & 0 & 0 \\ L_{rr,r1} & 0 & 0 & L_{rr,rr1} & 0 & 0 & 0\\ L_{rz,r1} & L_{rz,z1} & 0 & 0 & L_{rz,rz1} & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & L_{\theta\theta,\theta\theta 1} & 0 \\ 0 & L_{zz,z1} & 0 & 0 & 0 & 0 & L_{zz,zz1} \\ \end{pmatrix} = \boldsymbol{L}_1^\circ \frac{\partial}{\partial z_1}, \end{equation}

where the non-zero elements are

(C7)\begin{equation} \left.\begin{gathered} L_{r,r1} ={-}U_{z} \frac{\partial }{\partial z_1} + \frac{\beta}{Re_c} \left( 2 \frac{\partial}{\partial z_0} \frac{\partial}{\partial z_1} \right),\quad L_{r,rz1}=\frac{1-\beta}{Re_c Wi} \frac{\partial }{\partial z_1}, \\ L_{z,z1}={-} U_{z} \frac{\partial }{\partial z_1} + \frac{\beta}{Re_c} \left(2 \frac{\partial}{\partial z_0} \frac{\partial}{\partial z_1} \right),\quad L_{z,p1}={-} \frac{\partial}{\partial z_1},\quad L_{z,zz1}=\frac{1-\beta}{Re_c Wi} \frac{\partial }{\partial z_1}, \\ L_{p,z1}=\frac{\partial }{\partial z_1},\quad L_{rr,r1}=2 C_{rz}\frac{\partial }{\partial z_1},\quad L_{rr,rr1}={-} U_{z} \frac{\partial }{\partial z_1}, \\ L_{rz,r1}=C_{zz}\frac{\partial }{\partial z_1},\quad L_{rz,z1}= C_{rz}\frac{\partial }{\partial z_1},\quad L_{rz,rz}={-} U_{z} \frac{\partial }{\partial z_1}, \\ L_{\theta\theta,\theta\theta}={-} U_{z} \frac{\partial }{\partial z_1},\quad L_{zz,z1}= 2 C_{zz}\frac{\partial }{\partial z_1},\quad L_{zz,zz1}={-} U_{z} \frac{\partial}{\partial z_1}. \end{gathered}\right\} \end{equation}

Operator $\boldsymbol {L}_1^\circ$ is constructed by moving ${\partial }/{\partial z_1}$ out of $\boldsymbol {L}_1$. Operator $\boldsymbol {L}_{2}$ is given as

(C8)\begin{equation} \boldsymbol{L}_2 = \begin{pmatrix} L_{r,r2} & 0 & 0 & L_{r,rr2} & L_{r,rz2} & L_{r,\theta\theta 2} & 0 \\ 0 & L_{z,z2} & 0 & 0 & L_{z,rz2} & 0 & L_{z,zz2} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \end{pmatrix} = \boldsymbol{L}_{2Re} + \boldsymbol{L}_{2}^\circ \frac{\partial^2 }{\partial z_1^2}, \end{equation}

where the non-zero elements are

(C9)\begin{align} \left.\begin{gathered} L_{r,r2} ={-}\frac{\beta}{Re_c^2} \left( \frac{\partial^{2}}{\partial r^{2}} + \frac{1}{r}\frac{\partial}{\partial r} + \frac{\partial^{2}}{\partial z_0^{2}} - \frac{1}{r^2}\right) + \frac{\beta}{Re_c} \frac{\partial^{2}}{\partial z_1^{2}}, \\ L_{r,rr2} ={-}\frac{1-\beta}{Re_c^2 Wi} \left( \frac{\partial}{\partial r} + \frac{1}{r} \right), \quad L_{r,rz2} ={-}\frac{1-\beta}{Re_c^2 Wi} \frac{\partial}{\partial z_0}, \quad L_{r,\theta\theta 2} = \frac{1-\beta}{Re_c^2 Wi} \left( \frac{1}{r} \right), \\ L_{z,z2} ={-}\frac{\beta}{Re_c^2} \left( \frac{\partial^{2}}{\partial r^{2}} + \frac{1}{r}\frac{\partial}{\partial r} + \frac{\partial^{2}}{\partial z_0^{2}} \right) + \frac{\beta}{Re_c} \frac{\partial^{2}}{\partial z_1^{2}}, \\ L_{r,rz2} ={-}\frac{1-\beta}{Re_c^2 Wi} \left( \frac{\partial}{\partial r} + \frac{1}{r} \right), \quad L_{r,zz2} ={-}\frac{1-\beta}{Re_c^2 Wi} \frac{\partial}{\partial z_0}, \end{gathered}\right\} \end{align}

and $\boldsymbol {L}_{2}^\circ$ is constructed with the terms having ${\partial ^2 }/{\partial z_1^2}$; and $\boldsymbol {L}_{2Re}$ consists of the remaining terms. In the nonlinear term $\boldsymbol {N}= \epsilon ^{2} \boldsymbol {N}_{2} + \epsilon ^{3} \boldsymbol {N}_{3} + O(\epsilon ^{4})$, the subscales are

(C10a)\begin{gather} \boldsymbol{N}_{2} ={-} (N_{u_{r}2},N_{u_{z}2},0,N_{c_{rr}2},N_{c_{rz}2},N_{c_{\theta \theta}2},N_{c_{zz}2})^\textrm{T}, \end{gather}
(C10b)\begin{gather}\boldsymbol{N}_{3} ={-} (N_{u_{r}3},N_{u_{z}3},0,N_{c_{rr}3},N_{c_{rz}3},N_{c_{\theta \theta}3},N_{c_{zz}3})^\textrm{T}, \end{gather}

where the non-zero elements are given as

(C11)\begin{equation} \left.\begin{gathered} \displaystyle N_{u_{r}2} = u_{r,1} \frac{\partial u_{r,1}}{\partial r} + u_{z,1} \frac{\partial u_{r,1}}{\partial z_0}, N_{u_{z}2} = u_{r,1} \frac{\partial u_{z,1}}{\partial r} + u_{z,1} \frac{\partial u_{z,1}}{\partial z_0}, \\ \displaystyle N_{c_{rr}2} = u_{r,1} \frac{\partial c_{rr,1}}{\partial r} + u_{z,1} \frac{\partial c_{rr,1}}{\partial z_0} - 2 c_{rr,1} \frac{\partial u_{r,1}}{\partial r} - 2 c_{rz,1} \frac{\partial u_{r,1}}{\partial z_0}, \\ N_{c_{rz}2} = u_{r,1} \frac{\partial c_{rz,1}}{\partial r} + u_{z,1} \frac{\partial c_{rz,1}}{\partial z_0} - c_{rr,1} \frac{\partial u_{z,1}}{\partial r} - c_{rz,1} \frac{\partial u_{z,1}}{\partial z_0} - c_{rz,1} \frac{\partial u_{r,1}}{\partial r} - c_{zz,1} \frac{\partial u_{r,1}}{\partial z_0},\\ N_{c_{\theta \theta}2} = u_{r,1} \frac{\partial c_{\theta \theta,1}}{\partial r} + u_{z,1} \frac{\partial c_{\theta \theta,1}}{\partial z_0} - \frac{2 u_{r,1}}{r} c_{\theta \theta,1},\\ N_{c_{zz}2} = u_{r,1} \frac{\partial c_{zz,1}}{\partial r} + u_{z,1} \frac{\partial c_{zz,1}}{\partial z_0} - 2 c_{rz,1} \frac{\partial u_{z,1}}{\partial r} - 2 c_{zz,1} \frac{\partial u_{z,1}}{\partial z_0}, \end{gathered}\right\} \end{equation}

and

(C12) \begin{align} \left.\begin{aligned} N_{u_{r}3} &= u_{r,1} \frac{\partial u_{r,2}}{\partial r} + u_{r,2} \frac{\partial u_{r,1}}{\partial r} + u_{z,1} \frac{\partial u_{r,2}}{\partial z_0} + u_{z,1} \frac{\partial u_{r,1}}{\partial z_1} + u_{z,2} \frac{\partial u_{r,1}}{\partial z_0},\\ N_{u_{z}3} &= u_{r,1} \frac{\partial u_{z,2}}{\partial r} + u_{r,2} \frac{\partial u_{z,1}}{\partial r} + u_{z,1} \frac{\partial u_{z,2}}{\partial z_0} + u_{z,1} \frac{\partial u_{z,1}}{\partial z_1} + u_{z,2} \frac{\partial u_{z,1}}{\partial z_0},\\ N_{c_{rr}3} &= u_{r,1} \frac{\partial c_{rr,2}}{\partial r} +u_{r,2} \frac{\partial c_{rr,1}}{\partial r} + u_{z,1} \frac{\partial c_{rr,2}}{\partial z_0} + u_{z,1} \frac{\partial c_{rr,1}}{\partial z_1} + u_{z,2} \frac{\partial c_{rr,1}}{\partial z_0}\\ &\quad - 2 c_{rr,1} \frac{\partial u_{r,2}}{\partial r} - 2 c_{rr,2} \frac{\partial u_{r,1}}{\partial r} - 2 c_{rz,1} \frac{\partial u_{r,2}}{\partial z_0} - 2 c_{rz,1} \frac{\partial u_{r,1}}{\partial z_1} - 2 c_{rz,2} \frac{\partial u_{r,1}}{\partial z_0},\\ N_{c_{rz}3} &= u_{r,1} \frac{\partial c_{rz,2}}{\partial r} +u_{r,2} \frac{\partial c_{rz,1}}{\partial r} + u_{z,1} \frac{\partial c_{rz,2}}{\partial z_0} + u_{z,1} \frac{\partial c_{rz,1}}{\partial z_1} + u_{z,2} \frac{\partial c_{rz,1}}{\partial z_0}\\ &\quad - c_{rr,1} \frac{\partial u_{z,2}}{\partial r} - c_{rr,2} \frac{\partial u_{z,1}}{\partial r} - c_{rz,1} \frac{\partial u_{z,2}}{\partial z_0} - c_{rz,1} \frac{\partial u_{z,1}}{\partial z_1} - c_{rz,2} \frac{\partial u_{z,1}}{\partial z_0} \\ &\quad - c_{rz,1} \frac{\partial u_{r,2}}{\partial r} - c_{rz,2} \frac{\partial u_{r,1}}{\partial r} - c_{zz,1} \frac{\partial u_{r,2}}{\partial z_0} - c_{zz,1} \frac{\partial u_{r,1}}{\partial z_1} - c_{zz,2} \frac{\partial u_{r,1}}{\partial z_0}, \\ N_{c_{\theta \theta}3} &= u_{r,1} \frac{\partial c_{\theta \theta,2}}{\partial r} +u_{r,2} \frac{\partial c_{\theta \theta,1}}{\partial r}+ u_{z,1} \frac{\partial c_{\theta \theta,2}}{\partial z_0} + u_{z,1} \frac{\partial c_{\theta \theta,1}}{\partial z_1} + u_{z,2} \frac{\partial c_{\theta \theta,1}}{\partial z_0}\\ &\quad - \frac{2 }{r} c_{\theta \theta,1} u_{r,2} - \frac{2 }{r} c_{\theta \theta,2} u_{r,1},\\ N_{c_{zz}3} &= u_{r,1} \frac{\partial c_{zz,2}}{\partial r} +u_{r,2} \frac{\partial c_{zz,1}}{\partial r} + u_{z,1} \frac{\partial c_{zz,2}}{\partial z_0} + u_{z,1} \frac{\partial c_{zz,1}}{\partial z_1} + u_{z,2} \frac{\partial c_{zz,1}}{\partial z_0}\\ &\quad - 2 c_{rz,1} \frac{\partial u_{z,2}}{\partial r} - 2 c_{rz,2} \frac{\partial u_{z,1}}{\partial r} - 2 c_{zz,1} \frac{\partial u_{z,2}}{\partial z_0} - 2 c_{zz,1} \frac{\partial u_{z,1}}{\partial z_1} - 2 c_{zz,2} \frac{\partial u_{z,1}}{\partial z_0}, \end{aligned}\right\} \end{align}

and $\tilde{\boldsymbol {N}}_{3}^{\circ }$ in (2.25) can be obtained by singling out the terms of the linear wave $e^{i\alpha z_0 + \mu t_0}$ and then moving the $A$-related coefficient out of the original expression, i.e.  in form of $|A|^2 A \tilde{\boldsymbol {N}}_{3}^{\circ } \exp ({\textrm {i}\alpha z_0 + \mu t_0})$.

C.3. Adjoint operators

For the $up$-$c$ formulation, we take the coefficient matrix in the inner product (2.17) to be the identity matrix $\boldsymbol {W}=\boldsymbol {I}$ and the adjoint operators can be then derived from (2.21). As a result, $\boldsymbol {M}_0=\boldsymbol {M}_0^{\dagger}$ is self-adjoint and $\boldsymbol {L}_0^{\dagger}$ is given as

(C13)\begin{equation} \boldsymbol{L}_0^{\dagger}= \begin{pmatrix} L_{r,r0}^{\dagger} & L_{r,z0}^{\dagger} & L_{r,p0}^{\dagger} & L_{r,rr0}^{\dagger} & L_{r,rz0}^{\dagger} & L_{r,\theta\theta 0}^{\dagger} & L_{r,zz0}^{\dagger} \\ 0 & L_{z,z0}^{\dagger} & L_{z,p0}^{\dagger} & 0 & L_{z,rz0}^{\dagger} & 0 & L_{z,zz0}^{\dagger} \\ L_{p,r0}^{\dagger} & L_{p,z0}^{\dagger} & 0 & 0 & 0 & 0 & 0 \\ L_{rr,r0}^{\dagger} & 0 & 0 & L_{rr,rr0}^{\dagger} & L_{rr,rz0}^{\dagger} & 0 & 0\\ L_{rz,r0}^{\dagger} & L_{rz,z0}^{\dagger} & 0 & 0 & L_{rz,rz0}^{\dagger} & 0 & L_{rz,zz0}^{\dagger}\\ L_{\theta\theta,r0}^{\dagger} & 0 & 0 & 0 & 0 & L_{\theta\theta,\theta\theta 0}^{\dagger} & 0 \\ 0 & L_{zz,z0}^{\dagger} & 0 & 0 & 0 & 0 & L_{zz,zz0}^{\dagger} \\ \end{pmatrix}, \end{equation}

where the non-zero elements are

(C14) \begin{equation} \left.\begin{gathered} L_{r,r0}^{\dagger} = U_{z} \frac{\partial }{\partial z_0} + \frac{\beta}{Re_c} \left( \frac{\partial^{2}}{\partial r^{2}} + \frac{1}{r}\frac{\partial}{\partial r} + \frac{\partial^{2}}{\partial z_0^{2}} - \frac{1}{r^2} \right),\quad L_{r,z0}^{\dagger}={-}U_{z}', \\ L_{r,p0}^{\dagger}={-}\frac{\partial }{\partial r},\quad L_{r,rr0}^{\dagger}={-}2\frac{\partial }{\partial r} - \frac{2}{r}- 2 C_{rz}\frac{\partial }{\partial z_0},\\ L_{r,rz0}^{\dagger}={-} 2 C_{rz}' - C_{rz}\frac{\partial }{\partial r} - \frac{C_{rz}}{ r} - C_{zz}\frac{\partial }{\partial z_0},\quad L_{r,\theta\theta 0}^{\dagger}= \frac{2}{r}, \quad L_{r,zz0}^{\dagger}={-} C_{zz}', \\ L_{z,z0}^{\dagger}= U_{z} \frac{\partial }{\partial z_0} + \frac{\beta}{Re_c} \left( \frac{\partial^{2}}{\partial r^{2}} + \frac{1}{r}\frac{\partial}{\partial r} + \frac{\partial^{2}}{\partial z_0^{2}} \right),\quad L_{z,p0}^{\dagger}={-}\frac{\partial }{\partial z_0}, \\ L_{z,rz0}^{\dagger}={-}\frac{\partial }{\partial r}-\frac{1}{r} - C_{rz}\frac{\partial }{\partial z_0}, \\ L_{z,zz0}^{\dagger}={-}2 C_{rz}\frac{\partial }{\partial r} -2 C_{rz}'- \frac{2 C_{rz}}{ r} - 2 C_{zz}\frac{\partial }{\partial z_0},\quad L_{p,r0}^{\dagger}=\frac{\partial}{\partial r} + \frac{1}{ r}, \\ L_{p,z0}^{\dagger}= \frac{\partial}{\partial z_0},\quad L_{rr,r0}^{\dagger}=\frac{1-\beta}{Re_c Wi} \left( -\frac{\partial }{\partial r} \right), \\ L_{rr,rr0}^{\dagger}= U_{z} \frac{\partial }{\partial z_0} - \frac{1}{Wi},\quad L_{rr,rz0}^{\dagger}= U_{z}', \\ L_{rz,r0}^{\dagger}=\frac{1-\beta}{Re_c Wi} \left( -\frac{\partial }{\partial z_0} \right),\quad L_{rz,z0}^{\dagger}=\frac{1-\beta}{Re_c Wi} \left( -\frac{\partial }{\partial r} \right), \\ L_{rz,rz0}^{\dagger}= U_{z} \frac{\partial }{\partial z_0} - \frac{1}{Wi},\quad L_{rz,zz0}^{\dagger}= 2 U_{z}', \\ L_{\theta\theta,r0}^{\dagger}=\frac{1-\beta}{Re_c Wi} \left(- \frac{1}{r} \right),\quad L_{\theta\theta,\theta\theta 0}^{\dagger}=U_{z} \frac{\partial }{\partial z_0} - \frac{1}{Wi}, \\ L_{zz,z0}^{\dagger}=\frac{1-\beta}{Re_c Wi} \left( -\frac{\partial }{\partial z_0} \right),\quad L_{zz,zz0}^{\dagger}=U_{z} \frac{\partial }{\partial z_0} - \frac{1}{Wi}. \end{gathered}\right\} \end{equation}

The boundary conditions for the adjoint problem are the same as those for the direct problem.

Appendix D. Validation of the Landau coefficients using DNS

D.1. Numerical methods for the DNS

We have conducted DNS to evaluate the Landau coefficient and validate our weakly nonlinear results. To simulate the pressure-driven viscoelastic pipe flow, an additional pressure gradient term $4(1+\alpha _f)/Re \, \hat {\boldsymbol {e}}_z$ is added to the right-hand side of the Navier–Stokes equation (2.1a,b), where $\alpha _f$ is a fluctuating coefficient adjusted at every time step to maintain a constant flux in the flow. The governing equation of the conformation tensor remains the same as in (2.2). This equation system is solved by a numerical solver developed based on an open-source code openpipeflow (Willis Reference Willis2017). A Fourier–Chebyshev pseudo-spectral method and a high-order finite-difference method is adopted in the homogeneous (azimuthal and axial) directions and in the radial direction, respectively. The no-slip boundary condition for velocity and the incompressibility condition for the flow field are enforced by a pressure-Poisson equation formulation and an influence-matrix technique is applied to enforce continuity. A second-order predictor–corrector scheme based on the Crank–Nicolson method is employed for the time-stepping. All the DNS are performed with an initial condition constructed with the eigenfunction obtained from the corresponding linear eigenvalue problem. The radial resolution is set by $N_r$ being the number of nodes in the radial direction and the time-step size $\delta t$ can be fixed or adaptive subject to a Courant–Friedrichs–Lewy number of 0.25. A more detailed description of this code has been reported in Sun, Wan & Zhang (Reference Sun, Wan and Zhang2021).

D.2. Convergence check of the DNS results

To provide an additional validation, we use the DNS method to investigate the evolution of a small disturbance, which starts from the linear phase, undergoes a weakly nonlinear phase and then saturates. As there is no diffusion term in the original governing equation of the conformation tensor, a sharp gradient can cause numerical instability. A widely adopted strategy to solve this problem is to add a small amount of global artificial diffusion to the equation (see Sureshkumar & Beris Reference Sureshkumar and Beris1995a; Lopez et al. Reference Lopez, Choueiri and Hof2019). However, we note that even a quite small amount of artificial diffusion can alter the linear and weakly nonlinear phases and, thus, we choose to not add any artificial diffusion. This strategy is feasible because only a few modes need to be resolved for the evaluation of the Landau coefficient $a_{3,r}$ (because the multiple-scale expansion is based on several linear and nonlinear modes and we do not need to resolve all the nonlinearities). In this case, increased spatial resolution in the radial direction and a smaller time-step size are needed to be applied: we find that the resolution of $N_k=2$ (modes in streamwise direction), $N_r=256$ and $\delta t=0.001$ suffices to accurately resolve the linear phase at $Re=800$, $\alpha =1$, $Wi=30$ and $\beta =0.65$ (in this case, the growth rate in the linear phase evaluated using DNS is $-0.01420$, which agrees well with the linear analysis prediction value $-0.01425$; and the relative error is about 0.35 %). The number of Fourier modes in the axial direction $N_k$ has also been examined as shown in figure 15. We found that inclusion of more modes indeed makes the computation more difficult. However, only the initial derivation from the linear phase is necessary to evaluate the Landau coefficient $a_{3,r}$. The result indicates that the essential nonlinearities can be faithfully resolved with $N_{k}\geqslant 4$.

Figure 15.

Evolution of disturbance kinetic energy (normalized by the kinetic energy of the corresponding Newtonian base flow $E_{k,b}$) at: (a) $Re=280$, $\alpha _{c}=0.515525$, $Wi=65$ and $\beta =0.65$; (b) $Re=90$, $\alpha _{c}=1.550722$, $Wi=65$ and $\beta =0.90$. These results are obtained from DNS and $N_k$ is the number of Fourier modes in the axial direction. Generally, $N_{k}\geqslant 4$ suffices to capture the weak nonlinearity for the evaluation of the Landau coefficients.

D.3. Comparison of the DNS results with weakly nonlinear predictions

As we are mainly interested in the global amplification of disturbance amplitudes, we will consider the following simplified Landau model:

(D1)\begin{equation} \frac{\textrm{d} |A|}{\textrm{d} t} = c_{1} |A| + c_{3} |A|^{3}. \end{equation}

Here $|A|$ is the amplitude of the perturbation, $c_{1}$ is the linear growth rate of the amplitude and $c_{3}$ determines the bifurcation type. The original Landau equation is in a series form and the truncation up to the third order is consistent with our theoretical derivation. In order to compare quantitatively with our weakly nonlinear analysis results, the global amplitude $|A|$ in DNS is defined based on the total energy of the disturbances contained in the mode ($k=1,n=0$, i.e.  the first axial Fourier mode in axisymmetric flows) as

(D2)\begin{equation} |A| = \sqrt{\frac{1}{2} \int_{0}^{1} \left( \left(|\tilde{u}_r|^2 + |\tilde{u}_z|^2 \right) + \frac{1-\beta}{Re Wi} \left( |\tilde{g}_{rr}|^{2}+2|\tilde{g}_{rz}|^{2}+|\tilde{g}_{\theta\theta }|^{2}+|\tilde{g}_{zz}|^2\right) \right) r\,\textrm{d}r}. \end{equation}

Based on the convergence check of the kinetic energy evolution curve, $N_k=5$ is applied in DNS for the evaluation of the Landau coefficient. The evolution of $|A|$ as defined in (D 2) is recorded as shown in figure 16(a) where $|A|$ grows and finally saturates. From the log plot shown in figure 16(b), an apparent linear growth stage can be observed followed by a nonlinear saturation phase. The Landau coefficients $c_{1}$ and $c_{3}$ can be calculated by recasting (D 1) in the form of ${\textrm {d} (\log |A|)}/{\textrm {d} t} = c_{1} + c_{3} |A|^{2}$, and plotting $\textrm {d} (\log |A|)/\textrm {d}t$ vs $|A|^{2}$ so that the vertical intercept point gives an estimation of $c_{1}$ and the gradient near this point is an approximation of $c_{3}$ (Thompson, Leweke & Provansal Reference Thompson, Leweke and Provansal2001). Such a plot is shown in figure 16(c) and the Landau coefficients can be determined as $c_{1} \approx 3.05 \times 10^{-4}$ and $c_{3} \approx -168$. It should be noted that the GLE (2.24) derived from the weakly nonlinear analysis is in a $t_{2}$ time scale, while the simplified Landau equation (D 1) used for the DNS data is in a $t$ time scale, and the coefficients are related by $c_{1}=a_{1}(Re-Re_{c})$ and $c_{3}=a_{3,r}$. For the present case with $Re=280$ and $Re_{c}=265.572335$, $a_{1}(Re-Re_{c}) \approx 3.13\times 10^{-4}$, which is close to $c_{1}$ with a relative error of about $2.6\,\%$; $a_{3,r} \approx -154$ is in agreement with $c_{3}$ with a relative error of about $9.1\,\%$ (see table 3 for the value of $a_{1}$ and also $a_{3,r}$). The sources of the error can be: (1) few $N_k$ modes, and (2) relative distance to the linear criticality condition. A further DNS test at $Re=270$ being closer to the critical value $Re_{c}=265.572335$ with other parameters unchanged, gives $c_{3} \approx -164$, which is also closer to the theoretical prediction $a_{3,r} \approx -154$ with a smaller relative error of $6.5\,\%$. If $Re$ is too close to the critical value, the exponential-growth phase will take a much long time. In addition to the supercritical cases, the weakly nonlinear analysis of a subcritical case at $Re_c=84.644523$, $\alpha _{c}=1.550722$, $Wi=65$, $\beta =0.90$ predicts $a_{3,r} \approx 1179$, which also agrees with the DNS results $c_{3} \approx 1241$ obtained at $Re=90$ with a relative error of $5.3\,\%$. For this subcritical case, although the DNS finally blows up as shown in figure 15(b), the weak nonlinearity relating to the initial derivation can still be captured, suggesting that the simplification by using few axial modes in DNS is reasonable. The DNS results double-check our theoretical prediction and verify the bifurcation types in different conditions.

Figure 16.

Evolution and post-processing of disturbance amplitude $|A|$ at $Re=280$, $\alpha _{c}=0.515525$, $Wi=65$ and $\beta =0.65$ obtained with $N_r=400$ and $N_k=5$ from the DNS: (a) original signal; (b) log plot of (a) demonstrating the linear growth and nonlinear saturation; (c) time derivative of $\log |A|$ vs $|A|^{2}$ with the vertical intercept giving the linear growth rate $c_{1} \approx 3.05 \times 10^{-4}$ and the gradient near this intercept point giving the Landau coefficient $c_{3} \approx -168$.

Appendix E. Linear instability in the UCM limit

In this appendix we provide more results on the viscoelastic pipe flows of UCM fluids. The eigenspectra as shown in figure 5 have a typical structure consisting of high-frequency Gorodtsov-Leonov (HFGL) modes (Gorodtsov & Leonov Reference Gorodtsov and Leonov1967), a continuous spectrum in a balloon-like shape due to finite spatial resolutions (Chaudhary et al. Reference Chaudhary, Garg, Shankar and Subramanian2019), and some discrete modes. Such a structure agrees well with that described in Chaudhary et al. (Reference Chaudhary, Garg, Subramanian and Shankar2021) for the eigenspectrum of pipe Poiseuille flows with UCM fluids. Theoretically, the HFGL modes should distribute around a horizontal line $\omega _i/\alpha = - 1/(2\alpha Wi)$, which predicts $\omega _i=-0.008$ in good agreement with that shown in panel figure 5(a) where $Re=1000$, $\alpha =0.2$ and $Wi=60$, and $\omega _i=-0.034$ in panel $(b)$ where $Re=260$, $\alpha =0.94$ and $Wi=14.5$. In addition, the continuous spectrum should locate around another horizontal line $\omega _i/\alpha = - 1/(\alpha Wi)$; and this theoretical prediction also agrees with the obtained continuous spectra where $\omega _i=-0.017$ for panel (a) and $\omega _i=-0.069$ for panel $(b)$. In addition to the convergence of the eigenspectra, the complex wave speed $c=\omega /\alpha$ is shown to be continuous with decreasing $\beta$ even to zero as plotted in figure 17. For the three sets of parameters explored, linear instability can be observed at small $\beta$ and it continues to exist when $\beta =0$ which is the UCM limit. Similar to the linear instability found by Garg et al. (Reference Garg, Chaudhary, Khalid, Shankar and Subramanian2018), this mode is also a centre mode. The eigenfunctions of the unstable mode presented in figure 5(a) are plotted in figure 18. The eigenfunctions are not localized near the pipe axis $r=0$. Instead, they go through the whole cross-section of the pipe and seem to be somehow like the eigenfunctions of wall modes which have large variations near the wall (see the variation of $\tilde {u}_{z,1}$ and $\tilde {c}_{zz,1}$).

Figure 17.

Continuous variation of $c=\omega /\alpha$ with viscosity ratio $\beta$ at three sets of parameters: (a) imaginary part $c_i$; (b) real part $c_r$. Linear instability exists ($c_i>0$) in the range of $\beta$ above the horizontal black line in (a).

Figure 18.

Eigenfunctions corresponding with the unstable mode shown in figure 5(a) at Reynolds number $Re=1000$, wavenumber $\alpha =0.2$ and Weissenberg number $Wi=60$: (a) velocity components $\tilde {u}_{z,1}$ and $\tilde {u}_{r,1}$; (b) conformation tensor components $\tilde {c}_{rr,1}$ and $\tilde {c}_{rz,1}$; (c) conformation tensor components $\tilde {c}_{\theta \theta,1}$ and $\tilde {c}_{zz,1}$.

Figure 19.

Neutral curves for pipe flow of UCM fluids: $(a)$ at varying Weissenberg numbers $Wi$; (b) at $Wi=19$ with the colourbar indicating the linear growth rate $\mu _r$. The two arrows in (b) point at the two critical points.

To examine the variation of instability with changing $Wi$, several neutral curves in the $(\alpha,Re)$ plane are plotted in figure 19. As expected, these neutral curves are also in the form of closed loops like those with finite viscosity ratios (see figure 3a), but in more complex shapes with multiple loops. At higher $Wi$, the linear instability is limited in a very narrow range of low wavenumbers, covering a broad range of $Re$; and the loop moves towards small $Re$ and high wavenumber as $Wi$ decreases. Notably, when $Wi$ is decreased, the loops may split and disappear. Such complex behaviour causes the discontinuity of the critical parameters and Landau coefficients at $Wi=19$ as described in § 4. Finally, at about $Wi=14.5$, the loop shrinks to almost a very small area (at $Re\approx 255$) and then the linear instability completely disappears (at least in our search of the linear instability for the current parameters). This type of behaviour of neutral curves (being closed, disconnected and may split) has been reported in other complex fluid systems, for instance, in stably stratified plane Poiseuille flows (Fujimura & Kelly Reference Fujimura and Kelly1997), convection in a multicomponent fluid layer (Terrones & Pearlstein Reference Terrones and Pearlstein1989) and the parallel shear flow between two inclined plates with different temperatures (Chen & Pearlstein Reference Chen and Pearlstein1989). In the present case, we can see that polymer elasticity (characterized by $Wi$) plays an important role in stabilizing UCM pipe flows, in contrast to the destabilizing effect in UCM channel flows (Porteous & Denn Reference Porteous and Denn1972; Chaudhary et al. Reference Chaudhary, Garg, Shankar and Subramanian2019). The lower bound of linear instability (reached at about $Wi=14.5$) aforementioned shows that the UCM pipe flow is linearly stable at low $Re$, analogous to the absence of linear instability at low $Re$ in UCM channel flows reported in some studies in early works (Ho & Denn Reference Ho and Denn1977; Lee & Finlayson Reference Lee and Finlayson1986).

References

REFERENCES

Bender, C.M. & Orszag, S.A. 1999 Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory. Springer.CrossRefGoogle Scholar
Bertola, V., Meulenbroek, B., Wagner, C., Storm, C., Morozov, A., van Saarloos, W. & Bonn, D. 2003 Experimental evidence for an intrinsic route to polymer melt fracture phenomena: a nonlinear instability of viscoelastic Poiseuille flow. Phys. Rev. Lett. 90 (11), 114502.CrossRefGoogle ScholarPubMed
Bird, R.B., Curtiss, C.F., Armstrong, R.C. & Hassager, O. 1987 Dynamics of Polymeric Liquids, Volume 2: Kinetic Theory. Wiley.Google Scholar
Bonn, D., Ingremeau, F., Amarouchene, Y. & Kellay, H. 2011 Large velocity fluctuations in small-Reynolds-number pipe flow of polymer solutions. Phys. Rev. E 84 (4), 045301.CrossRefGoogle ScholarPubMed
Bouteraa, M. & Nouar, C. 2015 Weakly nonlinear analysis of Rayleigh–Bénard convection in a non-Newtonian fluid between plates of finite conductivity: influence of shear-thinning effects. Phys. Rev. E 92 (6), 063017.CrossRefGoogle Scholar
Bouteraa, M., Nouar, C., Plaut, E., Métivier, C. & Kalck, A. 2015 Weakly nonlinear analysis of Rayleigh–Bénard convection in shear-thinning fluids: nature of the bifurcation and pattern selection. J. Fluid Mech. 767, 696–734.CrossRefGoogle Scholar
Buza, G., Page, J. & Kerswell, R.R. 2021 Weakly nonlinear analysis of the viscoelastic instability in channel flow for finite and vanishing Reynolds numbers. arXiv:2107.06191v1.Google Scholar
Chandra, B., Shankar, V. & Das, D. 2018 Onset of transition in the flow of polymer solutions through microtubes. J. Fluid Mech. 844, 1052–1083.CrossRefGoogle Scholar
Chandra, B., Shankar, V. & Das, D. 2020 Early transition, relaminarization and drag reduction in the flow of polymer solutions through microtubes. J. Fluid Mech. 885, A47.CrossRefGoogle Scholar
Chaudhary, I., Garg, P., Shankar, V. & Subramanian, G. 2019 Elasto-inertial wall mode instabilities in viscoelastic plane Poiseuille flow. J. Fluid Mech. 881, 119–163.CrossRefGoogle Scholar
Chaudhary, I., Garg, P., Subramanian, G. & Shankar, V. 2021 Linear instability of viscoelastic pipe flow. J. Fluid Mech. 908, A11.CrossRefGoogle Scholar
Chen, Y.-M. & Pearlstein, A.J. 1989 Stability of free-convection flows of variable-viscosity fluids in vertical and inclined slots. J. Fluid Mech. 198, 513–541.CrossRefGoogle Scholar
Chokshi, P. & Kumaran, V. 2009 Stability of the plane shear flow of dilute polymeric solutions. Phys. Fluids 21 (1), 014109.CrossRefGoogle Scholar
Choueiri, G.H., Lopez, J.M. & Hof, B. 2018 Exceeding the asymptotic limit of polymer drag reduction. Phys. Rev. Lett. 120, 124501.CrossRefGoogle ScholarPubMed
Choueiri, G.H., Lopez, J.M., Varshney, A., Sankar, S. & Hof, B. 2021 Experimental observation of the origin and structure of elasto-inertial turbulence. arXiv:2103.00023.Google Scholar
Cross, M.C. & Hohenberg, P.C. 1993 Pattern formation outside of equilibrium. Rev. Mod. Phys. 65, 851–1112.CrossRefGoogle Scholar
Davey, A. 1978 On Itoh's finite amplitude stability theory for pipe flow. J. Fluid Mech. 86 (4), 695–703.CrossRefGoogle Scholar
Davey, A. & Drazin, P.G. 1969 The stability of Poiseuille flow in a pipe. J. Fluid Mech. 36 (2), 209–218.CrossRefGoogle Scholar
Davey, A. & Nguyen, H. 1971 Finite-amplitude stability of pipe flow. J. Fluid Mech. 45 (4), 701–720.CrossRefGoogle Scholar
Dubief, Y., Terrapon, V.E. & Soria, J. 2013 On the mechanism of elasto-inertial turbulence. Phys. Fluids 25 (11), 110817.CrossRefGoogle ScholarPubMed
Forame, P.C., Hansen, R.J. & Little, R.C. 1972 Observations of early turbulence in the pipe flow of drag reducing polymer solutions. AIChE J. 18 (1), 213–217.CrossRefGoogle Scholar
Fujimura, K. 1989 The equivalence between two perturbation methods in weakly nonlinear stability theory for parallel shear flows. Proc. R. Soc. Lond. A 424 (1867), 373–392.Google Scholar
Fujimura, K. & Kelly, R. 1997 Degenerate bifurcation in stably stratified plane Poiseuille flow. J. Fluid Mech. 331, 261–282.CrossRefGoogle Scholar
Gao, Z., Sergent, A., Podvin, B., Xin, S., Le Quéré, P. & Tuckerman, L.S. 2013 Transition to chaos of natural convection between two infinite differentially heated vertical plates. Phys. Rev. E 88 (2), 023010.CrossRefGoogle ScholarPubMed
Garg, P., Chaudhary, I., Khalid, M., Shankar, V. & Subramanian, G. 2018 Viscoelastic pipe flow is linearly unstable. Phys. Rev. Lett. 121 (2), 024502.CrossRefGoogle ScholarPubMed
Gorodtsov, V.A. & Leonov, A.I. 1967 On a linear instability of a plane parallel Couette flow of viscoelastic fluid. J. Appl. Maths Mech. 31 (2), 310–319.CrossRefGoogle Scholar
Graham, M. 1998 Effect of axial flow on viscoelastic Taylor–Couette instability. J. Fluid Mech. 360, 341–374.CrossRefGoogle Scholar
Graham, M.D. 2014 Drag reduction and the dynamics of turbulence in simple and complex fluids. Phys. Fluids 26 (10), 625–656.CrossRefGoogle Scholar
Hameduddin, I., Gayme, D.F. & Zaki, T.A. 2019 Perturbative expansions of the conformation tensor in viscoelastic flows. J. Fluid Mech. 858, 377–406.CrossRefGoogle Scholar
Hameduddin, I., Meneveau, C., Zaki, T.A. & Gayme, D.F. 2018 Geometric decomposition of the conformation tensor in viscoelastic turbulence. J. Fluid Mech. 842, 395–427.CrossRefGoogle Scholar
Hameduddin, I. & Zaki, T.A. 2019 The mean conformation tensor in viscoelastic turbulence. J. Fluid Mech. 865, 363–380.CrossRefGoogle Scholar
Hansen, R. 1973 Stability of laminar pipe flows of drag reducing polymer solutions in the presence of high-phase-velocity disturbances. AIChE J. 19 (2), 298–304.CrossRefGoogle Scholar
Herbert, T. 1980 Nonlinear stability of parallel flows by high-ordered amplitude expansions. AIAA J. 18 (3), 243–248.CrossRefGoogle Scholar
Herbert, T. 1983 On perturbation methods in nonlinear stability theory. J. Fluid Mech. 126, 167–186.CrossRefGoogle Scholar
Ho, T.C. & Denn, M.M. 1977 Stability of plane Poiseuille flow of a highly elastic liquid. J. Non-Newtonian Fluid Mech. 3 (2), 179–195.CrossRefGoogle Scholar
Huerre, P. & Rossi, M. 1998 Hydrodynamics and nonlinear instabilities. In Hydrodynamic Instabilities in Open Flows (ed. C. Godrèche & P. Manneville), pp. 81–294. Cambridge University Press.CrossRefGoogle Scholar
Itoh, N. 1977 Nonlinear stability of parallel flows with subcritical Reynolds numbers. Part 2. Stability of pipe Poiseuille flow to finite axisymmetric disturbances. J. Fluid Mech. 82 (3), 469–479.CrossRefGoogle Scholar
Khalid, M., Chaudhary, I., Garg, P., Shankar, V. & Subramanian, G. 2021 The centre-mode instability of viscoelastic plane Poiseuille flow. J. Fluid Mech. 915, A43.CrossRefGoogle Scholar
Kolyshkin, A.A. & Ghidaoui, M.S. 2003 Stability analysis of shallow wake flows. J. Fluid Mech. 494, 355–377.CrossRefGoogle Scholar
Landau, L.D. 1944 On the problem of turbulence. C. R. Acad. Sci. URSS 44 (31), 1–314.Google Scholar
Lee, K.-C. & Finlayson, B.A. 1986 Stability of plane Poiseuille and Couette flow of a Maxwell fluid. J. Non-Newtonian Fluid Mech. 21 (1), 65–78.CrossRefGoogle Scholar
Lopez, J.M., Choueiri, G.H. & Hof, B. 2019 Dynamics of viscoelastic pipe flow at low Reynolds numbers in the maximum drag reduction limit. J. Fluid Mech. 874, 699–719.CrossRefGoogle Scholar
Luchini, P. & Bottaro, A. 2014 Adjoint equations in stability analysis. Annu. Rev. Fluid Mech. 46 (1), 493–517.CrossRefGoogle Scholar
Meseguer, A. & Trefethen, L.N. 2003 Linearized pipe flow to Reynolds number $10^{7}$. J. Comput. Phys. 186 (1), 178–197.CrossRefGoogle Scholar
Meulenbroek, B., Storm, C., Bertola, V., Wagner, C., Bonn, D. & van Saarloos, W. 2003 Intrinsic route to melt fracture in polymer extrusion: a weakly nonlinear subcritical instability of viscoelastic Poiseuille flow. Phys. Rev. Lett. 90 (2), 024502.CrossRefGoogle ScholarPubMed
Meulenbroek, B., Storm, C., Morozov, A.N. & van Saarloos, W. 2004 Weakly nonlinear subcritical instability of visco-elastic Poiseuille flow. J. Non-Newtonian Fluid Mech. 116 (2–3), 235–268.CrossRefGoogle Scholar
Mohseni, K. & Colonius, T. 2000 Numerical treatment of polar coordinate singularities. J. Comput. Phys. 157 (2), 787–795.CrossRefGoogle Scholar
Morozov, A. & van Saarloos, W. 2019 Subcritical instabilities in plane Poiseuille flow of an Oldroyd-B fluid. J. Stat. Phys. 175 (3), 554–577.CrossRefGoogle Scholar
Morozov, A.N. & van Saarloos, W. 2005 Subcritical finite-amplitude solutions for plane Couette flow of viscoelastic fluids. Phys. Rev. Lett. 95, 024501.CrossRefGoogle ScholarPubMed
Morozov, A.N. & van Saarloos, W. 2007 An introductory essay on subcritical instabilities and the transition to turbulence in visco-elastic parallel shear flows. Phys. Rep. 447 (3), 112–143.CrossRefGoogle Scholar
Orlandi, P. 2012 Fluid Flow Phenomena: A Numerical Toolkit. Springer Science & Business Media.Google Scholar
Page, J., Dubief, Y. & Kerswell, R.R. 2020 Exact traveling wave solutions in viscoelastic channel flow. Phys. Rev. Lett. 125, 154501.CrossRefGoogle ScholarPubMed
Patera, A.T. & Orszag, S.A. 1981 Finite-amplitude stability of axisymmetric pipe flow. J. Fluid Mech. 112, 467–474.CrossRefGoogle Scholar
Porteous, K.C. & Denn, M.M. 1972 Linear stability of plane Poiseuille flow of viscoelastic liquids. Trans. Soc. Rheol. 16 (2), 295–308.CrossRefGoogle Scholar
Ram, A. & Tamir, A. 1964 Structural turbulence in polymer solutions. J. Appl. Polym. Sci. 8 (6), 2751–2762.CrossRefGoogle Scholar
Ramanan, V., Kumar, K. & Graham, M. 1999 Stability of viscoelastic shear flows subjected to steady or oscillatory transverse flow. J. Fluid Mech. 379, 255–277.CrossRefGoogle Scholar
Renardy, M. 1992 A rigorous stability proof for plane Couette flow of an upper convected Maxwell fluid at zero Reynolds number. Eur. J. Mech. (B/Fluids) 11 (4), 511–516.Google Scholar
Reynolds, W. & Potter, M.C. 1967 Finite-amplitude instability of parallel shear flows. J. Fluid Mech. 27 (3), 465–492.CrossRefGoogle Scholar
Sadanandan, B. & Sureshkumar, R. 2002 Viscoelastic effects on the stability of wall-bounded shear flows. Phys. Fluids 14 (1), 41–48.CrossRefGoogle Scholar
Samanta, D., Dubief, Y., Holzner, M., Schäfer, C., Morozov, A.N., Wagner, C. & Hof, B. 2013 Elasto-inertial turbulence. Proc. Natl Acad. Sci. USA 110 (26), 10557–10562.CrossRefGoogle ScholarPubMed
Sen, P.K. & Venkateswarlu, D. 1983 On the stability of plane Poiseuille flow to finite-amplitude disturbances, considering the higher-order Landau coefficients. J. Fluid Mech. 133, 179–206.CrossRefGoogle Scholar
Shaqfeh, E.S. 1996 Purely elastic instabilities in viscometric flows. Annu. Rev. Fluid Mech. 28, 129–185.CrossRefGoogle Scholar
Shekar, A., McMullen, R.M., McKeon, B.J. & Graham, M.D. 2020 Self-sustained elastoinertial Tollmien–Schlichting waves. J. Fluid Mech. 897, A3.CrossRefGoogle Scholar
Shekar, A., McMullen, R.M., McKeon, B.J. & Graham, M.D. 2021 Tollmien–Schlichting route to elastoinertial turbulence in channel flow. arXiv:2104.10257.CrossRefGoogle Scholar
Shekar, A., McMullen, R.M., Wang, S.-N., McKeon, B.J. & Graham, M.D. 2019 Critical-layer structures and mechanisms in elastoinertial turbulence. Phys. Rev. Lett. 122, 124503.CrossRefGoogle ScholarPubMed
Sid, S., Terrapon, V.E. & Dubief, Y. 2018 Two-dimensional dynamics of elasto-inertial turbulence and its role in polymer drag reduction. Phys. Rev. Fluids 3 (1), 011301(R).CrossRefGoogle Scholar
Srinivas, S.S. & Kumaran, V. 2017 Effect of viscoelasticity on the soft-wall transition and turbulence in a microchannel. J. Fluid Mech. 812, 1076–1118.CrossRefGoogle Scholar
Stewartson, K. & Stuart, J.T. 1971 A non-linear instability theory for a wave system in plane Poiseuille flow. J. Fluid Mech. 48 (3), 529–545.CrossRefGoogle Scholar
Stuart, J.T. 1960 On the non-linear mechanics of wave disturbances in stable and unstable parallel flows. Part 1. The basic behaviour in plane Poiseuille flow. J. Fluid Mech. 9 (3), 353–370.CrossRefGoogle Scholar
Su, Z.-G., Zhang, Y.-M., Luo, K. & Yi, H.-L. 2020 Instability of electroconvection in viscoelastic fluids subjected to unipolar injection. Phys. Fluids 32 (10), 104102.CrossRefGoogle Scholar
Sun, G., Wan, D. & Zhang, M. 2021 Nonlinear evolutions of streaky structures in viscoelastic pipe flows. J. Non-Newtonian Fluid Mech. 295, 104622.CrossRefGoogle Scholar
Sureshkumar, R. & Beris, A.N. 1995 a Effect of artificial stress diffusivity on the stability of numerical calculations and the flow dynamics of time-dependent viscoelastic flows. J. Non-Newtonian Fluid Mech. 60 (1), 53–80.CrossRefGoogle Scholar
Sureshkumar, R. & Beris, A.N. 1995 b Linear stability analysis of viscoelastic Poiseuille flow using an Arnoldi-based orthogonalization algorithm. J. Non-Newtonian Fluid Mech. 56 (2), 151–182.CrossRefGoogle Scholar
Terrones, G. & Pearlstein, A.J. 1989 The onset of convection in a multicomponent fluid layer. Phys. Fluids A 1 (5), 845–853.CrossRefGoogle Scholar
Thompson, M.C., Leweke, T. & Provansal, M. 2001 Kinematics and dynamics of sphere wake transition. J. Fluids Struct. 15 (3–4), 575–585.CrossRefGoogle Scholar
Trefethen, L.N. 2000 Spectral Methods in MATLAB. SIAM.CrossRefGoogle Scholar
Virk, P.S. 1975 Drag reduction fundamentals. AIChE J. 21 (4), 625–656.CrossRefGoogle Scholar
Willis, A.P. 2017 The openpipeflow Navier–Stokes solver. SoftwareX 6, 124–127.CrossRefGoogle Scholar
Wu, J., Traoré, P., Pérez, A.T. & Zhang, M. 2015 Numerical analysis of the subcritical feature of electro-thermo-convection in a plane layer of dielectric liquid. Physica D 311, 45–57.CrossRefGoogle Scholar
Zhang, M. 2016 Weakly nonlinear stability analysis of subcritical electrohydrodynamic flow subject to strong unipolar injection. J. Fluid Mech. 792, 328–363.CrossRefGoogle Scholar
Zhang, M. 2021 Energy growth in subcritical viscoelastic pipe flows. J. Non-Newtonian Fluid Mech. 294, 104581.CrossRefGoogle Scholar
Zhang, M., Lashgari, I., Zaki, T.A. & Brandt, L. 2013 Linear stability analysis of channel flow of viscoelastic Oldroyd-B and FENE-P fluids. J. Fluid Mech. 737, 249–279.CrossRefGoogle Scholar
Figure 0

Figure 1.(a) Comparison of the eigenspectra of viscoelastic pipe flow of Oldroyd-B fluids at $Re=800$, $\alpha =1$, $Wi=65$, $\beta =0.65$; the inset is an enlarged view near the unstable mode. (b) Effect of $E$ on the eigenvalue of the least stable centre mode ($c=\omega /\alpha$ is used to facilitate the comparison and $c_i$ is the imaginary part and $c_r$ the real part) for UCM fluids at $Re=6000$ and $\alpha =2$.

Figure 1

Table 1.Convergence of the eigenvalue of the single unstable mode in two different formulations ($up$-$c$ and $\psi$-$g$) of the viscoelastic pipe flow at $Re=800$, $\alpha =1$, $Wi=65$ and $\beta =0.65$.

Figure 2

Table 2.Convergence of the Landau coefficients of the viscoelastic pipe flow of Oldroyd-B fluids at $Wi=65$, $\beta =0.65$, ${Re}_{c}=265.572335$ and $\alpha _{c}=0.515525$ using the $\psi$-$g$ formulation. Note that $a_{1}^* = a_{1} \times {10}^{5}$.

Figure 3

Table 3.Convergence of the Landau coefficients of the viscoelastic pipe flow of Oldroyd-B fluids at $Wi=65$, $\beta =0.65$, ${Re}_{c}=265.572335$ and $\alpha _{c}=0.515525$ using the $up$-$c$ formulation. Note that $a_{1}^* = a_{1} \times {10}^{5}$.

Figure 4

Figure 2.Dispersion relation at $Re_{c}=265.572335$, $\alpha _{c}=0.515525$, $Wi=65$ and $\beta =0.65$: (a) $\mu _r$ and (b) $\mu _i$ as functions of $\alpha$.

Figure 5

Figure 3.(a) Neutral curve in a $(\alpha,Re)$ plane with the colourbar for levels of growth rate and the area inside the loop is linear unstable. (b) Amplitudes of velocity eigenfunctions, where solid and dash–dotted curves are for the left-end critical point $Re_{c}=265.572335$, $\alpha _{c}=0.515525$ and dashed and dotted curves for the right-end critical point $Re_{c}=3553.575$, $\alpha _{c}=4.7337$. (c) Neutral curves in a log-log plot at a fixed viscosity ratio $\beta =0.65$ for various elasticity numbers $E$ and Weissenberg numbers $Wi$, where the discrete data points labelled by blue squares are manually digitized from figure 18(a) in Chaudhary et al. (2021). (d) Neutral curves in a linear-linear plot at various Weissenberg numbers $Wi$ for $\beta =0.65$.

Figure 6

Figure 4.Adjoint and direct eigenfunctions at $Re_{c}=265.572335$, $\alpha _{c}=0.515525$, $Wi=65$ and $\beta =0.65$, corresponding to the left-end critical point in figure 3(a): (a) adjoint velocity components $\tilde {u}_{z,1}^{{\dagger} }$ and $\tilde {u}_{r,1}^{{\dagger} }$; (b) adjoint conformation tensor components $\tilde {g}_{rz,1}^{{\dagger} }$ and $\tilde {c}_{rz,1}^{{\dagger} }$; (c) adjoint conformation tensor components $\tilde {g}_{zz,1}^{{\dagger} }$ and $\tilde {c}_{zz,1}^{{\dagger} }$; (d) velocity components $\tilde {u}_{z,1}$ and $\tilde {u}_{r,1}$; (e) conformation tensor components $\tilde {c}_{rr,1}$ and $\tilde {c}_{rz,1}$; (f) conformation tensor components $\tilde {c}_{\theta \theta,1}$ and $\tilde {c}_{zz,1}$. The conformation tensor components $\tilde {c}_{rr,1}$, $\tilde {c}_{rz,1}$, $\tilde {c}_{\theta \theta,1}$, $\tilde {c}_{zz,1}$ in the $\psi$-$g$ formulation are transformed from $\tilde {g}_{rr,1}$, $\tilde {g}_{rz,1}$, $\tilde {g}_{\theta \theta,1}$, $\tilde {g}_{zz,1}$ according to (2.7).

Figure 7

Figure 5.Eigenspectra obtained from a linear stability analysis of pipe Poiseuille flows of UCM fluids: (a) at $Re=1000$, $\alpha =0.2$ and $Wi=60$; (b) at $Re=260$, $\alpha =0.94$ and $Wi=14.5$.

Figure 8

Figure 6.Critical conditions and the Landau coefficient for pipe flow of Oldroyd-B fluids at various $Wi$ and $\beta$ (small and moderate values): (a) critical Reynolds number $Re_{c}$; (b) critical wavenumber $\alpha _{c}$; (c) linear growth rate $\mu _{r}$; (d) Landau coefficient $a_{3,r}$. All the cases here are supercritical as $a_{3,r}$ is negative.

Figure 9

Figure 7.Critical conditions and the Landau coefficient at various $Wi$ for pipe flow of UCM fluids ($\beta =0$): (a) critical Reynolds number $Re_{c}$; (b) critical wavenumber $\alpha _{c}$; (c) group velocity $c_{g,r}$; (d) Landau coefficient $a_{3,r}$. All the cases here are supercritical as $a_{3,r}$ is negative.

Figure 10

Figure 8.Flow pattern of the second-order solution at $Re_c = 1054.35$, $\alpha _{c}=0.1359$, $Wi=203$ and $\beta =0.56$. Results for (a) $\tilde {u}_{z,21}/U_z$; (b) $\textrm {tr}(\tilde {\boldsymbol {c}}_{21})$ along with the blue curves for streamlines.

Figure 11

Figure 9.Critical conditions and the Landau coefficient for pipe flow of Oldroyd-B fluids at various $Wi$ and $\beta$ (large values): (a) critical Reynolds number $Re_{c}$; (b) critical wavenumber $\alpha _{c}$; (c) Landau coefficient $a_{3,r}$. All the cases here are subcritical as $a_{3,r}$ is positive.

Figure 12

Figure 10.Continuous variation of the Landau coefficient with viscosity ratio $\beta$ at various Weissenberg numbers $Wi$: (a) critical Reynolds number $Re_{c}$; (b) critical wavenumber $\alpha _{c}$; (c) Landau coefficient $a_{3,r}$; (d) Landau coefficient $a_{3,r}$ near the bifurcation boundary. Data points below the horizontal dashed line in (c,d) denote supercritical bifurcations, while those above the line denote subcritical bifurcations.

Figure 13

Figure 11.Critical conditions in a large $(Wi,\beta )$ plane. (a,b) Critical Reynolds number $Re_c$; (c,d) critical wavenumber $\alpha _c$; (e,f) elasticity number at the critical condition $E_c=Wi/Re_c$. The right column is presented because for some values of $Wi,\beta$, there is no linear instability. The dashed red line labels the bifurcation boundary, see the discussion in figure 12.

Figure 14

Figure 12.Landau coefficient in a $(Wi,\beta )$ plane at linear critical conditions for high $Wi$ (a) and in a curve plot for low $Wi$ (b). The dashed red line representing $\beta \approx 0.785$ in (a) marks the bifurcation boundary, the same as those in figure 11; this boundary almost coincides with the contour of $E_c \approx 0.46$ in figure 11(e).

Figure 15

Figure 13.A scaling law of the Landau coefficient $a_3$: (a) collapse of the $a_{3,r}Wi$-$\beta$ curves for different $Wi$; (b) values of $a_{3,r}Wi$ as a function of $Wi$ for some $\beta$. The horizontal dashed line in (a) marks the bifurcation boundary.

Figure 16

Figure 14.Illustration of different bifurcation types ($y$-axis: disturbance amplitude $|A|$; $x$-axis: the control parameter $Re$) near the linear critical points (marked with grey dotted circles). Solid curves are stable solutions and dashed curves unstable solutions; arrows indicate the evolution direction of solutions. (a) Supercritical bifurcation ($a_{3,r}<0$); (b) subcritical bifurcation ($a_{3,r}>0$); (c) a subcritical bifurcation but may appear to be supercritical in experiments (black dots). The red dotted line indicates the systematic disturbance, which cannot be entirely removed. The black dots represent experiment measurements with arbitrary error bars.

Figure 17

Figure 15.Evolution of disturbance kinetic energy (normalized by the kinetic energy of the corresponding Newtonian base flow $E_{k,b}$) at: (a) $Re=280$, $\alpha _{c}=0.515525$, $Wi=65$ and $\beta =0.65$; (b) $Re=90$, $\alpha _{c}=1.550722$, $Wi=65$ and $\beta =0.90$. These results are obtained from DNS and $N_k$ is the number of Fourier modes in the axial direction. Generally, $N_{k}\geqslant 4$ suffices to capture the weak nonlinearity for the evaluation of the Landau coefficients.

Figure 18

Figure 16.Evolution and post-processing of disturbance amplitude $|A|$ at $Re=280$, $\alpha _{c}=0.515525$, $Wi=65$ and $\beta =0.65$ obtained with $N_r=400$ and $N_k=5$ from the DNS: (a) original signal; (b) log plot of (a) demonstrating the linear growth and nonlinear saturation; (c) time derivative of $\log |A|$ vs $|A|^{2}$ with the vertical intercept giving the linear growth rate $c_{1} \approx 3.05 \times 10^{-4}$ and the gradient near this intercept point giving the Landau coefficient $c_{3} \approx -168$.

Figure 19

Figure 17.Continuous variation of $c=\omega /\alpha$ with viscosity ratio $\beta$ at three sets of parameters: (a) imaginary part $c_i$; (b) real part $c_r$. Linear instability exists ($c_i>0$) in the range of $\beta$ above the horizontal black line in (a).

Figure 20

Figure 18.Eigenfunctions corresponding with the unstable mode shown in figure 5(a) at Reynolds number $Re=1000$, wavenumber $\alpha =0.2$ and Weissenberg number $Wi=60$: (a) velocity components $\tilde {u}_{z,1}$ and $\tilde {u}_{r,1}$; (b) conformation tensor components $\tilde {c}_{rr,1}$ and $\tilde {c}_{rz,1}$; (c) conformation tensor components $\tilde {c}_{\theta \theta,1}$ and $\tilde {c}_{zz,1}$.

Figure 21

Figure 19.Neutral curves for pipe flow of UCM fluids: $(a)$ at varying Weissenberg numbers $Wi$; (b) at $Wi=19$ with the colourbar indicating the linear growth rate $\mu _r$. The two arrows in (b) point at the two critical points.