2. Formulation of the problem

2.1. Governing equations

Consider the flow of an incompressible, shear-thinning fluid through an infinitely long circular pipe. We work in cylindrical coordinates $(r,\theta ,z)$ , where the radial, azimuthal and axial components of the velocity vector are denoted as $u$ , $v$ and $w$ , respectively. The velocity $\boldsymbol {u}=[u,v,w](r,\theta ,z,t)$ and the pressure $p(r,\theta ,z,t)$ are assumed to be governed by the non-dimensional Cauchy momentum equation and incompressibility condition

(2.1a) \begin{align} \frac {\partial \boldsymbol {u}}{\partial t} + (\boldsymbol {u} \boldsymbol{\cdot }\boldsymbol{\nabla }) \boldsymbol {u} &= -\boldsymbol{\nabla }\left(p-\frac {q}{{\textit{Re}}} z\right) + \frac {1}{{\textit{Re}}}\boldsymbol{\nabla }\boldsymbol{\cdot }(2\mu D), \end{align}

(2.1b) \begin{align}\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol {u} &= 0, \end{align}

with the strain rate tensor $D=(\boldsymbol{\nabla }\boldsymbol {u}+(\boldsymbol{\nabla }\boldsymbol {u})^{\text{T}})/2$ , where $T$ denotes transpose, and normalised dynamic viscosity $\mu (r,\theta ,z,t)$ . The length scale is the radius of the pipe, $R^*$ , the velocity scale is the centreline velocity of the laminar base flow, $U_c^*$ , and the pressure scale is $\rho ^* U_c^{*2}$ , where $ \rho ^*$ is the density of the fluid. The scaled constant pressure gradient $q\gt 0$ drives the flow. The no-slip conditions $u=v=w=0$ are imposed on the pipe wall $r=1$ .

We adopt the Carreau–Yasuda model (Carreau Reference Carreau1972; Yasuda, Armstrong & Cohen Reference Yasuda, Armstrong and Cohen1981)

(2.2) \begin{equation} \mu =\mu _\infty +(1-\mu _\infty )\{1+(\lambda \dot {\gamma })^a\}^{(n-1)/a}, \end{equation}

where the strain rate magnitude $\dot {\gamma }=(2D\negthinspace :\negthinspace D)^{1/2}$ is also called the shear rate in the GNF community. The Reynolds number is defined by ${\textit{Re}}=\rho ^* R^*U_c^*/\mu ^*_0$ , using the dimensional viscosity at zero shear rate, $\mu ^*_0$ . In (2.2), $\mu _\infty =\mu _\infty ^*/\mu ^*_0$ , the ratio of viscosities at infinite and zero shear rates, typically ranges from $ 10^{-3}$ to $ 10^{-4}$ for shear-thinning fluids in experiments (Escudier et al. Reference Escudier, Poole, Presti, Dales, Nouar, Desaubry, Graham and Pullum2005, Reference Escudier, Nickson and Poole2009). The Carreau number, $\lambda =(U_c^*/R^*)\lambda ^*$ , can be determined from the ‘time constant’ $\lambda ^*$ , which represents an inverse shear rate marking the onset of shear thinning. When the Yasuda parameter $a=2$ , the constitutive relation reduces to that in the Carreau model, and our non-dimensional formulation coincides with that used in Liu & Liu (Reference Liu and Liu2012).

2.2. Base flow

The constant $q$ is determined such that the laminar base flow has a centreline velocity of unity. Substituting $(u,v,w,p)=(0,0,\overline {w}(r),\overline {p}(r))$ into the governing equations, we find that $\overline {w}$ and $q$ can be determined by solving

(2.3a) \begin{align}& \qquad \qquad\quad r^{-1}(r\overline {\mu }\,\overline {w}^{\prime})'=-q, \end{align}

(2.3b) \begin{align}&\overline {\mu }(r)=\mu _\infty +(1-\mu _\infty )\{1+|\lambda \overline {w}^{\prime}|^a\}^{(n-1)/a}, \end{align}

subject to the boundary conditions $\overline {w}(1)=0$ , $\overline {w}(0)=1$ . The primes denote the ordinary differentiation, and overlines indicate laminar base flow quantities. We solve (2.3) for $\overline {w}(r)$ using a numerical method (see Appendix A). Outcomes agree with the analytical solution involving the Gauss hypergeometric function recently found by Wang (Reference Wang2022). For our purposes, numerical computation is more convenient, as it facilitates the easy evaluation of higher-order derivatives of the profile.

Figure 1. (a) Base flow profile for $n=0.5$ and 0.2. The blue symbols are $\overline {w}$ found by the numerical computation of (2.3a ) with $ \mu _\infty = 0$ , $a=2$ and $ \lambda =100$ . The black lines denote the power-law approximation $\bar {w}=1-r^{1+1/n}$ . The dashed line represents the parabolic base flow profile of the Hagen–Poiseuille flow ( $n=1$ ). (b) Complex growth rate at the parameters $ (\mu _\infty ,a,n,\lambda ,{\textit{Re}}) = (0,2,0.8,10,1831.5)$ and the wavenumbers $ k=1$ , $ m=0$ . The black squares are taken from Liu & Liu (Reference Liu and Liu2012) for comparison.

It is common to set $\mu _{\infty }=0$ for simplicity in numerical computations (e.g. Liu & Liu (Reference Liu and Liu2012), Plaut et al. (Reference Plaut, Roland and Nouar2017)), and we will also examine this idealised case in § 3. In this case, the large- $\lambda$ limit of the base flow and the viscosity can be approximated by the power law as

