2.1. Base flow

The system admits an axisymmetric solution (see figure 2 a) that can be obtained by specifying the parameters $(n,\mu _{\infty },\lambda ,a,Re,\epsilon ,\alpha )$ , where $\alpha =2\pi /L$ . This solution, denoted by $[u,v,w,p]=[\overline {u},0,\overline {w},\overline {p}](r,z)$ , is taken as the base flow for the linear stability analysis. Note that $q$ in (2.1a ) must be adjusted to ensure that $\overline {w}=1$ at the centreline $r=0$ and is therefore not treated as an independent parameter.

We are primarily interested in the case where the amplitude of the wall undulation $\epsilon$ is small. In this regime, the base flow can be expressed as the sum of the smooth-wall base flow $[u,v,w,p]=[0,0,\overline {w}_0(r),0]$ and small modulation terms. The function $\overline {w}_0$ can be obtained by solving the streamwise momentum equation

(2.4a) \begin{align} r^{-1}\big(r\overline {\mu }_0\,\overline {w}_0^{\prime}\big)'&=-q, \end{align}

(2.4b) \begin{align} \text{where}\qquad \overline {\mu }_0(r)&=\mu _\infty +(1-\mu _\infty )\big\{1+\big|\lambda \overline {w}_0^{\prime}\big|^a\big\}^{(n-1)/a}, \\[10pt] \nonumber \end{align}

subject to the boundary conditions $\overline {w}_0(1)=0$ and $\overline {w}_0(0)=1$ . The prime denotes differentiation with respect to $r$ .

The power-law fluid approximation becomes valid when $\lambda$ is large and $\mu _{\infty }\ll O(\lambda ^{n-1})$ :

(2.5) \begin{align} \overline {w}_0\approx 1-r^{1+1/n}, \qquad \overline {\mu }_0\approx \big(\lambda |\overline {w}_0^{\prime}| \big)^{n-1}\approx \lambda ^{n-1}(1+1/n)^{n-1}r^{1-1/n}. \end{align}

Note that this approximation does not hold near the pipe centre. This reflects the well-known singular behaviour of the power-law model, which is one reason the Carreau model ( $a=2$ ) is more preferred in numerical analysis. The Carreau–Yasuda model extends the latter by allowing a general $a$ , which improves the fitting to experimental measurements.

2.2. Perturbation flow

We now introduce three-dimensional perturbations of infinitesimal amplitude $\Delta$ to the base flow:

(2.6a) \begin{align} u&=\overline {u}(r,z)+\Delta \tilde {u}(r,z)e^{\sigma t}\cos (m\theta ), \end{align}

(2.6b) \begin{align} v&=\Delta \tilde {v}(r,z)e^{\sigma t}\sin (m\theta ), \end{align}

(2.6c) \begin{align} w&=\overline {w}(r,z)+\Delta \tilde {w}(r,z)e^{\sigma t}\cos (m\theta ), \end{align}

(2.6d) \begin{align} p&=\overline {p}(r,z)+\Delta \tilde {p}(r,z)e^{\sigma t}\cos (m\theta ), \end{align}

(2.6e) \begin{align} \mu &=\overline {\mu }(r,z)+\Delta \tilde {\mu }(r,z)e^{\sigma t}\cos (m\theta ). \end{align}

Here, $m$ and $\sigma$ are the azimuthal wavenumber and the growth rate of the perturbation, respectively. The perturbations are also assumed to be periodic with length $L$ in the axial direction. We focus on modes with zero frequency, which are usually the most dangerous; (2.6) represents their general form. Note that subharmonic perturbations are not considered in this study. For the cases to be examined in § 5.1, we have verified at a few selected parameter values that they remain stable even if subharmonic perturbations are allowed.

For all flow visualisations, the perturbation is normalised so that the maximum of $|w-\overline {w}|$ becomes unity. Figure 2(b) shows a typical flow structure of the perturbation obtained as a saturated state from DNS of the full GNF equations. The streamwise velocity clearly exhibits a streaky pattern (blue), whereas the azimuthal velocity (green) is wave-like and concentrated near the wall; we note the general similarity of the structures shown in figure 2(b) to the Newtonian-fluid instability modes shown in figure 3 of Loh & Blackburn (Reference Loh and Blackburn2011) and also that the magnitude of the (streaky) streamwise velocity perturbations are generally found to be much larger than those of the azimuthal/radial (roll) perturbations.

2.3. Stability analysis

Our goal is to trace the neutral stability curves in parameter space, defined by where the real part of $\sigma$ is zero. To this end, we employ three different methods, each with its own advantages and limitations.

The first method is DNS of the full GNF equations (2.1) and (2.3) using a spectral element solver (Blackburn, Rudman & Singh Reference Blackburn, Rudman and Singh2025). This approach provides an accurate means of identifying neutral points in parameter space. However, it is computationally expensive and may become impractical for high Reynolds numbers or strongly shear-thinning fluids. The spatial discretisation employs nodal spectral elements, in which the velocity field is represented by high-order Lagrange interpolants defined on Gauss–Lobatto–Legendre points within each element. Our default computational mesh consists of 10 elements in the axial direction and 14 elements in the radial direction, with mesh refinement applied towards the wall to adequately resolve near-wall flow structures. A polynomial order of 7 is used within each element. For high-Reynolds-number computations ( ${\textit{Re}}=10^5$ ), a refined mesh with 17 elements in the radial direction is employed to ensure sufficient resolution.

To compute the growth rate using DNS, one first obtains the base flow by imposing axisymmetry and numerically time integrating the governing equations until a steady state is reached. Stability at a given set of parameters is then determined by time integrating the system from an initial condition given by the base flow with small random perturbations of $O(10^{-4})$ . This velocity amplitude is sufficiently small that nonlinear effects are negligible until the modal growth of the perturbation becomes clearly observable. We verified that DNS yields stability results consistent with those obtained from linear stability analysis for Newtonian fluids (the latter approach was used by Loh & Blackburn Reference Loh and Blackburn2011).

The second method determines the neutral curve by extending the large-Reynolds-number asymptotic theory of Hall & Ozcakir (Reference Hall and Ozcakir2021) to shear-thinning fluids. The resulting approximate formula for the neutral curve serves as the fastest tool for exploring the parameter space. As we shall see in §§ 3 and 4, a relatively simple analytic formula can be obtained when the power-law approximation is valid. However, for more general rheologies, such a clean formula does not exist. Moreover, as we shall see later (just following (3.32)), the asymptotic analysis drops certain terms from the governing equations under the assumption that $O(1/\ln Re)$ is small, although these terms remain significant unless ${\textit{Re}}$ is unrealistically large; this issue also arises in the Newtonian case (Hall & Ozcakir Reference Hall and Ozcakir2021).

The third method, the asymptotic preserved reduction (APR), represents a compromise between the first two approaches. The main idea is to simplify the linearised equations while retaining terms that are important in the high-Reynolds-number asymptotic analysis. Although solving the resulting reduced equations to obtain an approximate neutral curve still requires numerical computation, it is significantly easier than solving the full system. Furthermore, our formulation retains the aforementioned terms involving the relative $O(1/\ln Re)$ factor in the asymptotic expansion, thereby improving the accuracy of the predicted neutral curve.

3.1. Reduction of the base flow equations

Let us write $E=e^{i\alpha z}$ and apply a Fourier expansion in $z$ to the base flow. Higher-order harmonics are expected to be negligible as the wall is perturbed by a single-mode trigonometric function. Thus, the following approximation may hold:

(3.1a) \begin{align} \overline {u}&\approx \epsilon \overline {u}_1(r)E+\text{c.c.}, \end{align}

(3.1b) \begin{align} \overline {w}&\approx \overline {w}_0(r)+\epsilon \overline {w}_1(r)E+\text{c.c.}, \end{align}

(3.1c) \begin{align} \overline {p}&\approx \epsilon \overline {p}_1(r)E+\text{c.c.}, \end{align}

(3.1d) \begin{align} \overline {\mu }&\approx \overline {\mu }_0(r)+\epsilon \overline {\mu }_1(r)E+\text{c.c.}, \end{align}

where c.c. denotes the complex conjugate. In the asymptotic limit the terms with subscript 0 satisfy (2.4). Terms carrying subscript 1 correspond to modulation of the base flow relative to the smooth wall case. Although they are scaled by $\epsilon$ for convenience, it will become evident later through matched asymptotic expansion that an additional correction involving ${\textit{Re}}$ is required to accurately estimate their magnitude.

Substituting (3.1) into the governing equations (2.1) and collecting the terms proportional to $E$ , we obtain

(3.2a) \begin{align} i\alpha \overline {u}_1 \overline {w}_0 +\overline {p}_1^{\prime}&= \frac {1}{\textit{Re}} \overline {\mu }_0 \overline {u}_1^{\prime\prime}, \end{align}