(2.4) \begin{eqnarray} \overline {w}=(1-r^{1+1/n})+\cdots , \qquad \overline {\mu }\approx (\lambda |\overline {w}^{\prime}| )^{n-1}+\cdots ,\end{eqnarray}

except for a small region around the centreline of the pipe where $r=O(\lambda ^{-n})$ (see Appendix A). Figure 1(a) compares the numerical solution of (2.3) for $\lambda =100$ with the power-law approximation. At this value of $\lambda$ , the Carreau fluid behaves almost like a power-law fluid.

2.3. Parameters used in experiments

The five flow parameters $(\mu _\infty ,a,n,\lambda )$ and ${\textit{Re}}$ defined above are useful for theoretical analysis but are not optimal for organising experimental data. We first note that to use the Carreau–Yasuda law, the constants $\mu _0^*$ , $\mu _\infty ^*$ , $\lambda ^*$ , $n$ and $a$ need to be found by fitting experimental data for the specific fluid in question. Since $\lambda ^*$ is a constant particular to the fluid, when the Reynolds number varies in the experiments, $\lambda$ is not a constant but rather a quantity proportional to ${\textit{Re}}$ . Therefore it is more convenient to specify $\varLambda =\lambda /{\textit{Re}}=(\lambda ^*\mu _0^*)/(\rho ^* R^{*2})$ instead of $\lambda$ . Note that $\varLambda$ depends on the pipe radius $R^*$ .

Another important consideration is the definition of the Reynolds number. In experiments, it is more convenient to fix the flow rate rather than the pressure gradient, and thus the bulk velocity $U^*_b$ is used as the velocity scale. For example, Escudier et al. (Reference Escudier, Poole, Presti, Dales, Nouar, Desaubry, Graham and Pullum2005) and Wen et al. (Reference Wen, Poole, Willis and Dennis2017) employed the Reynolds number

(2.5) \begin{equation} {\textit{Re}}_b=\frac {\rho ^* (2R^*) U_b^*}{\langle \mu_\textit{wall}^{*}\rangle} \end{equation}

using the dimensional viscosity at the wall, $\mu _{{wall}}^*$ . Angle brackets denote the average over $\theta$ , $z$ and $t$ . Recalling the non-dimensionalisation introduced in § 2.1, we have $\langle \mu _{{wall}}^* \rangle =\mu _0^*\langle \mu \rangle |_{r=1}$ and $U_b^*=(\int ^1_0U_c^*\langle w \rangle r {\rm d}r)/(\int ^1_0 r {\rm d}r)$ . Therefore, the conversion recipe for the two Reynolds numbers can be obtained as

(2.6) \begin{equation} {\textit{Re}}_b=\frac {4{\textit{Re}}}{\langle \mu \rangle |_{r=1}}\int ^1_0 \langle w \rangle r{\rm d}r. \end{equation}

For the base flow, the ratio $(\textit{Re}_b/\textit{Re})=({4}/{\overline{\mu}|_{r=1}})\int ^1_0 \overline {w}r{\rm d}r$ can be easily computed numerically. The power-law fluid approximation (2.4) suggests that for large $\lambda$

(2.7) \begin{equation} \frac {{\textit{Re}}_b}{{\textit{Re}}}\approx \frac {2n\big(1+\frac {1}{n}\big)^{2-n} }{(3n+1)}\lambda ^{1-n}. \end{equation}

2.4. Linear stability analysis

In the linear stability analysis, we assume that the perturbation $[\tilde {u},\tilde {v},\tilde {w},\tilde {p}]=[u,v,w-\overline {w},p-\overline {p}]$ can be written in the form $[\hat {u}(r),\hat {v}(r),\hat {w}(r),\hat {p}(r)]\exp (im\theta +ik(z-ct))+\text{c.c.}$ , where c.c. is the complex conjugate. For given flow parameters, the azimuthal wavenumber $m$ and the axial wavenumber $k$ , we can numerically solve the linearised governing equations for the complex phase speed $c=c_r+ic_i$ . The resultant stability problem is identical to those described in Liu & Liu (Reference Liu and Liu2012), with a correction to the obvious typo in their (23)

(2.8a) \begin{align} &\qquad\qquad\qquad\qquad\qquad\quad \hat {u}' + \frac {\hat {u}}{r} + \frac {im}{r} \hat {v} + ik \hat {w} = 0, \end{align}

(2.8b) \begin{align}&\quad ik(\bar {w}-c)\hat {u} +\hat {p}' = \frac {1}{{\textit{Re}}} \left \{ \bar {\mu } \left [ \mathcal{L} \hat {u} - \frac {\hat {u}}{r^2} - \frac {2im}{r^2} \hat {v} \right ] + 2 \bar {\mu }' \hat {u}' + \check {\mu } \big(ik\hat {w}^{\prime} -k^2 \hat {u}\big)\right \}\!, \end{align}

(2.8c) \begin{align}&\quad ik(\bar {w}-c)\hat {v} +\frac {im}{r} \hat {p}= \frac {1}{{\textit{Re}}} \left \{ \bar {\mu } \left [ \mathcal{L} \hat {v} - \frac {\hat {v}}{r^2} + \frac {2im}{r^2} \hat {u}\right ] + \bar {\mu }' \left ( \frac {im}{r} \hat {u} + \hat {v}' - \frac {\hat {v}}{r} \right )\! \right \}\!, \end{align}