(3.2b) \begin{align} i\alpha \overline {w}_0 \overline {w}_1+ \overline {u}_1 \overline {w}_0^{\prime} + i\alpha \overline {p}_1&= \frac {1}{\textit{Re}} \overline {\mu }_0 \overline {w}_1^{\prime\prime}, \end{align}

(3.2c) \begin{align} \overline {u}^{\prime}_1+\frac {1}{r}\overline {u}_1+i\alpha \overline {w}_1&=0. \end{align}

Here, in the viscous terms, only those involving second derivatives in the radial direction have been retained, as they become important in the near-wall critical layer. Equations (3.2) together with the boundary conditions

(3.3) \begin{align} \overline {u}_1=0, {\quad \overline {u}_1^{\prime}=-i\alpha \gamma }\quad \text{at} \quad r=1, \end{align}

where $\gamma =-\overline {w}_0^{\prime}|_{r=1}$ , constitute the APR system for the base flow modulation. Note that $\overline {w}_0$ and $\overline {\mu }_0$ are assumed to satisfy (2.4). The derivation of (3.3) can be found in Appendix A; this expression captures the leading-order effect of the wavy wall.

As shown in figure 1, in the matched asymptotic expansion analysis, it is necessary to consider both the near-wall critical layer of thickness $O({\textit{Re}}^{-1/3})$ and the wavy-wall layer of thickness $O(\epsilon )$ . Both layers are located near the pipe wall, with the latter layer lying beneath the former layer (it will be shown later that $\epsilon \ll {\textit{Re}}^{-1/3}$ ). The analysis of the wavy wall layer is only needed to find the conditions (3.3) and does not explicitly appear in the subsequent analysis. In the bulk region away from the wall, the base flow modulation satisfies inviscid equations to leading order, whereas within the critical layer, the viscous terms retained in (3.2) become important. That is, the APR solves the flow in both the core region and the critical layer simultaneously.

We now turn our attention to the high-Reynolds-number asymptotic analysis. A consistent expansion in the bulk region takes the following form:

(3.4) \begin{align} \big[\overline {u}_1,\overline {w}_1,\overline {p}_1\big]={\textit{Re}}^{-1/3}\big[\overline {U}_{1},\overline {W}_{1},\overline {P}_{1}\big](r)+{\cdots} . \end{align}

The $O(\epsilon )$ wall undulation induces a streamwise velocity perturbation of the same order within the near-wall critical layer. The flow in this region generates an $O({\textit{Re}}^{-1/3}\epsilon )$ radial perturbation that drives the bulk perturbation; this is why the factor ${\textit{Re}}^{-1/3}$ is introduced in (3.4). Substituting this expansion into the APR equation (3.2), it is straightforward to verify that the viscous terms vanish at leading order, and the Rayleigh equation

(3.5) \begin{align} \overline {w}_0 \left(\overline {U}_1^{\prime\prime}+ \frac {\overline {U}_1^{\prime}}{r} - \frac {\overline {U}_1}{r^2} -\alpha ^2 \overline {U}_1 \right) + \frac {\overline {w}_0^{\prime}}{r}\overline {U}_1 - \overline {w}_0^{\prime\prime}\overline {U}_1=0 \end{align}

is obtained through the standard reduction of the governing equation to a single equation for the wall-normal perturbation velocity (see e.g. Drazin Reference Drazin2002). From (2.3), we have the expansion

(3.6) \begin{align} \overline {\mu }_1= {\textit{Re}}^{-1/3}\frac {\check {\mu }}{\overline {w}_0^{\prime}}\overline {W}_1^{\prime}+{\cdots} , \end{align}

where

(3.7) \begin{align} \check {\mu }(r)=(n-1)(1-\mu _\infty )\big\{1+\big|\lambda \overline {w}_0^{\prime}\big|^a\big\}^{(n-1-a)/a}\big|\lambda \overline {w}_0^{\prime}\big|^a. \end{align}

The need for the ${\textit{Re}}^{-1/3}$ factor in (3.4) becomes clearer upon analysing the near wall critical layer introducing the stretched variable

(3.8) \begin{align} Y={\textit{Re}}^{1/3}(1-r). \end{align}

The base flow modulation in this region has an asymptotic expansion

(3.9a) \begin{align} \overline {u}_1&=-{\textit{Re}}^{-1/3}\overline {U}_{c1}(Y)+{\cdots} , \end{align}

(3.9b) \begin{align} \overline {w}_1&=\overline {W}_{c1}(Y)+{\cdots} , \end{align}

(3.9c) \begin{align} \overline {p}_1&={\textit{Re}}^{-1/3}\overline {P}_{c1}(Y)+{\cdots} . \end{align}

Noting that the Taylor expansions $\overline {w}_0={\textit{Re}}^{-1/3}\gamma Y+{\cdots}$ and $\overline {\mu }_0=b+{\cdots}$ hold with constants $\gamma =-\overline {w}_0^{\prime}|_{r=1}$ and $b=\overline {\mu }_0|_{r=1}$ , the leading order equations can be found as

(3.10a) \begin{align} \partial _Y\overline {P}_{c1} &=0 , \end{align}

(3.10b) \begin{align} i\alpha \gamma Y \overline {W}_{c1}+ \gamma \overline {U}_{c1} +i\alpha \overline {P}_{c1}&= b \partial _Y^2\overline {W}_{c1}, \end{align}

(3.10c) \begin{align} \partial _Y\overline {U}_{c1}+i\alpha \overline {W}_{c1}&=0. \\[10pt] \nonumber \end{align}

At the bottom of the critical layer, we must impose the wavy-wall boundary conditions (3.3), which now translate into

(3.11) \begin{align} \overline {U}_{c1}=0,\quad {\partial _Y\overline {U}_{c1}=-i\alpha \gamma } \quad \text{at} \quad Y=0. \end{align}

After some algebra, the solution of (3.10)–(3.11) can be found as

(3.12) \begin{align} {\overline {U}_{c1}(Y)=-3(i\alpha \gamma )^{2/3}b^{1/3}\int _0^{\xi }\int ^{\infty }_{s}\text{Ai}(\tau )\,{\rm d}\tau\,{\rm d}s},\quad \xi =\left(\frac {i\alpha \gamma }{b}\right)^{1/3}Y, \end{align}

where $\text{Ai}$ is the Airy function (see Appendix B for derivation). As $Y\rightarrow \infty$ ,

(3.13) \begin{align} \overline {U}_{c1}\sim i^{2/3}u_{\infty }, {\qquad u_{\infty }= 3(\alpha \gamma )^{2/3}b^{1/3}\text{Ai}'(0)}. \end{align}

This implies that (3.5) must be solved with the boundary condition $\overline {U}_{1}=-i^{2/3}u_{\infty }$ at $r=1$ ; the matching here is successful owing to the coefficient ${\textit{Re}}^{-1/3}$ in (3.4).

A composite solution can be constructed using the leading-order solutions in the bulk region and the critical layer. For example, an approximation of the radial velocity can be obtained by substituting

(3.14) \begin{align} \overline {u}_1={\textit{Re}}^{-1/3}\big\{\overline {U}_{1}(r)-\overline {U}_{c1}(Y)+i^{2/3}u_{\infty }\big\} \end{align}

into (3.1a ). The Rayleigh equation (3.5) does not admit an analytic solution, but $\overline {U}_1(r)$ can be computed numerically, e.g. using a Chebyshev collocation method.

Note that, from the previous discussion, it is clear which asymptotic expansion should be used for asymptotic analysis of the full governing equations (2.1). For example, the expansion for the bulk flow can be obtained by combining (3.1) and (3.4):

(3.15a) \begin{align} \overline {w}&=\overline {w}_0(r)+\epsilon {\textit{Re}}^{-1/3}\overline {W}_1(r)E+\text{c.c.}+{\cdots} , \end{align}

(3.15b) \begin{align} \overline {u}&=\epsilon {\textit{Re}}^{-1/3}\overline {U}_1(r)E+\text{c.c.}+{\cdots} , \end{align}

(3.15c) \begin{align} \overline {p}&=\epsilon {\textit{Re}}^{-1/3}\overline {P}_1(r)E+\text{c.c.}+{\cdots} . \end{align}

3.2. Reduction of the perturbation equations

Fourier expanding the perturbation fields in the $z$ -direction and discarding the harmonics yields

(3.16a) \begin{align} \tilde {u}&\approx {\textit{Re}}^{-1}\tilde {u}_0(r)+\epsilon \tilde {u}_1(r)E+\text{c.c.}, \end{align}

(3.16b) \begin{align} \tilde {v}&\approx {\textit{Re}}^{-1}\tilde {v}_0(r)+\epsilon \tilde {v}_1(r)E+\text{c.c.}, \end{align}

(3.16c) \begin{align} \tilde {w}&\approx \tilde {w}_0(r)+\epsilon \tilde {w}_1(r)E+\text{c.c.}, \end{align}

(3.16d) \begin{align} \tilde {p}&\approx {\textit{Re}}^{-2}\tilde {p}_0(r)+\epsilon \tilde {p}_1(r)E+\text{c.c.}, \end{align}

(3.16e) \begin{align} \tilde {\mu }&\approx \tilde {\mu }_0(r)+\epsilon \tilde {\mu }_1(r)E+\text{c.c.}, \end{align}

where the subscript 0 denotes the vortex component and the subscript 1 denotes the wave component. The coefficients ${\textit{Re}}^{-1}$ and $\epsilon$ reflect the approximate magnitude of the coefficients; however, as in the case of the base flow, a precise estimation requires matched asymptotic expansion analysis. Among the vortex velocity components, the streamwise component is referred to as the streak, while the cross-streamwise components are called the roll.

We substitute the decomposition (2.6) into the governing equations (2.1), and then further substitute (3.1) and (3.16). Terms proportional to $E$ , $\Delta E$ , $\varDelta$ correspond to the equations for the base flow modulation, the wave and the vortex, respectively. We scale the growth rate as $\sigma =\sigma _0/Re$ so that the unsteady term balances in the vortex part.

The wave equations are simplified similarly to those for the base flow modulation; higher-order Fourier harmonics are discarded and only the viscous terms involving second-order radial derivatives are retained. Furthermore, terms involving $\tilde {u}_0$ , $\tilde {v}_0$ and $\sigma _0$ are neglected, as they are expected to be small owing to the scaling factor $1/Re$ . Under these assumptions, the wave part of the APR equations can be obtained as

(3.17a) \begin{align} i \alpha \overline {w}_0 \tilde {u}_1 +i \alpha \tilde {w}_0 \overline {u}_1 +\tilde {p}_1^{\prime}&= \frac {1}{\textit{Re}} \bigg ( \overline {\mu }_0\tilde {u}_1^{\prime\prime} +\tilde {\mu }_0\overline {u}_1^{\prime\prime} \bigg ), \end{align}

(3.17b) \begin{align} i\alpha \overline {w}_0 \tilde {v}_1 -\frac {m \tilde {p}_1}{r}&= \frac {1}{\textit{Re}} \overline {\mu }_0 \tilde {v}_1^{\prime\prime} , \end{align}

(3.17c) \begin{align} i \alpha \big(\tilde {w}_0 \overline {w}_1 + \tilde {w}_1 \overline {w}_0\big) + \tilde {u}_1 \overline {w}_0^{\prime} + \overline {u}_1 \tilde {w}_0^{\prime} +i\alpha \tilde {p}_1 &= \frac {1}{\textit{Re}} \Big ( \overline {\mu }_0 \tilde {w}_1^{\prime\prime}+\tilde {\mu }_0 \overline {w}_1^{\prime\prime} \Big ), \end{align}

(3.17d) \begin{align} i\alpha \tilde {w}_1+\tilde {u}_1^{\prime}+r^{-1}\tilde {u}_1+mr^{-1}\tilde {v}_1&=0, \end{align}

where $ \tilde {\mu }_0=(\tilde {w}^{\prime}_0/\overline {w}_0^{\prime})\check {\mu }$ . The boundary conditions for these equations can be readily obtained as

(3.18) \begin{align} \tilde {u}_1=\tilde {v}_1=0, \qquad {\tilde {u}_1^{\prime}=-i\alpha \varGamma }\qquad \text{at}\ r=1, \end{align}

where $\varGamma =-\tilde {w}_0^{\prime}|_{r=1}$ , using an argument similar to that in Appendix A.

The roll and streak parts of the APR equations are obtained as

(3.19a) \begin{align} &\sigma _0 \tilde {u}_0 -\overline {\mu }_0\left(\tilde {u}_0^{\prime\prime}+r^{-1}\tilde {u}_0^{\prime}-r^{-2}\tilde {u}_0-r^{-2}m^2 \tilde {u}_0-2r^{-2}m\tilde {v}_0\right)\nonumber \\&\quad -2\overline {\mu }_0^{\prime}\tilde {u}_0^{\prime} +\tilde {p}_0^{\prime}=-{\textit{Re}}^2\epsilon ^2f_u, \end{align}

(3.19b) \begin{align} &\sigma _0 \tilde {v}_0 -\overline {\mu }_0\left(\tilde {v}_0^{\prime\prime}+r^{-1}\tilde {v}_0^{\prime} -r^{-2}\tilde {v}_0 -r^{-2}m^2 \tilde {v}_0-2r^{-2}m\tilde {u}_0\right)\nonumber \\&\quad -\overline {\mu }_0^{\prime}\left(\tilde {v}_0^{\prime}-r^{-1}\tilde {v}_0-mr^{-1}\tilde {u}_0\right)-r^{-1}m\tilde {p}_0=-{\textit{Re}}^2\epsilon ^2f_v, \end{align}

(3.19c) \begin{align} &\qquad\qquad\qquad\tilde {u}_0^{\prime}+r^{-1}\tilde {u}_0+r^{-1}m\tilde {v}_0=0 \end{align}

and

(3.19d) \begin{align} \sigma _0 \tilde {w}_0 -\overline {\mu }_0\big(\tilde {w}_0^{\prime\prime}+r^{-1}\tilde {w}_0^{\prime} -r^{-2}m^2 \tilde {w}_0\big)-\overline {\mu }^{\prime}_0\tilde {w}_0^{\prime} -\tilde {\mu }_0\big(\overline {w}_0^{\prime\prime}+r^{-1}\overline {w}_0^{\prime}\big)-\tilde {\mu }_0^{\prime}\overline {w}_0^{\prime} +\overline {w}_0^{\prime} \tilde {u}_0=0, \end{align}

respectively, where $f_u$ and $f_v$ represent the Reynolds stress terms. By retaining only the formally leading-order terms, we find

(3.19e) \begin{align} f_u &=\left(2r^{-1} \left(r\overline {u}_1^{\dagger } \tilde {u}_{1}\right)^{\prime}+r^{-1}m \overline {u}_{1}^{\dagger }\tilde {v}_1 \right)+\text{c.c.}, \end{align}

(3.19f) \begin{align} f_v&=r^{-2}\left(r^2\overline {u}_1^{\dagger }\tilde {v}_1\right)^{\prime}+\text{c.c.}, \end{align}

indicating that the wave and vortex are coupled. The boundary conditions are

(3.20) \begin{align} \tilde {u}_0=\tilde {u}_0^{\prime}=\tilde {v}_0=0 \quad \text{at} \quad r=1. \end{align}

The asymptotic analysis proceeds as follows. In the bulk region, the wave fields are expanded as

(3.21) \begin{align} \big[\tilde {u}_1,\tilde {v}_1,\tilde {w}_1,\tilde {p}_1\big]={\textit{Re}}^{-1/3}\big[\tilde {U}_{1},\tilde {V}_{1},\tilde {W}_{1},\tilde {P}_{1}\big](r)+{\cdots} . \end{align}

On the other hand, since the magnitude of the vortex field is already accurately captured in (3.16), we express the leading-order terms directly using $\tilde {u}_0,\tilde {v}_0,\tilde {w}_0,\tilde {p}_0$ . The viscous terms in (3.17) can be neglected, yielding

(3.22a) \begin{align} \overline {w}_0 \left(\tilde {U}_1^{\prime\prime} + \frac {\tilde {U}_1^{\prime}}{r}- \frac {\tilde {U}_1}{r^2}-L\tilde {U}_1 - \frac {2m}{r^2} \tilde {V}_1\right) + \overline {w}_0^{\prime} \frac {\tilde {U}_1}{r} -\overline {w}_0^{\prime\prime} \tilde {U}_1 \nonumber \\ +\tilde {w}_0 \left( \overline {U}_1^{\prime\prime}+ \frac {\overline {U}_1^{\prime}}{r} - \frac {\overline {U}_1}{r^2} -L\overline {U}_1\right) + \tilde {w}_0^{\prime} \frac {\overline {U}_1}{r}-\tilde {w}_0^{\prime\prime} \overline {U}_1=0 , \end{align}

(3.22b) \begin{align} \overline {w}_0 \left (\tilde {U}_1^{\prime} + \frac {\tilde {U}_1}{r}+\frac {rL}{k} \tilde {V}_1\right ) -\overline {w}_0^{\prime}\tilde {U}_1 + \tilde {w}_0 \left(\overline {U}_1^{\prime}+\frac {\overline {U}_1}{r} \right) - \tilde {w}_0^{\prime}\overline {U}_1 =0, \end{align}

where $L=\alpha ^2+m^2/r^2$ . They are the inviscid linear stability equations for the unmodulated base flow, with additional inhomogeneous terms. As in the Rayleigh–Craik equation studied in oceanography (see e.g. Craik Reference Craik1982; Phillips & Shen Reference Phillips and Shen1996), the inhomogeneous terms arise from the interaction of the base flow modulation and the streak perturbation.

Inside the near-wall critical layer, we switch the wave expansions to

(3.23a) \begin{align} \tilde {u}_1&=-{\textit{Re}}^{-1/3}\tilde {U}_{c1}(Y)+{\cdots} , \end{align}

(3.23b) \begin{align} \tilde {v}_1&=\tilde {V}_{c1}(Y)+{\cdots} , \end{align}

(3.23c) \begin{align} \tilde {w}_1&=\tilde {W}_{c1}(Y)+{\cdots} , \end{align}

(3.23d) \begin{align} \tilde {p}_1&={\textit{Re}}^{-1/3}\tilde {P}_{c1}(Y)+{\cdots} , \end{align}

and write $\overline {w}_0={\textit{Re}}^{-1/3}\gamma Y+{\cdots}$ , $\overline {\mu }_0=b+{\cdots}$ , $\tilde {w}_0={\textit{Re}}^{-1/3}\varGamma Y+{\cdots}$ and $\tilde {\mu }_0=({\varGamma }/{\gamma })\beta b+{\cdots}$ with $\beta =(\check {\mu }/\overline {\mu }_0)|_{r=1}$ . The leading order equations are then obtained as

(3.24a) \begin{align} \partial _Y\tilde {P}_{c1}&=0, \end{align}

(3.24b) \begin{align} i\alpha \gamma Y \tilde {V}_{c1} -m \tilde {P}_{c1}&=b \partial _Y^2\tilde {V}_{c1}, \end{align}

(3.24c) \begin{align} i \alpha \big(\varGamma Y \overline {W}_{c1} + \gamma Y \tilde {W}_{c1}\big) + \gamma \tilde {U}_{c1} + \varGamma \overline {U}_{c1} +i\alpha \tilde {P}_{c1} &= b \partial _Y^2\tilde {W}_{c1}+\frac {\varGamma }{\gamma }\beta b \partial _Y^2\overline {W}_{c1}, \end{align}

(3.24d) \begin{align} \partial _Y\tilde {U}_{c1}+m\tilde {V}_{c1}+i\alpha \tilde {W}_{c1}&=0, \end{align}

subject to the boundary conditions

(3.25) \begin{align} \tilde {U}_{c1}=\tilde {V}_{c1}=0,\qquad {\partial _Y\tilde {U}_{c1}=-i\alpha \varGamma } \qquad \text{at} \ Y=0. \end{align}

Solving (3.24)–(3.25) yields (see Appendix B)

(3.26) \begin{align} \tilde {U}_{c1}(Y)&=-(i\alpha \gamma )^{2/3}b^{1/3} \frac {\varGamma (1-\beta )}{\gamma }\left (3\int _0^{\xi } \int ^{\infty }_{s}\text{Ai}(\tau )\,{\rm d}\tau\, {\rm d}s-\text{Ai}'(\xi )+\text{Ai}'(0)\right )\!, \end{align}

(3.27) \begin{align} \tilde {V}_{c1}(Y)&=-\left (\frac {b}{i\alpha \gamma }\right )^{2/3}\frac {m\tilde {P}_{c1}}{b}S(\xi ), \end{align}

with the asymptotic behaviour

(3.28) \begin{align} \tilde {U}_{c1}&\sim i^{2/3}\frac {2\varGamma }{3\gamma } (1-\beta ) u_{\infty }, \end{align}

(3.29) \begin{align} \tilde {V}_{c1} &\sim \frac {m\tilde {P}_{c1}}{i\alpha \gamma Y}, \\[10pt] \nonumber \end{align}

as $Y\rightarrow \infty$ . Here, $S(\xi )$ denotes the Scorer function (see e.g. Olver Reference Olver1974). Using equations (B3), (3.13), (3.28), the pressure can be found as

(3.30) \begin{align} \tilde {P}_{c1}=\frac {i\alpha }{\alpha ^2+m^2}\lim _{Y\rightarrow \infty }\big(\varGamma \overline {U}_{c1}+\gamma \tilde {U}_{c1}\big)= \frac {i^{5/3}\alpha \varGamma }{\alpha ^2+m^2} \frac {5-2\beta }{3}u_{\infty }. \end{align}

As (3.29) indicates, the wave amplitude grows as the near-wall critical layer is approached, and so does the Reynolds stress appearing in the vortex equations. This effect manifests as a non-zero azimuthal velocity of the roll component at the wall, which we denote hereafter by $\mathcal{V}=\tilde {v}_0|_{r=1}$ .

The constant $\mathcal{V}$ can be found as follows. From the previous analysis, it is evident that the Reynolds stress terms in the vortex equations (3.19) become significant within the critical layer. Therefore, $\epsilon$ must be selected so that a proper asymptotic balance is maintained in that region. Using (3.13), (3.29) and (3.19f ), we obtain $f_v=F_v(Y)+{\cdots}$ in the critical layer, with asymptotic behaviour

(3.31) \begin{align} F_v(Y)\sim -\frac {G}{Y^2}, \qquad G\equiv \frac {m}{ \gamma }\frac {\varGamma }{\alpha ^2+m^2} \frac {10-4\beta }{3}u_{\infty }^2 \end{align}

as $Y\rightarrow \infty$ . Azimuthal component of the roll equation (3.19b ) reduces to $b\partial _Y^2\tilde {V}_{0c}=F_v$ , where we used the expansion $\tilde {v}_0={\textit{Re}}^{-2/3}(\epsilon Re)^2\tilde {V}_{0c}(Y)+{\cdots}$ . This implies that $\tilde {V}_{0c}$ behaves like $(G/b)\ln Y+\textrm {const.}=(G/b)\ln {\textit{Re}}^{1/3}+{\cdots}$ for large $Y$ . Hence, we conclude

(3.32) \begin{align} \mathcal{V}=\mathscr{R}\frac {G}{3b},\quad \text{where} \ \mathscr{R}\equiv \epsilon ^2 {\textit{Re}}^{4/3}\ln Re, \end{align}

discarding the terms that are smaller by a relative factor of $O(1/\ln Re)$ . Note also that the combined parameter $\mathscr{R}$ (denoted by $R^*$ by Hall & Ozcakir Reference Hall and Ozcakir2021) must be $O(1)$ and this requirement determines

(3.33) \begin{align} \epsilon =O({\textit{Re}}^{-2/3}/\sqrt {\ln \textit{Re}})\ll O\big({\textit{Re}}^{-1/3}\big). \end{align}

In the bulk region, the leading-order vortex problem is given by

(3.34a) \begin{align} \sigma _0 \tilde {u}_0 -\overline {\mu }_0\left(\tilde {u}_0^{\prime\prime}+r^{-1}\tilde {u}_0^{\prime}-r^{-2}\tilde {u}_0-r^{-2}m^2 \tilde {u}_0-2r^{-2}m\tilde {v}_0\right)\nonumber \\ -2\overline {\mu }_0^{\prime}\tilde {u}_0^{\prime} +\tilde {p}_0^{\prime}=0, \end{align}

(3.34b) \begin{align} \sigma _0 \tilde {v}_0 -\overline {\mu }_0\left(\tilde {v}_0^{\prime\prime}+r^{-1}\tilde {v}_0^{\prime} -r^{-2}\tilde {v}_0 -r^{-2}m^2 \tilde {v}_0-2r^{-2}m\tilde {u}_0\right)\nonumber \\ -\overline {\mu }_0^{\prime}\left(\tilde {v}_0^{\prime}-r^{-1}\tilde {v}_0-mr^{-1}\tilde {u}_0\right)-r^{-1}m\tilde {p}_0=0, \end{align}

(3.34c) \begin{align} \tilde {u}_0^{\prime}+r^{-1}\tilde {u}_0+r^{-1}m\tilde {v}_0=0, \end{align}

(3.34d) \begin{align} \sigma _0 \tilde {w}_0 -\overline {\mu }_0\left(\tilde {w}_0^{\prime\prime}+r^{-1}\tilde {w}_0^{\prime} -r^{-2}m^2 \tilde {w}_0\right)-\overline {\mu }^{\prime}_0\tilde {w}_0^{\prime}\nonumber \\ -\tilde {\mu }_0\left(\overline {w}_0^{\prime\prime}+r^{-1}\overline {w}_0^{\prime}\right)-\tilde {\mu }_0^{\prime}\overline {w}_0^{\prime} +\overline {w}_0^{\prime} \tilde {u}_0=0, \end{align}

subject to the boundary conditions

(3.35) \begin{align} \tilde {u}_0=\tilde {w}_0=0,\qquad \tilde {v}_0=\mathcal{V}\qquad \text{at} \ r=1. \end{align}

Solving this problem with the slip velocity $\mathcal{V}$ given by (3.32) determines the growth rate  $\sigma _0$ .

4.1. Numerical method for APR

Recall that the APR system consists of equations for the base flow modulation (3.2)–(3.3), wave (3.17)–(3.18) and vortex (3.19)–(3.20) parts. These form an eigenvalue problem with $\sigma _0$ as the eigenvalue, which can be solved once the parameters $(n,\mu _{\infty },\lambda ,a,Re,\epsilon ,\alpha ,m)$ are specified. However, if one is only interested in computing the neutral stability curve, it is far more efficient to solve for $\epsilon$ that yields $\sigma _0=0$ while fixing other parameters, as described later.

We first note that (3.2) governing the base flow modulation can be combined into a single equation. Using the linear operator $\mathscr{L}$ defined in (C1), this equation can be compactly written in the form

(4.1) \begin{align} \mathscr{L}\,\overline {u}_1=0. \end{align}

Substituting the expansion

(4.2) \begin{align} {\overline {u}_1= \frac {i\alpha \gamma }{2} r\big(1-r^2\big)+\sum _{l=0}^N \hat {X}_lr\big(1-r^2\big)^2T_{2l}(r)} \end{align}

into (4.1) and evaluating at the collocation points $r_k=\cos ((k+1)\pi /(2N+4)), k=0,1,\ldots ,N$ , we have a linear algebraic problem that can be solved numerically for $\hat {X}_l$ . Here, $T_l(r)$ is the $l$ th Chebyshev polynomial. Note that (4.2) is designed to automatically satisfy the boundary conditions (3.3). There is freedom in the choice of the function used in the first term (i.e. $r(1-r^2)$ ); other functions could be acceptable provided they satisfy the regularity condition at $r=0$ and the boundary conditions at $r=1$ .

Having determined the base flow modulation $\overline {u}_1$ , we now proceed to the perturbation analysis. The wave, roll and streak parts of the equations can be written in the form

(4.3a) \begin{align} \left [\,\, \begin{array}{c} \tilde{ \kern-9pt\mathscr{L}}_{1u}\tilde {u}_{1}+ \, \, \, \tilde{ \kern-9pt\mathscr{L}}_{1v}\tilde {v}_1\\ \tilde{ \kern-9pt\mathscr{L}}_{2u}\tilde {u}_1+ \,\,\, \tilde{ \kern-9pt\mathscr{L}}_{2v}\tilde {v}_1 \end{array} \!\!\right ] &= \left [ \!\!\begin{array}{c} \mathscr{B}_1(\overline {u}_1,\tilde {w}_0)\\ \mathscr{B}_2(\overline {u}_1,\tilde {w}_0) \end{array}\!\! \right ]\! , \end{align}

(4.3b) \begin{align} \tilde{ \kern-9pt\mathscr{L}}_r\tilde {u}_0&= \epsilon ^2\big\{\mathscr{B}_u\big(\overline {u}_1,\tilde {u}_1\big)+\mathscr{B}_v\big(\overline {u}_1,\tilde {v}_1\big)\big\}, \end{align}

(4.3c) \begin{align} \tilde{ \kern-9pt\mathscr{L}}_s\tilde {w}_0&=\overline {w}_0^{\prime}\tilde {u}_0, \end{align}

using the linear operators $\ \tilde{ \kern-9pt\mathscr{L}}_{1u}, \ \tilde{ \kern-9pt\mathscr{L}}_{1v}, \ \tilde{ \kern-9pt\mathscr{L}}_{2u}, \ \tilde{ \kern-9pt\mathscr{L}}_{2v}, \ \tilde{ \kern-9pt\mathscr{L}}_r, \ \tilde{ \kern-9pt\mathscr{L}}_s$ and the bilinear operators $\mathscr{B}_1,\mathscr{B}_2,\mathscr{B}_u,\mathscr{B}_v$ (see Appendix C). The abstract formulation presented here highlights the physical mechanism of the streaky instability. Wall undulation creates a base flow modulation, which interacts with the streak and influences the inviscid stability of the wave. Then, the interaction between the wave and the base flow modulation generates Reynolds stresses that drive the rolls, which in turn amplify the streaks via the lift-up mechanism.

For computational convenience, we assume that the eigenfunction is normalised such that $\tilde {w}_0^{\prime}|_{r=1}=-\varGamma$ is unity. The boundary conditions (3.18) and (3.20) are satisfied using the expansions

(4.4a) \begin{align} \tilde {u}_1&= -\frac {i\alpha }{2}r^{j-1}\big(1-r^2\big)+\sum _{l=0}^N \hat {Y}_lr^{j-1}\big(1-r^2\big)^2T_{2l}(r), \end{align}

(4.4b) \begin{align} \tilde {v}_1&=\sum _{l=0}^N \hat {Z}_lr^{j-1}\big(1-r^2\big)T_{2l}(r), \end{align}

(4.4c) \begin{align} \epsilon ^{-2}\tilde {u}_0&=\sum _{l=0}^L \hat {U}_lr^{j-1}\big(1-r^2\big)^2T_{2l}(r), \end{align}

(4.4d) \begin{align} \epsilon ^{-2}\tilde {w}_0&=\sum _{l=0}^L \hat {W}_lr^{j}\big(1-r^2\big)T_{2l}(r), \end{align}

where $j=1$ if $m=1$ ; $j=2$ if $m$ is even; and $j=3$ if $m$ is odd and $m\neq 1$ . The first few terms in the Taylor expansion of the previous basis functions are compatible with the regularity conditions at the pipe centreline, similar to those used by Deguchi & Walton (Reference Deguchi and Walton2013). Note that the assumption $\tilde {w}_0^{\prime}|_{r=1}=1$ implies that

(4.5) \begin{align} \epsilon ^{-2}=-2\sum _{l=0}^L \hat {W}_l. \end{align}

We apply the following iterative scheme to the discretised system to determine the value of $\epsilon$ that satisfies (4.5).

1. (i) Assume that $\epsilon ^2$ and the coefficients $\hat {W}_l$ are given. Then, the right-hand side of (4.3a ) is known, and the resulting linear algebraic system obtained after discretisation can be solved for $\hat {Y}_l$ and $\hat {Z}_l$ .

2. (ii) Solve the discretised version of (4.3b ) for $\hat {U}_l$ , using $\hat {Y}_l$ and $\hat {Z}_l$ found in step (i).

3. (iii) Solve the discretised version of (4.3c ) for $\hat {W}_l$ , using $\hat {U}_l$ found in step (ii). Evaluate $\epsilon ^{-2}$ from (4.5).

4. (iv) Compare the $L^2$ norm of the difference between the $\hat {W}_l$ assumed in step (i) and that obtained in step (iii). If the norm is less than $10^{-9}$ , the computation is terminated. Otherwise, return to step (i) updating the values of $\hat {W}_l$ and $\epsilon ^{-2}$ .

The converged value of $\epsilon$ corresponds to the linear critical point and the associated eigenfunction can be constructed by $\hat {Y}_l,\hat {Z}_l,\hat {U}_l,\hat {W}_l$ .

4.2. Analytic neutral curve for the power-law fluids

Combining (3.31), (3.13) and (3.32), we obtain the slip velocity in (3.35) as determined by the wave part of the asymptotic analysis:

(4.6) \begin{align} \frac {\mathcal{V}}{\varGamma }=\frac {m\alpha ^{4/3} \gamma ^{1/3}(10-4\beta )(\text{Ai}'(0))^2}{b^{1/3}(\alpha ^2+m^2)}\mathscr{R}. \end{align}

Here, $\text{Ai}'(0)\approx -0.2588194$ . The vortex equations (3.34), with $\sigma _0=0$ and subject to the boundary conditions (3.35), admit a solution only when $\mathcal{V}/\varGamma$ takes a specific value. Substituting this value into (4.6) determines $\mathscr{R}$ , from which we can obtain the neutral value of $\epsilon$ .

In general, the vortex equations are ordinary differential equations with complicated coefficients and therefore require numerical treatment. However, in the power-law fluid limit (see § 2.1), analytical progress can be made, since the viscosity coefficients are simplified as

(4.7) \begin{align} \overline {\mu }_0={\textit{Kr}}^{2s},\qquad \check {\mu }=(n-1){\textit{Kr}}^{2s},\qquad \tilde {\mu }_0=K\frac {(1-n)}{(1+1/n)}r^{1-2/n} \frac {{\rm d}\tilde {w}_0}{{\rm d}r}, \end{align}

where

(4.8) \begin{align} s= (1-1/n)/2,\qquad K= \lambda ^{n-1}(1+1/n)^{n-1}. \end{align}

Noting $\overline {w}_0=1-r^{1+1/n}$ , the neutral solution of (3.34) that remains finite at the channel centre and satisfy the boundary conditions (3.35) can be found to be

(4.9a) \begin{align} \tilde {u}_0(r)&=\frac {\mathcal{V}}{\zeta -\eta } m\big(r^{\eta }-r^{\zeta }\big), \tilde {v}_0(r)=\frac {\mathcal{V}}{\zeta -\eta } \big\{(\zeta +1)r^{\zeta }-(\eta +1)r^{\eta }\big\}, \end{align}

(4.9b) \begin{align} \tilde {w}_0(r)&=\frac {\mathcal{V}}{K} \frac {m(1+1/n)}{\zeta -\eta } \left\{\frac {r^{\zeta +1+2/n}-r^{\chi }}{J_{\zeta }}-\frac {r^{\eta +1+2/n}-r^{\chi }}{J_{\eta }}\right\}\!, \end{align}

where $\zeta ,\eta ,\chi , J_{\zeta }$ and $J_{\eta }$ are functions of $m$ and $n$ as given in Appendix D. In the pipe centre region, where $r=O(\lambda ^{-n})$ , the power-law approximation is not valid (see Appendix B of He et al. Reference He, Deguchi, Song and Blackburn2025) and thus requires separate treatment. However, this region is passive and does not affect the leading order solution.

Using (4.9b ) to compute $\varGamma =-\tilde {w}_0^{\prime}|_{r=1}$ yields $\mathcal{V}/\varGamma$ in the form given in (D6). Equating this with the right-hand side of (4.6) and substituting $\beta =(n-1)$ , $b=K$ , $\gamma =(1+1/n)$ yields a formula for determining $\mathscr{R}$ as a function of $m$ , $n$ , $\alpha$ and $\lambda$ . Using the expression for $\mathscr{R}$ in (3.32), the neutral value of $\epsilon$ can be obtained as a function of $m,n,\alpha , \lambda ,$ and ${\textit{Re}}$ .

For power-law fluids, the neutral value of $\epsilon$ can be expressed in terms of the reduced set of parameters ${\textit{Re}}/K$ , $m,n,$ and $\alpha$ , since the viscosity in (4.7) is proportional to $K$ and this implies that the effective Reynolds number of the problem is ${\textit{Re}}/K$ . To obtain the desired expression for the stability result, we use the definition

(4.10) \begin{align} \mathscr{R}=\epsilon ^2 {\textit{Re}}^{4/3}\ln ({\textit{Re}}/K) \end{align}

instead of that given in (3.32). Adopting this $\mathscr{R}$ implies that we are partially restoring the higher-order terms omitted just before (3.32); note that this does not compromise the accuracy of the leading order approximation. Equating the two different expressions for $\mathcal{V}/\varGamma$ then gives

(4.11) \begin{align} \epsilon ^2 ({\textit{Re}}/K)^{4/3}\ln ({\textit{Re}}/K)=\mathscr{F}(m,n)\frac {\alpha ^2+m^2}{m^2\alpha ^{4/3}}, \end{align}

where

(4.12) \begin{align} \mathscr{F}(m,n)\equiv \frac {\big(\text{Ai}'(0)\big)^{-2}(\zeta -\eta )J_{\eta }J_{\zeta }}{(1+1/n)^{4/3}(14-4n)\big\{\eta J_{\zeta }-\zeta J_{\eta }+(1-\chi +2/n)(J_{\zeta }-J_{\eta })\big\}}. \end{align}

5.1. Validation of the theoretical analysis

The black squares in figure 3(a) show the neutral curve for the full Navier–Stokes equations obtained by DNS. A corrugated pipe with wavelength $L=1$ is used (i.e. $\alpha =2\pi$ ). The results up to ${\textit{Re}}\approx 6000$ coincide with those of Loh & Blackburn (Reference Loh and Blackburn2011). The vertical axis uses $2\epsilon$ to match the wall amplitude definition used in that paper (where it is denoted by $\alpha$ ). The black dashed line shows the asymptotic result of Hall & Ozcakir (Reference Hall and Ozcakir2021); the same result is recovered by setting $n=1$ in our analytic expression (4.11). The additional high- ${\textit{Re}}$ DNS data we computed clearly show that the asymptotic analysis captures the overall trend of the neutral curve reasonably well. An even better approximation is provided by the solid line, which represents the APR result.

Figure 3.

Neutral points in the ${\textit{Re}}$ – $2\epsilon$ parameter plane for an almost power-law fluid $(\mu _{\infty },\lambda ,a)=(0,100,2)$ . The wavelength of the wall undulation is set to $L=1$ . (a) Comparison between the Newtonian results ( $n=1$ , black) and those for a mildly shear-thinning fluid ( $n=0.8$ , red). The squares indicate the DNS results, the solid lines show the APR results and the dashed lines represent the analytic results for power-law fluids. All results are computed with $m=3$ . (b) Neutral points of the $n=0.8$ shear-thinning fluid at ${\textit{Re}}=10^4$ for various $m$ . Squares indicate the DNS results, crosses the APR results and open circles the analytic results for power-law fluids.

Interestingly, the asymptotic results approximate the neutral curve from above, whereas the APR results approach it from below. The same trend is observed for the $n=0.8$ shear-thinning fluid, as indicated by the red data in figure 3(a). In both DNS and APR, we used the Carreau–Yasuda parameters $(\mu _{\infty },\lambda ,a)=(0,100,2)$ . The behaviour of this fluid can be well approximated by power-law rheology (see He et al. Reference He, Deguchi, Song and Blackburn2025) so the comparison to the analytic result (4.11) is meaningful.

All the results in figure 3(a) use the azimuthal wavenumber $m=3$ since Loh & Blackburn (Reference Loh and Blackburn2011) showed this wavenumber to be the most unstable mode, except in a region close to the onset ${\textit{Re}}$ of the instability ( ${\textit{Re}}\approx 2280$ ). We confirmed that this conclusion also holds for the $n=0.8$ case by performing the computations shown in figure 3(b). The APR and analytic results exhibit the same qualitative behaviour as the DNS computations. However, while the APR correctly predicts that $m=3$ is the most dangerous mode, the asymptotic result slightly mispredicts it as $m=4$ . In figure 3(a), the APR and asymptotic results for $n=0.8$ provide good predictions when ${\textit{Re}}$ is large, $\epsilon$ is small and (3.33) is satisfied. For larger corrugations, the approximations are no longer valid, which explains the discrepancy with the DNS near the onset ${\textit{Re}}$ ( ${\textit{Re}} \approx 747$ ).

Figure 4.

Comparison of the radial velocity of the base flow at the neutral value of $\epsilon$ at $m=3$ in figure 3(b): (a) DNS ( $C_{\textit{max}}=5.99\times 10^{-4}$ ), (b) APR ( $C_{\textit{max}}=7.03\times 10^{-4}$ ) and (c) analytical asymptotic solution ( $C_{\textit{max}}=3.50\times 10^{-4}$ ).

At $m=3$ of figure 3(b) (i.e. $(n,\mu _{\infty },\lambda ,a,Re,\alpha )=(0.8,0,100,2,10^4,2\pi )$ ), we compare the neutral flow fields in figures 4–6. Figure 4 presents comparison of the base flow fields. Panel (a) shows the radial velocity component $\overline {u}$ , obtained from the DNS computation. Panel (b) is the APR result, where the right-hand side of (3.1a ) is evaluated using $\overline {u}_1$ obtained by solving (4.1). Panel (c) corresponds to the flow field from the asymptotic analysis, where $\overline {u}_1$ is taken as the composite solution (3.14). All panels show qualitatively the same results. Note that we determined $C_{max }$ in the colourbar from the maximum absolute value of the plotted quantities. The slight difference in $C_{max }$ arises from the use of the neutral value of $\epsilon$ specific to each case. A closer look at the figures shows that DNS results treat the wavy wall as is, while the other two methods use the region $r \in [0,1]$ .

Figure 5.

Comparison of the streamwise perturbation velocity. The parameters are the same as figure 4. (a) DNS, (b) APR and (c) analytical asymptotic solution.

Figure 6.

Comparison of the azimuthal perturbation velocity. The parameters are the same as figure 4. Slices at $\theta =0$ are shown. Here and hereafter, $C_{\textit{max}}$ denotes the maximum absolute value of the physical quantity shown in the colour map. (a) DNS ( $C_{\textit{max}}=8.19\times 10^{-3}$ ), (b) APR ( $C_{\textit{max}}=8.22\times 10^{-3}$ ) and (c) analytical asymptotic solution ( $C_{\textit{max}}=6.90\times 10^{-3}$ ).

Figure 5 illustrates the streamwise component of the perturbation velocity, $w-\overline {w}$ . Panel (a) shows the DNS result; although this is the slice at $z=0$ , the flow field is nearly identical for any slice, as also seen in figure 2. This observation is consistent with the analyses in §§ 2 and 3, since combining (2.6c ) and (3.16c ) to obtain the approximation

(5.1) \begin{align} w-\overline {w}\approx \varDelta \{\tilde {w}_0(r)+\epsilon \tilde {w}_1(r)E+\text{c.c.}\}\cos (m\theta ) \end{align}

reveals that the streak component dominates. In panels (b) and (c), we plot the right-hand side of (5.1) with the wave component omitted. Panel (b) is the APR result where $\tilde {w}_0$ is obtained by solving (4.3) numerically. Panel (c) is the analytic solution (4.9b ) found in the asymptotic analysis. Results shown in panels agree extremely well and are almost indistinguishable.

The most convenient quantity for observing the wave component of the perturbation is the azimuthal velocity $v$ . Combining (2.6b ) and (3.16b ), we obtain

(5.2) \begin{align} v\approx \varDelta \big\{{\textit{Re}}^{-1}\tilde {v}_0(r)+\epsilon \tilde {v}_1(r)E+\text{c.c.}\big\}\sin (m\theta ). \end{align}

Here, recall that, ignoring the logarithmic factor, $\epsilon$ is approximately $O({\textit{Re}}^{-2/3})$ ; see (3.33). Furthermore, we know that $\tilde {v}_1$ is $O({\textit{Re}}^{-1/3})$ in the bulk region (see (3.21)) and $O(1)$ in the near wall critical layer (see (3.21)). The plot of $v$ shown in figure 6(a) thus exhibits wave-like perturbations in the near-wall region; this figure corresponds to the DNS result at $\theta =0$ . Panels (b) and (c) present the APR and asymptotic results, respectively, obtained by retaining only the wave component on the right-hand side of (5.2). The $\tilde {v}_1$ in the equation is obtained by (4.4b ) for APR, while for the asymptotic results, the analytical solution (3.27) is used via the leading order approximation of (3.23b ). Note that a composite solution is not required, as the critical layer solution decays towards the bulk region. Overall, all panels are in good agreement. The value of $C_{\textit{max}}$ in the colourbar for the APR result matches very well with that of the DNS result, while the asymptotic result does not. This is expected because the APR includes higher-order correction terms.

5.2. Towards highly shear-thinning fluid flow

In what follows, we shift our focus to uncovering the general stability characteristics. Accordingly, the values of $m$ and $\alpha$ are chosen so that the instability appears at the smallest possible $\epsilon$ .

Figure 7.

Stability prediction by the asymptotic analysis, where $ \mathscr{R}^{{pw}}$ denotes the right-hand side of (4.11). (a) Results for $n=0.8$ . The curves are computed using $m=1,2,\ldots , 12$ , among which $m=3$ yields the smallest possible $ \mathscr{R}^{{pw}}$ . (b) Optimised $\mathscr{R}^{{pw}}$ as a function of $n$ . The values of $m$ that give the optimum are indicated by the numbers (they switch at the circles).

This optimisation is easily done in the analytical asymptotic result. Let us denote the right-hand side of (4.11) by $\mathscr{R}^{{pw}}$ ; note that a small value of $\mathscr{R}^{{pw}}$ indicates that instability occurs at a small $\epsilon$ for the same ${\textit{Re}}/K$ . For all $m$ , the minimum of $\mathscr{R}^{{pw}}$ occurs at $\alpha =2^{1/2}m$ , just as Hall & Ozcakir (Reference Hall and Ozcakir2021) pointed out for the Newtonian case. This is verified graphically in figure 7(a) using $n=0.8$ . The asymptotic result predicts that the most dangerous mode corresponds to $m=3$ and $\alpha = 2^{1/2}\times 3\approx 4.24$ .

The same analysis is repeated for various power-law indices $n$ and the optimised values of $\mathscr{R}^{{pw}}$ are summarised in figure 7(b). The numbers shown in the figure indicate the value of $m$ that gives the optimal condition. As shown in figure 8, the optimal neutral perturbation in highly shear-thinning fluids develops small-scale eddies and streaks within the near-wall region. This is consistent with physical intuition, as the corresponding base flow develops a pronounced near-wall boundary layer (see figure 1a of He et al. Reference He, Deguchi, Song and Blackburn2025, for example). The streak constitutes the leading order component of the streamwise velocity perturbation, while the cross-sectional velocity perturbation components are dominated by the near-wall wave. Thus, in figure 9, we plot $U_{c1}$ and $V_{c1}$ using (3.26) and (3.27), respectively. We again see that the near-wall length scale decreases with reducing $n$ , which explains why instability appears at smaller corrugation amplitudes in figure 7(b) for strongly shear-thinning flows.

Figure 8.

Perturbation streaks obtained at the neutral points in figure 7(b). The format of the plot is the same as figure 5. (a) $n=0.5$ ( $m=4$ ); (b) $n=0.4$ ( $m=5$ ); (c) $n=0.3$ ( $m=6$ ); (d) $n=0.25$ ( $m=7$ ); (e) $n=0.2$ ( $m=8$ ).

Figure 9.

Magnitude of the wave velocity components $U_{c1}$ and $V_{c1}$ as a function of $Y={\textit{Re}}^{1/3}(1-r)$ . The five curves use the same parameter values and normalisation as in figure 8.

For highly shear-thinning fluids, the optimal value of $\mathscr{R}^{{pw}}$ becomes smaller, which implies that the streaky instability is more easily triggered by wall undulations. Of course, this argument implicitly assumes that ${\textit{Re}}/K$ is fixed. This is effectively equivalent to fixing the Reynolds number used in the experimental community (e.g. Escudier et al. Reference Escudier, Poole, Presti, Dales, Nouar, Desaubry, Graham and Pullum2005; Wen et al. Reference Wen, Poole, Willis and Dennis2017),

(5.3) \begin{equation} {\textit{Re}}_b\equiv \frac {\rho ^* (2R^*) U_b^*}{ \mu _{\textit{wall}}^* }\approx \frac {4\textit{Re}}{\overline {\mu }_0|_{r=1}}\int ^1_0 \overline {w}_0r\,{\rm d}r, \end{equation}

based on the bulk velocity $U_b^*$ , the mean dimensional viscosity at the wall $\mu _{\textit{wall}}^*$ and the fluid density $\rho ^*$ . The approximation on the right-hand side of (5.3) neglects the base flow modulation. In the power-law limit, it can be further simplified as

(5.4) \begin{equation} {\textit{Re}}_b\approx \frac {2n\left(1+\frac {1}{n}\right)^{2-n} }{(3n+1)}\lambda ^{1-n}\textit{Re}=\frac {2(n+1) }{(3n+1)}\frac {\textit{Re}}{K} \end{equation}

and the factor $(n+1)/(3n+1)$ depends only weakly on $n$ .

Figure 10.

Comparison of the neutral wall amplitude, where the value of $\alpha$ is optimised for each fixed $m$ : (a) $m=1$ and (b) $m=5$ . The dashed line corresponds to the analytic result (4.11); $\epsilon$ is obtained using $\lambda =100$ and ${\textit{Re}}=10^4$ . The solid line is the APR result with $(\mu _{\infty },\lambda ,a)=(0,100,2)$ and ${\textit{Re}}=10^4$ .

Figure 11.

Perturbation streaks obtained at the neutral point for $n=0.2$ in figure 10(b). The format of the plot is the same as figure 5. (a) Asymptotic result; (b) APR result.

Figure 10 shows a verification of the asymptotic results (dashed) using the APR (solid), with parameters $(\mu _{\infty },\lambda ,a)=(0,100,2)$ and ${\textit{Re}}=10^4$ used for the latter. Panels (a) and (b) correspond to $m=1$ and $m=5$ , respectively; similar qualitative agreement is observed for other values of $m$ as well. The exponents $\zeta$ and $\eta$ , used in the analytic roll-streak flow field (4.9), become complex for certain parameter values. For example, when $m=5$ , this occurs for $n\lt 6^{1/2}/(6^{1/2}+5)\approx 0.3288$ . Nevertheless, the flow fields from both of the analyses agree reasonably well, as shown in figure 11. The optimum value of $\alpha$ for $m=5$ is $\alpha =5\times 2^{1/2}\approx 7.07$ for the asymptotic result (panel a) and $\alpha \approx 8.50$ for the APR result (panel b).

Figure 12(a) examines the effect of $\mu _\infty$ on the optimal $m$ using the APR. The other parameters are fixed at $\lambda =100,\ a=2$ and ${\textit{Re}}=10^4$ , with $m$ and $\alpha$ optimised to minimise $\epsilon$ . Increasing $\mu _\infty$ makes the stability less sensitive to changes in the power-law index $n$ . As remarked in § 2.1, the validity of the power-law approximation requires $\mu _{\infty }\ll O(\lambda ^{n-1})$ . This condition is easily satisfied when $n$ is close to unity (e.g. $100^{n-1}\approx 0.63$ for $n=0.9$ ) but not when $n$ is small (e.g. $100^{n-1}\approx 0.016$ for $n=0.1$ ). For the black line in figure 12(a), we set $\mu _{\infty }=10^{-4}$ which safely satisfies the condition for $n \in (0.1,1)$ . We confirmed that the result is indeed almost identical to the $\mu _{\infty }=0$ case. However, when $\mu _{\infty }$ is $10^{-3}$ or larger, the results depart from the $\mu _{\infty }=0$ case for small $n$ .

The optimal value of $m$ in the black line of figure 12(a) differs slightly from that in figure 7(b) reflecting the limited accuracy of the asymptotic theory. For example, in the Newtonian case $(n=1)$ , the asymptotic result predicts the most dangerous mode to be $m=3$ , while the APR yields $m=2$ ; see figure 12(b) also.

5.3. Results in an experimentally feasible setup

The fluid data used in this section are summarised in table 1. The first row is the 7 % aluminium soap (AS) data reported by Myers (Reference Myers2005), which exhibits very strong shear-thinning behaviour. The second and third rows correspond to the 0.2 % polyacrylamide (PAA) and 0.2 % xanthan gum (XG) solutions, respectively, used by Escudier et al. (Reference Escudier, Poole, Presti, Dales, Nouar, Desaubry, Graham and Pullum2005). In all of these fluids, viscoelasticity is weak and shear-thinning behaviour is dominant.

Table 1.

Carreau–Yasuda model parameters used in § 5.3. The first row corresponds to the 7 % aluminium soap (AS) in decalin and m-cresol reported by Myers (Reference Myers2005); the second and third rows to the aqueous solutions of $0.2\,\%$ polyacrylamide (PAA) and $0.2\,\%$ xanthan gum (XG) listed in table 1 of Escudier et al. (Reference Escudier, Poole, Presti, Dales, Nouar, Desaubry, Graham and Pullum2005). The values of $\rho ^*$ are estimated from the densities of the solvent and solute.

Figure 12.

Stability analysis using APR with ${\textit{Re}}=10^4, \lambda =100$ and $a=2$ . (a) Neutral wall amplitudes $2\epsilon$ for various values of $\mu _{\infty }$ : $10^{-4}$ (black), $10^{-3}$ (purple) and $10^{-2}$ (blue). Both $m$ and $\alpha$ are optimised, and the numbers in the figure indicate the optimal value of $m$ . (b) Crosses show the results with $n = 0.8$ , with $\alpha$ optimised for each $m$ . Circles are the corresponding asymptotic analytic results.

Figure 12. Long description

Panel A: A line graph shows neutral wall amplitudes for various values of n. The x-axis represents n, ranging from 0.1 to 1.0, and the y-axis represents 2 times epsilon, ranging from 10^-4 to 10^-2. Different colored lines represent different values of m: black for m = 2, purple for m = 4, and blue for m = 8. The numbers on the graph indicate the optimal value of m. Panel B: A scatter plot shows crosses representing results with m = 1, with optimised for each n. Circles represent the corresponding asymptotic analytic results. The x-axis represents m, ranging from 1 to 6, and the y-axis represents 2 times epsilon, ranging from 0.002 to 0.007.

When comparing experiments with the stability theory, it is important to note that the Carreau number $\lambda =(U_c^*/R^*)\lambda ^*$ does not remain constant when the base flow is varied. We therefore introduce a new parameter $\varLambda =\lambda /Re=(\lambda ^*\mu _0^*)/(\rho ^* R^{*2})$ . All our computations assume a radius of $R^*=0.05$ m, as used in the experiments of Escudier et al. (Reference Escudier, Poole, Presti, Dales, Nouar, Desaubry, Graham and Pullum2005). Also, following experimental convention, we employ the Reynolds number ${\textit{Re}}_b$ defined in (5.3) to organise the data.

Figure 13.

Stability diagram for 7 % AS flowing through a pipe of radius $R^*=0.05$ m. The APR is used to find the neutral value of $\epsilon$ . The value of $\alpha$ is optimised. The circles indicate the point at which the most dangerous value of $m$ switches. The grey band shows the range ${\textit{Re}}_b \in [1.621\times 10^{5},5.814\times 10^5]$ in which a spiral instability exists for a smooth wall pipe (see figure 7b of He et al. Reference He, Deguchi, Song and Blackburn2025).

Figure 13 presents the stability results obtained for 7 % AS ( $\varLambda =55.17$ ; other parameters are listed in table 1). Here and hereafter, stability results are obtained using the APR, with $\alpha$ optimised to minimise $\epsilon$ . The most dangerous azimuthal wavenumber $m$ varies with ${\textit{Re}}_b$ . The value of $\lambda$ increases with ${\textit{Re}}_b$ , implying that the flow eventually approaches the Newtonian limit. However, since the value of $\mu _{\infty }$ is small for 7 % AS, this approach occurs extremely slowly. In the figure, the wall amplitude $2\epsilon$ is compensated by ${\textit{Re}}_b^{2/3}$ , as the resulting quantity is expected to vary weakly in the Newtonian case. The grey band indicates the range of ${\textit{Re}}_b$ where He et al. (Reference He, Deguchi, Song and Blackburn2025) found linear instability with $m=1$ spiral perturbations for the smooth wall case (i.e. $2\epsilon =0$ ). Within this band, the values of $2\epsilon$ at the onset of streaky instability are sufficiently small ( $2\epsilon \in [6.48\times 10^{-4},1.79\times 10^{-3}]$ ), so the spiral instability may persist. The spiral mode has a relatively long wavelength ( $L=2\pi /\alpha \in [17.167,55.116]$ ) compared with the streak mode ( $L\in [2.856,3.342]$ for $m=1$ within the band).

Figure 14.

Stability diagram using the parameters reported by Escudier et al. (Reference Escudier, Poole, Presti, Dales, Nouar, Desaubry, Graham and Pullum2005). The computational procedure is the same as in figure 13. The vertical dash-dotted line indicates the values of ${\textit{Re}}_b$ at which Escudier et al. (Reference Escudier, Poole, Presti, Dales, Nouar, Desaubry, Graham and Pullum2005) observed the asymmetric mean flow. (a) 0.2 % PAA; (b) 0.2 % XG.

Panels (a) and (b) in figure 14 show the similar stability results for 0.2 % PAA ( $\varLambda =13.05$ ) and 0.2 % XG ( $\varLambda =0.134$ ), respectively. In those cases, the spiral instability found in the smooth wall case does not occur, because $\mu _{\infty }$ and $n$ are not small enough. However, Escudier et al. (Reference Escudier, Poole, Presti, Dales, Nouar, Desaubry, Graham and Pullum2005) reported another phenomenon not seen in Newtonian fluids: the emergence of a strongly asymmetric mean flow. We examine the magnitude of wall undulation required for the instability at the parameter values used in the experiments. If the predicted value of $\epsilon$ is sufficiently large, this may rule out the possibility that the streaky instability mechanism contributes to the generation of the asymmetric mean flow.

An asymmetric mean flow was observed at ${\textit{Re}}_b=10950$ in the experiments of Escudier et al. (Reference Escudier, Poole, Presti, Dales, Nouar, Desaubry, Graham and Pullum2005) using 0.2 % PAA. According to figure 14(a), at this ${\textit{Re}}_b$ , the $m=3$ instability emerges at $2\epsilon \approx 0.0124$ with an optimal wavelength $L\approx 0.962$ . In dimensional forms, the wall amplitude and wavelength become $0.0124\times 0.05$ m $=$ 0.62 mm and $0.962 \times 0.05$ m = 4.81 cm, respectively. Since the asymmetric mean flow can arise directly from the $m=1$ mode, we performed a similar estimate and obtained 13 % larger critical wall amplitude ( $2\epsilon \approx 0.0140$ ) associated with a wavelength 2.4 times longer. Unless the pipe is intentionally manufactured with such corrugations, it is not realistic for this type of wall geometry to arise in a typical smooth experimental pipe.

For 0.2 % XG, analysed in figure 14(b), Escudier et al. (Reference Escudier, Poole, Presti, Dales, Nouar, Desaubry, Graham and Pullum2005) observed an asymmetric mean flow at a lower ${\textit{Re}}_b=1620$ . At this Reynolds number, triggering the streaky instability requires a larger wall amplitude $2\epsilon \approx 0.0718$ . In dimensional terms, this corresponds to 3.59 mm, and the $m=3$ mode emerges with a wavelength of 3.51 cm. To excite the $m=1$ mode, a wall amplitude approximately 7 % larger (3.85 mm) is required, with the optimal wavelength of 7.08 cm. Again, such large wall undulations are unlikely to exist in the experimental set-up, so we conclude that it is highly unlikely that the streaky instability contributes to the generation of the asymmetric mean flow.