(2.8d) \begin{align}& ik(\bar {w}-c)\hat {w} + \bar {w}^{\prime} \hat {u}+ik \hat {p} = \frac {1}{{\textit{Re}}} \left \{ \bar {\mu } \mathcal{L} \hat {w} + \left(\bar {\mu }'+\check {\mu }'+\frac {\check {\mu }}{r}\right) \!\left ( \hat {w}^{\prime} + ik\hat {u} \right ) + \check {\mu } \!\left ( \hat {w}^{\prime\prime} + ik\hat {u}' \right )\! \right \}\!. \end{align}

Here, we used the shorthand notations

(2.9) \begin{align}& \check {\mu }(r)=(n-1)(1-\mu _\infty )\{1+|\lambda \bar {w}^{\prime}|^a\}^{(n-1-a)/a}|\lambda \bar {w}^{\prime}|^a, \end{align}

(2.10) \begin{align}&\qquad\qquad\quad \mathcal{L} = \frac {\partial ^2}{\partial r^2} + \frac {1}{r} \frac {\partial }{\partial r} - \frac {m^2}{r^2} - k^2. \end{align}

The boundary conditions become $\hat {u}(1)=\hat {v}(1)=\hat {w}(1)=0$ .

Our numerical code is based on the method described in Deguchi & Nagata (Reference Deguchi and Nagata2011), where the poloidal–toroidal potential approach is employed. Spatial discretisation is performed using Fourier–Galerkin and Chebyshev-collocation methods. The radial basis functions follow those used in Deguchi & Walton (Reference Deguchi and Walton2013). This code has been repeatedly validated by successfully reproducing known travelling waves in Newtonian fluids (see Song & Deguchi (Reference Song and Deguchi2025), for example).

Here, we present only the details of the linear stability analysis; for the nonlinear computations, refer to the papers cited above. In our flow domain, infinitesimally small velocity perturbations are written as $[\tilde {u},\tilde {v},\tilde {w},\tilde {p}]=\boldsymbol{\nabla }\times \boldsymbol{\nabla }\times (\phi \textbf{e}_r)+\boldsymbol{\nabla }\times (\psi \textbf{e}_r)$ , where $\phi$ and $\psi$ are poloidal and toroidal potentials, respectively. The two equations for them can be obtained by operating $\textbf{e}_r\boldsymbol{\cdot }\boldsymbol{\nabla }\times \boldsymbol{\nabla }\times$ and $\textbf{e}_r\boldsymbol{\cdot }\boldsymbol{\nabla }\times$ on the momentum equations. The potentials are then approximated using a finite set of basis functions as follows:

(2.11a) \begin{align} & \phi =\sum _{l=0}^N a_l \varPhi _l^{(m)}(r)\exp (im\theta +ik(z-ct))+\text{c.c.}, \end{align}

(2.11b) \begin{align} & \psi =\sum _{l=0}^N b_l \varPsi _l^{(m)}(r)\exp (im\theta +ik(z-ct))+\text{c.c.}, \end{align}

where c.c. denotes complex conjugate and

(2.12a) \begin{align} \varPhi _l^{(m)}(r)&= \left \{\!\! \begin{array}{c} r(1-r^2)^2T_{2l}(r),\qquad \text{if $m=0$},\qquad \qquad \qquad \qquad \\ r(1-r^2)^2T_{2l+1}(r),\qquad \text{if $m$ is odd}, \qquad \qquad \qquad \\ r^3(1-r^2)^2T_{2l}(r),\qquad \text{if $m$ is even and $m\neq 0$}, \end{array} \right . \end{align}

(2.12b) \begin{align} \varPsi _l^{(m)}(r)&= \left \{\!\! \begin{array}{c} r(1-r^2)T_{2l}(r),\qquad \text{if $m$ is even},\qquad \\ r(1-r^2)T_{2l+1}(r),\qquad \text{if $m$ is odd}. \end{array} \right . \\[12pt] \nonumber \end{align}

Here, $T_l(r)$ is the $l$ th Chebyshev polynomial. Evaluating the governing equations at the collocation points $r_k=\cos ((k+1)\pi /(2N+4)), k=0,1,\ldots ,N$ , we find that non-trivial coefficients $a_l$ and $b_l$ exist only when the growth rate $\sigma =-ikc$ is an eigenvalue of the resulting algebraic eigenvalue problem. This problem can be solved using LAPACK solvers or Rayleigh quotient iteration scheme. For Newtonian fluids, we verified that the computed eigenvalues match those listed in Schmid & Henningson (Reference Schmid and Henningson1994) and Meseguer & Trefethen (Reference Meseguer and Trefethen2003) to at least nine decimal places. Shear-thinning effects on the eigenvalues were assessed through comparison with Liu & Liu (Reference Liu and Liu2012); see figure 1(b). For most parameters in this paper, the code using $N=200$ achieves excellent convergence without any spurious eigenvalues.

3. Linear stability results for $\mu _{\infty }=0$

To obtain a general understanding of the stability characteristics of shear-thinning fluids, we first focus on the case where $a=2$ and $\mu _{\infty }=0$ . All the stability results presented in this paper use $m=1$ unless otherwise stated.

Figure 2(a) shows the stability diagram in the ${\textit{Re}}$ – $k$ plane for $n=0.2$ . At this value of $n$ , instability does indeed occur for $m=1$ . The dashed and solid lines are the neutral curves for $\lambda =5$ and 100, respectively. The region enclosed by the curve is unstable. For example, at $\lambda =100$ , an unstable mode is obtained at a point marked by the symbol. The resolution test for this mode is presented in table 1, confirming that the typically chosen value of $N=200$ is more than sufficient. Figure 3 illustrates the corresponding eigenfunction, which is also well converged. For other values of $m$ , no instability is observed.

Figure 2. Stability results for $n = 0.2$ and $\mu _{\infty }=0$ . The Carreau model ((2.2) with $a=2$ ) is used except for the magenta open circles. All the instabilities presented in the figures correspond to $m=1$ . (a) Neutral curves in the $ {\textit{Re}}-k$ plane at $\lambda =100$ (solid) and $5$ (dashed). The red lines are the long-wavelength asymptotic results for $\lambda =100$ . The star symbol represents the parameter at which the unstable eigenvalue in table 1 is computed. (b) Neutral curves in the $ \lambda ^{1-n}{\textit{Re}}-k$ plane at $\lambda =50$ (black filled circles) and 100 (solid line). The magenta open circles indicate the neutral curve for the Cross model (3.1) with $\lambda = 100$ .

Figure 3. Eigenfunction of the unstable mode found at the symbol in figure 2(a). See table 1 for the parameters used. The solid and dashed lines are the real and imaginary parts, respectively.

The unstable mode occurs universally in flows that can be approximated by power-law fluids. Under the approximation introduced in § 2.2, the viscosity functions behave like $\overline {\mu }\approx (\lambda |\overline {w}^{\prime}| )^{n-1}$ and $\check {\mu }\approx (n-1)(\lambda |\overline {w}^{\prime}| )^{n-1}$ . Therefore, the effective Reynolds number in the stability problem (2.8) is $\lambda ^{1-n}{\textit{Re}}$ . In fact, when $\lambda$ is sufficiently large, the neutral curves can be better organised using the rescaled Reynolds number $\lambda ^{1-n}{\textit{Re}}$ (equivalently ${\textit{Re}}_b$ , see (2.7)), as shown in figure 2(b). Similar results can be obtained for the Cross model:

(3.1) \begin{equation} \mu =\mu _\infty +\frac {1-\mu _\infty }{1+(\lambda \dot {\gamma })^{1-n}}. \end{equation}

The open magenta circles in figure 2(b) represent the neutral points for this model with $\mu _{\infty }=0,\lambda =100$ . Apart from the viscosity model, the computational method remains the same as in the Carreau fluid cases. The viscosity defined by (3.1) also behaves like a power-law fluid when $\lambda$ is sufficiently large. Consequently, the stability results align with the same universal curve as in the Carreau fluid case when the rescaled Reynolds number $\lambda ^{1-n}{\textit{Re}}$ is used.

Table 1. The most unstable complex growth rate $\sigma =-ikc$ found at the point marked in figure 2(a). The parameters are $(a,\mu _{\infty },n,\lambda ,{\textit{Re}},k)=(2,0,0.2,100,2800,0.4)$ .

The origin of the instability cannot be an inviscid mechanism, as the sufficient condition for stability established by Batchelor & Gill (Reference Batchelor and Gill1962) is fulfilled; see Appendix B. Therefore, the unstable mode is of the viscous type. Akin to the Tollmien–Schlichting (TS) wave, in figure 3, the critical layer $r=0.969$ , computed from the phase speed $c=0.171$ , is located near the wall. It appears that the thin critical layer lies inside a thicker one. This second layer is consistent in thickness with the near-wall boundary layer developed by the base flow when $n$ is small (see figure 1).

In figure 2, at large Reynolds numbers, the wavenumber $k$ along the neutral curves is inversely proportional to ${\textit{Re}}$ , suggesting that the limiting case can be described by asymptotic analysis. The derivation of the leading-order problem, hereafter referred to as the long-wavelength limit problem, is straightforward. It can be obtained by substituting the regular expansion $[\hat {u},\hat {v},\hat {w},\hat {p}]=[{\textit{Re}}^{-1}\hat {u}_0,{\textit{Re}}^{-1}\hat {v}_0,\hat {w}_0,{\textit{Re}}^{-2}\hat {p}_0]+\cdots$ into the linearised Navier–Stokes equations, keeping $k_0={\textit{Re}} \,k$ and $c$ as $O(1)$ , and discarding small terms

(3.2a) \begin{align} &\qquad\qquad\qquad\qquad\quad \hat {u}_0' + \frac {\hat {u}_0}{r} + \frac {im}{r} \hat {v}_0 + ik_0 \hat {w}_0 = 0, \end{align}

(3.2b) \begin{align} &\qquad \begin{aligned} ik_0(\bar {w}-c)\hat {u}_0 +\hat {p}_0' &= \bar {\mu } \left [ \mathcal{L}_0 \hat {u}_0 - \frac {\hat {u}_0}{r^2} - \frac {2im}{r^2} \hat {v}_0 \right ]+ 2 \bar {\mu }' \hat {u}_0' + \check {\mu } ik_0\hat {w}_0' , \end{aligned} \end{align}

(3.2c) \begin{align} &\begin{aligned} ik_0(\bar {w}-c)\hat {v}_0 +\frac {im}{r} \hat {p}_0&= \bar {\mu } \left [ \mathcal{L}_0 \hat {v}_0 - \frac {\hat {v}_0}{r^2} + \frac {2im}{r^2} \hat {u}_0\right ] + \bar {\mu }' \left ( \frac {im}{r} \hat {u}_0 + \hat {v}_0' - \frac {\hat {v}_0}{r} \right ), \end{aligned} \end{align}

(3.2d) \begin{align} &\qquad\quad\begin{aligned} ik_0(\bar {w}-c)\hat {w}_0 + \bar {w}^{\prime} \hat {u}_0 &= \bar {\mu } \mathcal{L}_0 \hat {w}_0 + \left(\bar {\mu }'+\check {\mu }'+\frac {\check {\mu }}{r}\right) \hat {w}_0' + \check {\mu } \hat {w}_0^{\prime\prime} , \end{aligned} \end{align}

where $\mathcal{L}_0= \partial ^2/\partial r^2 + r^{-1}\partial /\partial r - m^2/r^2$ . The boundary conditions are $\hat {u}_0(1)=\hat {v}_0(1)=\hat {w}_0(1)=0$ , so the basis functions (2.12) can be employed. The above equations resemble a linearised version of Prandtl’s boundary layer equations, but they apply across the entire flow region. Using the parameters from figure 2(a) with $\lambda =100$ , we can solve (3.2) numerically and found that there are two rescaled wavenumbers $k_0=449.6$ and $9635.34$ that make the real part of $\sigma _0$ zero for the leading mode. Those values determine the red lines seen in figure 2(a), which gives a good approximation of the neutral curve when ${\textit{Re}}$ is large. As this result suggests, the long-wavelength problem serves as a useful tool for investigating the existence of instability.

Of particular interest for practical applications is determining the values of $n$ for which instability occurs. Figure 4(a) shows the neutral stability curves in the $n$ – ${\textit{Re}}$ plane for $\lambda =100$ and various values of $k$ . The envelope of the curves shown in the figure gives the stability boundary. This boundary seems to exhibit a well-defined cutoff value of $n$ at large Reynolds numbers. Approaching the cutoff, the optimum values of $k$ that define the stability boundary decrease. This behaviour of $k$ is typical when a cutoff in shear-flow instability occurs, as first observed in plane Couette–Poiseuille flow by Cowley & Smith (Reference Cowley and Smith1985). Another example can be found in the study of annular Poiseuille flow by Heaton (Reference Heaton2008), where a long-wavelength limit system similar to (3.2) was derived.

Figure 4. The dependence of the stability results on the power-law index $n$ . (a) Neutral curves in the ${\textit{Re}}$ – $n$ plane for $ \lambda =100$ . The dashed line indicates the cutoff value of $n$ for the unstable region. (b) Neutral curves in the $k_0$ – $n$ plane at $ \lambda =100$ , where $k_0=k{\textit{Re}}$ . The black dashed curves are the results for ${\textit{Re}}=10^4, 5\times 10^4$ and $10^5$ . The red solid curve is the long-wavelength asymptotic result. The circle indicates the point used in figure 5(a).

The three dashed curves in figure 4(b) are the neutral curves for ${\textit{Re}}=10^4, 5\times 10^4$ and $10^5$ . Here, the vertical axis is the rescaled wavenumber $k_0={\textit{Re}}\,k$ . As ${\textit{Re}}$ increases, the neutral curves asymptote to the red curve, which is computed using the reduced (3.2). The calculation of the turning point of this curve, indicated by the black circle, provides a convenient way for estimating the cutoff value of $n$ for instability. The value $n\approx 0.35$ found at the black circle, indicated by the vertical dashed line in figure 4(a), indeed represents the critical threshold beyond which the instability no longer occurs at $\lambda =100$ .

Figure 5. The stability analysis based on the long-wavelength limit problem (3.2). (a) Neutral curves in the $n$ – $\lambda$ plane with optimised $k_0$ . (b–d) Streamwise velocity of the neutral eigenfunction in a pipe cross-section, scaled by its local maximum (marked by a star). The parameters used are (b) $(m,\lambda ,n,k_0)= (1,1000, 0.3518, 358)$ ; (c) $(m,\lambda ,n,k_0)= (2,1000, 0.1407, 163)$ ; and (d) $(m,\lambda ,n,k_0)= (3,1000, 0.08816, 167)$ .

The long-wavelength limit system (3.2) allows us to estimate the thresholds for the existence of the $m=1$ instability for each $\lambda$ , as indicated by the black circles in figure 5(a). The threshold value $n=0.35$ is robust with respect to increases in $\lambda$ . This threshold recedes toward smaller values of $n$ as $m$ increases. Panels (b), (c) and (d) present a comparison of the neutral modes corresponding to $m=1$ , $m=2$ and $m=3$ at $\lambda =1000$ .

4. Linear stability results based on experimentally measured parameters

Let us now examine physically meaningful parameters. We first use the long-wavelength limit to narrow down the parameters for which instability may arise. The influence of $\mu _{\infty }$ and $a$ , which were held fixed in the previous section, will also be clarified. We then proceed by lifting the long-wavelength assumption to pinpoint the critical Reynolds number.

The open circles in figure 6(a) present stability results derived from a long-wavelength asymptotic analysis, analogous to those in figure 5(a) but for $(\mu _\infty ,a)=(1.116\times 10^{-4},2.0)$ . These values are taken from the fitting parameters for 7 % aluminium soap (AS); see the first row of table 2. This fluid has unusually small $\mu _{\infty }$ , and its neutral curve closely resembles that for $\mu _{\infty }=0$ (thick grey curve). However, noticeable changes emerge as $\mu _{\infty }$ increases. The filled circles show the stability results for $(\mu _\infty ,a)=(1.207\times 10^{-3},2.01)$ , corresponding to $0.2\%$ polyacrylamide (PAA) used in Escudier et al. (Reference Escudier, Poole, Presti, Dales, Nouar, Desaubry, Graham and Pullum2005); see the second row of table 2.

Figure 6. Similar plots to figure 5(a), using different parameters $ (\mu _\infty ,a)$ . (a) The open and filled circles are results for $ (\mu _\infty ,a) = (1.116\times 10^{-4} , 2.0)$ and $ (\mu _\infty ,a) = (1.207\times 10^{-3}, 2.01)$ , respectively. For comparison, the thick curves with light colour represent the results with $(\mu _\infty ,a) =(0,2)$ (same as figure 5 a). The red dash-dot-dot and dash-dot lines are the value of $n$ for 7 % AS and 0.2 % PAA (see table 2). (b) The solid curves with symbols are computed using $ (\mu _\infty ,a) = (2.1875\times 10^{-2},2.01)$ . For the two dashed curves with a lighter colour, $a$ is set to 1.3 and 1.0, respectively. The red dash-dot line indicates the value of $n$ for blood (see table 2).

Table 2. Carreau–Yasuda model parameters for various fluids. The first row corresponds to the 7 % AS in decalin and m-cresol reported in Myers (Reference Myers2005); the second row to the aqueous solutions of $0.2\,\%$ PAA listed in table 1 of Escudier et al. (Reference Escudier, Poole, Presti, Dales, Nouar, Desaubry, Graham and Pullum2005); and the third row to the blood parameters taken from Boyd et al. (Reference Boyd, Buick and Green2007). The values of $\rho ^*$ are estimated from the densities of the solvent and solute.

By increasing the parameter $\mu _{\infty }$ of the latter neutral curve to $2.1875\times 10^{-2}$ , we obtain the black solid curve in figure 6(b). It is evident that, as $\mu _{\infty }$ increases – even while remaining significantly less than unity – stabilisation occurs in the large $\lambda$ parameter region. Repeating the discussion in § 2.2 while keeping $\mu _{\infty }$ clarifies that $\mu _{\infty }=O(\lambda ^{n-1})$ is large enough to influence the behaviour of viscosity (2.3b ). In fact, for $\mu _{\infty }\gg O(\lambda ^{n-1})$ , the viscosity is nearly Newtonian, leading to the absence of instability. The dashed curves illustrate that a decrease in the value of $a$ also contributes to stabilisation. While the instability exists at $a=1$ , they disappear when $a$ is reduced to 0.64. This value, along with $\mu _{\infty } = 2.1875 \times 10^{-2}$ and $n=0.2128$ , corresponds to the Carreau–Yasuda parameters for blood as reported in Boyd et al. (Reference Boyd, Buick and Green2007).

It should be noted that, while the long-wavelength approximation employed in figure 6 is useful for assessing whether instability occurs for certain choices of $k,{\textit{Re}}$ and $\lambda$ , identifying the physically relevant critical Reynolds number requires solving the full stability problem. Figure 7(a) presents the numerical results for the 0.2 % PAA parameters. Here, following the remark in § 2.3, we fix $\varLambda$ and vary ${\textit{Re}}_b$ . The parameter $\varLambda$ is determined not only by the fluid properties but also by the pipe radius $R^*$ , and it has a significant influence on the flow stability. To observe instability, $\varLambda$ needs to be as small as $5.5\times 10^{-4}$ , which corresponds to a pipe radius of approximately $R^*=8$ m. Escudier et al. (Reference Escudier, Poole, Presti, Dales, Nouar, Desaubry, Graham and Pullum2005) employed a much smaller radius $R^*=$ 0.05 m in their experimental apparatus. The use of this parameter yields $\varLambda = 13.05$ , which unfortunately does not exhibit instability. Theoretically, this result is not surprising, given that figure 6(a) indicates the potential for instability only near $\lambda = O(10)$ . The predicted critical Reynolds number of $O(1)$ , obtained from the relation $\varLambda = \lambda /{\textit{Re}}$ , is too low. In contrast, using 7 % AS allows instability to occur with the experimentally feasible pipe radius of $R^* = 0.05$ m ( $\varLambda =55.17$ ), as shown by figure 7(b). In all cases, the instability forms an isolated region in the $k-{\textit{Re}}_b$ plane. Note that increasing ${\textit{Re}}_b$ leads to flow stabilisation, since the value of $\lambda$ eventually exceeds the Newtonian recovery threshold, $\mu _{\infty }^{1/(n-1)}$ .

Figure 7. Neutral curves in the $k-{\textit{Re}}_b$ plane. Parameter values are listed in table 2. (a) 0.2 % PAA. The values of $\varLambda$ used are $5.5\times 10^{-4}$ , $5\times 10^{-4}$ and $4\times 10^{-4}$ ; (b) $7\,\%$ AS with $\varLambda =55.17$ . The red bullet indicates the critical point $(k_c,{\textit{Re}}_{b,c})\approx (0.366,1.621\times 10^5)$ .

We also carried out stability analyses using the parameters for blood. The most practically relevant values of $\varLambda$ are $8.13$ and $345.04$ , corresponding to the aorta ( $2R^*=2.54\times 10^{-2}$ m) and brachial artery ( $2R^*=3.90\times 10^{-3}$ m), respectively, as listed in table I of Boyd et al. (Reference Boyd, Buick and Green2007); however, no instability was detected for either case. Additional tests with even smaller $\varLambda$ values yielded the same result. This outcome is consistent with the long-wavelength limit analysis shown in figure 6(b), where no sign of instability was observed for $n=0.2128, a=0.64$ .

5. Bifurcation analysis

Bifurcation theory suggests that finite-amplitude travelling-wave solutions emerge from the linear critical points identified so far. This section targets the critical point indicated by the red bullet in figure 7(b).

A naive approach to obtaining nonlinear solutions is to integrate the governing equations in time near the critical parameter, hoping for convergence to a finite-amplitude equilibrium state. Thus, we first attempted direct numerical simulation near the linear critical parameters using the spectral element solver, Semtex (Blackburn et al. Reference Blackburn, Lee, Albrecht and Singh2019, Reference Blackburn, Rudman and Singh2025). However, for $n\lt 0.5$ , numerical instability arose unless a very small time step was chosen. Similar computational challenges for small $n$ have also been reported in Plaut et al. (Reference Plaut, Roland and Nouar2017) and Wen et al. (Reference Wen, Poole, Willis and Dennis2017). Furthermore, the wavelength of the perturbation that yields the critical value ${\textit{Re}}_{b,c}$ in figure 7(b) is fairly long, while high resolution is required in the radial direction. Therefore, we conclude that obtaining a travelling-wave state through direct numerical simulation is not presently feasible.

Our computation here instead utilises Newton’s method implemented in the code by Deguchi & Nagata (Reference Deguchi and Nagata2011). This code employs an analytically derived Jacobian matrix, which eliminates the need for time integration and thus avoids the aforementioned numerical instability. However, deriving the analytic Jacobian matrix is a cumbersome task as $\mu$ depends on the perturbation. This motivates us to consider the Taylor expansion of the viscosity

(5.1) \begin{eqnarray} \mu =\overline {\mu }+\frac {\check {\mu }}{\overline {w}^{\prime}}\frac {\partial \tilde {w}}{\partial r}+\cdots , \end{eqnarray}

which is valid when the perturbation $\tilde {\boldsymbol {u}}=(u,v,w-\overline {w})$ is smaller than the base flow $\overline {w}$ . In our numerical computations, we adopt an approximation where the expansion is truncated at the second term. This approximation is justified for bifurcating solutions near the linear critical point and ensures accurate reproduction of weakly nonlinear analysis results, such as the Landau coefficient. Retaining Fourier modes up to the fourth harmonic is sufficient, since we only investigate the vicinity of the bifurcation point.

The amplitude equations derived from the weakly nonlinear analysis suggest the existence of two types of solutions. This is confirmed by our computations, as shown in figure 8(a). We start Newton’s method by using a single neutral eigenfunction as an initial guess. With an appropriate choice of amplitude, the iterations converge, resulting in the black curve in figure 8(a). In this bifurcation diagram, we measured the solutions by the normalised energy norm of the velocity perturbation

(5.2) \begin{equation} \delta =\left [\int ^1_0\langle \tilde {u}^2+\tilde {v}^2+\tilde {w}^2\rangle r{\rm d}r\,\Big /\int ^1_0r{\rm d}r\right ]^{1/2}. \end{equation}

Figure 8. Bifurcation analysis from the neutral point in figure 7(b). The case of $7\,\%$ AS with $\varLambda =55.17$ . The axial wavenumber is fixed at $k=0.366$ . (a) Bifurcation diagram based on the energy norm of the velocity perturbation, $\delta$ . (b) Same results, shown the ratio of the shear stress from fluctuations to that of the mean flow. (c) Same results, shown in terms of the friction coefficient.

The angle bracket denotes $\theta$ – $z$ average. The structure observed in the isosurface of $\tilde {u}$ , shown in figure 9(a), displays characteristics that prompt us to refer to this solution as the ‘spiral solution’. The helical invariance of the solution clearly arises from the linear neutral eigenfunction with $m=1$ . This invariance is also evident in the streamwise vorticity, $\tilde {\omega }=r^{-1}\partial _{\theta }\tilde {w}-\partial _z \tilde {v}$ (figure 9 b).

Figure 9. Perturbation flow field of the solutions at ${\textit{Re}}_b=1.52\times 10^5$ . The phase is defined by $\varphi =k(z-ct)$ , where $c$ is the phase speed of the travelling wave. Panel (a) shows the streamwise velocity $\tilde {u}$ of the spiral solution. The yellow/blue surfaces depict the positive/negative isosurfaces at 88 % of the maximum magnitude. The red/green surfaces represent the positive/negative isosurfaces at 10 % of the maximum magnitude for the viscosity variation $\mu -\bar {\mu }$ . Panel (b) is the vorticity $\tilde {\omega }$ of the spiral solution, with the isosurfaces plotted at 88 % of the maximum magnitude. Panels (c) and (d) show plots similar to panels (a) and (b), but for the mirror-symmetric solution.

The symmetry of the system indicates that, when a helical neutral mode with a specific pitch exists, there is always another helical mode with the opposite pitch. Using a superposition of the symmetric pair of neutral modes as the initial condition for Newton’s method leads to convergence to a ‘mirror-symmetric solution’ (the blue curve in figure 8 a). The symmetry of the solution is evident from the isosurfaces shown in figure 9(c,d).

As seen in figure 8(a), for both solution branches, the bifurcation is subcritical. Given that $\delta \ll 1$ , the approximation based on the Taylor expansion (5.1) is considered reasonable. The circles show a square-root fit, indicating the weakly nonlinear regime. The weakly nonlinear approximation for the mirror-symmetric solution is only valid very close to the bifurcation point, as shown in the inset. Note that the control parameter in our numerical code is ${\textit{Re}}$ . Thus to draw figure 8(a), we first compute the values of ${\textit{Re}}_b({\textit{Re}})$ and $\delta ({\textit{Re}})$ for each solution, and then plot the parametric curve in the ${\textit{Re}}_b$ – $\delta$ plane.

The numerically obtained solutions exhibit features that are not observed in the travelling waves of Newtonian fluids. In Newtonian pipe flow, the pressure gradient $q/{\textit{Re}}$ can be readily determined from the mean flow alone. However, for GNF, it is also necessary to account for the fluctuating components. This can be found by taking the spatial average of the $z$ -component of the momentum equation (2.1a )

(5.3) \begin{eqnarray} \frac {1}{2}\frac {q}{{\textit{Re}}}=-\frac {1}{{\textit{Re}}}\langle \mu \partial _r w \rangle |_{r=1}=\tau _m+\tau _f. \end{eqnarray}

On the right-hand side of this global momentum balance, $\tau _m=-({1}/{{\textit{Re}}})\langle \mu \rangle \langle w \rangle ' |_{r=1}$ represents the shear stress from the mean flow, while $\tau _f=-({1}/{{\textit{Re}}})\langle (\mu -\langle \mu \rangle )\partial _r(w- \langle w \rangle )\rangle |_{r=1}$ corresponds to the contribution from fluctuations, which is absent in the case of Newtonian fluids. As seen in figure 8(b), the signs of $\tau _m$ and $\tau _f$ are opposite in our solution.

Figure 8(c) shows the variation of

(5.4) \begin{eqnarray} \Delta C_{\!f}=\frac {C_{\!f}-\overline {C}_{\!f}}{\overline {C}_{\!f}}. \end{eqnarray}

Here, $C_{\!f}$ is the friction factor, which is commonly defined by

(5.5) \begin{eqnarray} C_{\!f}=\frac {4R^*\Delta p^*}{(U_b^*)^2}=\frac {4q}{{\textit{Re}} U_b^2}, \end{eqnarray}

using the dimensional pressure gradient driving the flow, $\Delta p^*$ , and the non-dimensional bulk velocity, $U_b=U_b^*/U_c^*=2\int ^1_0\langle w\rangle r{\rm d}r$ . Following the above remark, the definition uses the pressure gradient instead of wall shear. In (5.4), $\overline {C}_{\!f}$ denotes the friction factor of the laminar flow with the same ${\textit{Re}}_b$ as the travelling-wave solution. Therefore, figure 8(c) suggests that, surprisingly, the resistance experienced by the pipe decreases from the laminar state.

Such ‘sublaminar drag’ is widely known to be impossible in Newtonian fluids (see Marusic et al. (Reference Marusic, Joseph and Mahesh2008) and Fukagata, Sugiyama & Kasagi (Reference Fukagata, Sugiyama and Kasagi2009), for example). An analysis of the energy equation in Appendix C reveals that the following inequality holds:

(5.6) \begin{eqnarray} 0\leqslant \frac {8}{U_b^2{\textit{Re}}_b \overline {C}_{\!f} \langle \mu \rangle |_{r=1}}\int ^1_0\langle (\overline {\mu }-\mu ) D\negthinspace :\negthinspace D\rangle r{\rm d}r +\left . \left (\frac {\langle \mu \rangle -\overline {\mu } }{\langle \mu \rangle }\right )\right |_{r=1}+\Delta C_{\!f}. \end{eqnarray}

The second term on the right-hand side is simply a consequence of using the wall viscosity in the definition of ${\textit{Re}}_b$ (see (2.5)), and is therefore of little physical significance. The first term in (5.4) is central to the possibility of sublaminar drag reduction in GNF. The importance of variation of viscosity, $\mu -\overline {\mu }$ , can also be found in the generalised Fukagata–Iwamoto–Kasagi (FIK) identity shown in Appendix D. Compared with the Newtonian version by Fukagata, Iwamoto & Kasagi (Reference Fukagata, Iwamoto and Kasagi2002), an additional term depending on the viscosity variation appears (see (D5)).

Figure 10. The curves indicate the value of $\langle \mu -\bar {\mu }\rangle$ as a function of $r$ for (a) spiral solution; (b) mirror-symmetric solution. Here, ${\textit{Re}}_b=1.52\times 10^5$ . The insets denote colour maps for $\mu -\bar {\mu }$ at $\varphi =0$ .

The viscosity variation $\mu -\overline {\mu }$ is concentrated near the pipe centre; see figure 9(a,c). This may seem counterintuitive at first. However, this behaviour is consistent with (5.1), where the second term is proportional to $\check {\mu }$ . Recall that for large $\lambda$ , the approximation $\check {\mu }\approx (n-1)(\lambda |\overline {w}^{\prime}| )^{n-1}=O(\lambda ^{n-1})\ll 1$ holds almost everywhere except near the pipe centre. The $\mu -\overline {\mu }$ field exhibits a complex structure, and it is not immediately clear how it contributes to the integral in (5.6). Nevertheless, its $\theta$ – $z$ average, $\langle \mu -\bar {\mu }\rangle =\langle \mu \rangle -\bar {\mu }$ , shown in figure 10, at least plays a role in reducing drag, as the figure suggests that $\int ^1_0(\overline {\mu }-\langle \mu \rangle ) \langle D\negthinspace :\negthinspace D\rangle r{\rm d}r\gt 0$ .

Since the definition of ${\textit{Re}}_b$ depends on $\langle \mu _{{wall}}^* \rangle$ (see (2.5)), one may wonder whether sublaminar drag can be observed using other choices of reference viscosity. We performed comparisons using both the average viscosity and $\mu _0^*$ , and found that the conclusion remains unchanged. Drag reduction can be more directly confirmed by keeping the pressure gradient fixed and comparing the resulting flow rates.