2. Interaction equations for Hagen–Poiseuille flow with transpiration

We use the centreline axial velocity $U_c$ and pipe radius $a$ as velocity and length scales to define the Reynolds number $Re=U_c a/\nu$, where $\nu$ is kinematic viscosity. The incompressible Navier–Stokes equations in coordinate-free form are then

(2.1)\begin{equation} \partial_t\hat{\boldsymbol{u}} + \hat{\boldsymbol{u}}\boldsymbol{\cdot}\boldsymbol{\nabla}{\hat{\boldsymbol{u}}} ={-}\boldsymbol{\nabla} \hat{p} + {Re}^{{-}1}\nabla^2\hat{\boldsymbol{u}}, \quad \boldsymbol{\nabla}\boldsymbol{\cdot} \hat{\boldsymbol{u}} = 0, \end{equation}

which we utilize in cylindrical polar coordinates $(r,\theta,z)$ with corresponding velocity $\hat {\boldsymbol {u}} = (\hat {u},\hat {v},\hat {w})$. The transpiration at the wall takes the form $\hat {u}=2\delta \cos (\beta z)$ with axial wavenumber $\beta$, and, as shown in Hall (Reference Hall2022), for $Re\gg 1$, the flow adjusts to the transpiration in a viscous wall layer of thickness $Re^{-1/3}$. Since the layer is thin, the structure in the wall layer is essentially identical to that given by Hall (Reference Hall2022) for channel flows; therefore, we give only the most important points of the analysis in the wall layer here and refer the reader to that work for full details.

The flow in the wall layer driven by the transpiration velocity leads to velocity components in the axial and radial directions of respective sizes $\delta Re ^{1/3}$ and $\delta$, respectively. That flow then drives a flow in the core region where $1-r=O(1)$ with radial and azimuthal velocities of size $\delta$. However, we assume that the leading-order flow in the core is a steady roll–streak flow perturbed at higher order by the effect of wall transpiration. It is well known for streamwise vortex structures that the roll velocity field is smaller by a factor $Re$ than the streak. Thus, we take the velocity in the core to be

(2.2)\begin{equation} \hat{\boldsymbol{u}}=( Re^{{-}1}u(r,\theta),Re^{{-}1}v(r,\theta),w(r,\theta)) + \text{small wave}, \end{equation}

with corresponding perturbation pressure $p(r,\theta )Re^{-2}$, where the flow is $2{\rm \pi} /N$-periodic in the azimuthal direction $\theta$ and the magnitude of the small wave is $O(\delta )$, where the relationship between $\delta$ and $Re$ will be fixed in what follows by the requirement that the roll–streak flow is sustained by the wave's Reynolds stresses in the wall layer. Inserting (2.2) into (2.1) one finds that within the core the leading-order approximation to the steady roll–streak flow for $Re \gg 1$ is

(2.3a)$$\begin{gather} \frac{\partial^2 u }{\partial r^2} + \frac{1}{r^2}\frac{\partial^2 u}{\partial\theta^2} + \frac{1}{r} \frac{\partial u}{\partial r} - \frac{2}{r^2} \frac{\partial v}{\partial\theta} - \frac{u}{r^2} ={-} \frac{\partial p}{\partial r} + u\frac{\partial u} {\partial r} + \frac{v}{r}\left(\frac{\partial u}{\partial\theta}-v\right), \end{gather}$$

(2.3b)$$\begin{gather}\frac{\partial^2 v}{\partial r^2} + \frac{1}{r^2} \frac{\partial^2 v}{\partial\theta^2} + \frac{1}{r} \frac{\partial v}{\partial r} + \frac{2}{r^2} \frac{\partial u} {\partial\theta} - \frac{v}{r^2} ={-} {\frac{1}{r}}\frac{\partial p}{\partial\theta} + u\frac{\partial v}{\partial r} +\frac{v}{r}\left(\frac{\partial v}{\partial\theta}+u\right), \end{gather}$$

(2.3c)$$\begin{gather}\frac{\partial^2 w}{\partial r^2} + \frac{1}{r^2} \frac{\partial^2 w}{\partial\theta^2} + \frac{1}{r} \frac{\partial w}{\partial r} = {4}+ u\frac{\partial w}{\partial r} + \frac{v}{r}\frac{\partial w}{\partial\theta}, \end{gather}$$

(2.3d)$$\begin{gather}\frac{\partial u} {\partial r} + \frac{u}{r} + \frac{1}{r} \frac{\partial v} {\partial\theta} = 0. \end{gather}$$

Note that the Reynolds number no longer appears directly owing to the scaling introduced in (2.2). In the absence of any coupling between the roll–streak flow and the wave, the solution of (2.3) reduces to Hagen–Poiseuille flow since that is stable to streamwise vortex perturbations. The coupling takes place in the viscous wall layer where the size of the transpiration velocity is increased until Reynolds stresses associated with that wave field produce an azimuthal velocity at the edge of the wall layer comparable in size to that of the roll–streak flow in the core. Within the viscous wall layer we introduce the critical-layer scaling variable $\zeta =Re^{1/3}(1-r)$ and represent the wall shear by $\lambda (\theta )=-w_r(1,\theta )$. Additionally defining $E=\mathrm {e}^{\mathrm {i} \beta z}$ to express axial waviness, we disturb the flow driven by the wave to an axially periodic disturbance of size $\delta$ so the velocity and pressure in the wall layer may be written as

(2.4a)$$\begin{gather} \hat{u}=\frac{u_m(\zeta, \theta)}{Re^{4/3}}+\cdots + {\delta} (U(\zeta,\theta)E+\mathrm{c.c.})+\cdots, \end{gather}$$

(2.4b)$$\begin{gather}\hat{v}=\frac{v_m(\zeta,\theta)}{Re }+\cdots +\delta Re^{1/3}(V(\zeta,\theta)E+\mathrm{c.c.})+\cdots, \end{gather}$$

(2.4c)$$\begin{gather}\hat{w}=\frac{\zeta \lambda(\theta)}{Re^{1/3}}+\cdots +\delta Re^{1/3}(W(\zeta,\theta)E+\mathrm{c.c.})+\cdots, \end{gather}$$

(2.4d)$$\begin{gather}\hat{p}={-\frac{4z}{Re}}+\frac{p_m(\zeta,\theta)}{Re^2 }+\cdots + {\delta} (P(\theta)E+\mathrm{c.c.})+\cdots, \end{gather}$$

where $u_m$, $v_m$ and $p_m$ represent two-dimensional streamwise averages. Here we have assumed that the wave component of the pressure is independent of $\zeta$ (reasonable, because the wall layer is thin) so that substituting the above expansions into the equations of motion we find that the $(\theta,z)$ momentum equations and continuity reduce to

(2.5a)$$\begin{gather} -U_\zeta+V_\theta+\mathrm{i} \beta W = 0, \end{gather}$$

(2.5b)$$\begin{gather}\mathrm{i}\beta \lambda \zeta V ={-}P^{\prime} +V_{\zeta\zeta}, \end{gather}$$

(2.5c)$$\begin{gather}\mathrm{i} \beta \lambda \zeta W -\lambda U + V\zeta \lambda^{\prime} ={-}\mathrm{i}\beta P +W_{\zeta\zeta}. \end{gather}$$

The above equations are to be solved subject to

(2.6a)$$\begin{gather} (U,V,W) = {(1,0,0)}\quad \text{at}\ \zeta=0, \end{gather}$$

(2.6b)$$\begin{gather}(U_\zeta, V, W)\rightarrow 0 \quad \text{as } \zeta \rightarrow \infty, \end{gather}$$

where (2.6a) represents the wall transpiration. Following Hall (Reference Hall2022) we find that

(2.7a)$$\begin{gather} U= 1-\frac{\lambda^{\prime}(\theta) P^{\prime}(\theta)}{2 \lambda(\theta) \gamma}[{\mathcal{L}}^{\prime\prime}(\xi)-1], \end{gather}$$

(2.7b)$$\begin{gather}V= \gamma^{-{2/3}} P^{\prime}(\theta){\mathcal{L}}(\xi), \end{gather}$$

where $\gamma =\mathrm {i} \beta \lambda (\theta )$, $\xi =\gamma ^{1/3}\zeta$, a dash denotes an ordinary derivative with respect to argument, and ${\mathcal {L}}$ satisfies

(2.8)\begin{equation} {\mathcal{L}}^{\prime\prime}-\xi{\mathcal{L}}=1, \quad{\mathcal{L}}=0\text{ at } \xi=0 \text{ and } \infty. \end{equation}

The wave component of pressure is then given by

(2.9)\begin{equation} P^{\prime\prime} - \frac{3\lambda^\prime}{2\lambda} P^{\prime}- \beta^2 P ={-}\mathrm{i}\beta\lambda. \end{equation}

In the absence of forcing from the wall the right-hand side of the above equation is zero and the equation reduces to that given in Hall & Smith (Reference Hall and Smith1991).

Now consider the $z$-independent part of the velocity field in the wall layer. In particular taking the mean in $z$ of the azimuthal momentum equation there we find that $v_m$ satisfies

(2.10)\begin{equation} \frac{\partial^2 v_m}{\partial \zeta^2}= \delta^2Re^2[-(UV^*)_\zeta+2VV_\theta^*]+\mathrm{c.c.}, \end{equation}

where $^\ast$ denotes complex conjugation. For large $\zeta$, the term in the square brackets is proportional to $\zeta ^{-2}$ so that $v_m\sim \log \zeta$ for $\zeta \gg 1$. Noting that $\zeta =Re^{1/3} (1-r)$ we see that the $z$-independent azimuthal flow driven by the Reynolds stress terms in the wall layer reinforces the azimuthal velocity component of the roll–streak flow in the core if ${\delta ^2 } Re^2\log Re =O(1)$. Thus, we introduce a scaled transpiration parameter $\varLambda$ by writing

(2.11)\begin{equation} \delta=\varLambda Re^{{-}1}[\log Re]^{{-}1/2} \end{equation}

and using the large-$\zeta$ form of the right-hand side of (2.10), we find that the matching condition with the core flow then becomes

(2.12)\begin{equation} (u,v,w)=\left[0,-\frac{\varLambda^2}{3 \lambda(\theta)^2} ( \vert P(\theta) \vert^2+\beta^{{-}2} \vert P^{\prime}(\theta) \vert^2 )^{\prime}, 0\right]. \end{equation}

The right-hand side of (2.9) is due to the transpiration forcing and setting it equal to zero gives the free-interaction problem VWIb. On the other hand, if we drop the $\theta$-dependence, the pressure equation gives the pressure field driven entirely by the transpiration. But then, the roll–streak flow is no longer driven by the wave and no interaction occurs. Thus, the roll–streak equations coupled with the boundary conditions and pressure allow for an exact coherent structure in pipe flow in the presence of transpiration. However, there is a significant difference between the VWIa and transpiration-augmented VWIb interaction problems. The difference is that in the latter problem nonlinear structures bifurcate from the transpiration-modified basic Hagen–Poiseuille flow whereas in VWIa nonlinear states arise from a saddle-node bifurcation.

The linear stability of transpiration-modified basic Hagen–Poiseuille flow is obtained from the interaction equations (2.3), (2.9), (2.12) by perturbing the flow to a small amplitude roll–streak flow in the manner given in Hall (Reference Hall2022). In the core region we write

(2.13)$$\begin{gather} {\boldsymbol{u}}=(0,0,1-r^2)+A(u_1(r)\cos N\theta, v_1(r) \sin N\theta, w_1(r)\cos N\theta)+\cdots, \end{gather}$$

(2.14)$$\begin{gather}p=Ap_1(r)\cos N\theta+\cdots, \end{gather}$$

where $N$ is an integer and $A$ is small. Linearizing (2.3) we obtain a set of linear equations for ${\boldsymbol {u}}_1,p_1$ which we solve subject to regularity at $r=0$ and $u_1,w_1$ vanishing at $r=1$. We obtain

(2.15)\begin{align} {\boldsymbol{u}} _1(r) &= \left(N[r^{N-1}-r^{N+1}],[N+2]r^{N+1}-Nr^{N-1}, \vphantom{\left[\frac{r^2+1}{N+2}\right]}\right.\nonumber\\ &\qquad \left.\frac{Nr^N(r^2-1)}{4}\left[\frac{r^2+1}{N+2}-\frac{2}{N+1}\right]\right), \end{align}

which gives

(2.16)\begin{equation} \frac{v_1(1)}{w'_1(1)}={-}2\frac{(N+1)(N+2)}{N}. \end{equation}

Next we find the small amplitude solution of the pressure equation (2.9) by first writing

(2.17)\begin{equation} P=P_0+AP_1\cos N\theta+\cdots, \end{equation}

and then equating terms $O(1), O(A)$ in turn we obtain

(2.18a,b)\begin{equation} P_0={-}\frac{2\mathrm{i}}{\beta},\quad P_1=\frac{\mathrm{i}\beta w_1'(1)}{\beta^2+N^2}. \end{equation}

We then substitute the above pressure expansion into (2.12) and equate terms $O(A)$ to give

(2.19)\begin{equation} \frac{v_1(1)}{w'_1(1)}={-}\frac{\varLambda^2N }{3[\beta^2+N^2]}, \end{equation}

so comparing the above two expressions for $v_1(1)/w_1'(1)$ we obtain

(2.20)\begin{equation} \varLambda^2=6(1+\beta^2/N^2)(N+1)(N+2). \end{equation}

Figure 1 shows the neutral curves for various wavenumbers $\beta$ where $\varLambda$ is shown as a function of $N$ in figure 1(a). For a given value of $\beta$, the value of $N$ that minimises $\varLambda$ satisfies $\beta ^2=N^3(3 + 2 N)/ (4 + 3 N)$. Since $N$ must be an integer, the most dangerous mode will correspond to $N=1$ when $\beta < 2^{1/2}$ and to an integer value number $N\approx [(3/2)\beta ^2]^{1/3}$ for large $N$. Figure 1(b) displays $\varLambda _c=\min _{N}\varLambda$ as a function of $\beta$. Numbers above solid/dashed curves indicate the value of $N$ that minimizes $\varLambda$ for a given $\beta$ which corresponds to the most dangerous mode.

Figure 1.

(a) Plot of $\varLambda$ vs $N$ for given $\beta$ values. Black markers show the minimum $\varLambda$ values to indicate the most dangerous modes $N=1$, 5, 15 for $\beta =1$, 10, 50, respectively. (b) Plot of $\varLambda _c$ vs $\beta$ where numbers over solid/dashed curves indicate the most dangerous mode. (c) Critical $\delta$ value as a function of $Re$ for $N=3$ (red) and $N=4$ (blue) in comparison to the expected asymptotic scaling given by (2.11) and (2.20) (for $N=3$, $\beta =10$) shown by dashed line.

Figure 1(c) displays independent calculations of wall transpiration amplitudes that lead to neutral stability of the base flow as a function of $Re$. These calculations were performed at $\beta =10$ for disturbances with $N=3$ and $N=4$ via full linearized Navier–Stokes calculations; see Blackburn (Reference Blackburn2002) and Barkley, Blackburn & Sherwin (Reference Barkley, Blackburn and Sherwin2008). We see that the asymptotic and numerical solutions are in progressively better agreement as $Re$ increases but it is interesting to note that the agreement is not as good as found in the planar case by Hall (Reference Hall2022). We expect that agreement would be improved by continuing the expansion to higher order but our emphasis here is on the nonlinear problem and so we do not do so.

3. The numerical solution of the full interaction equations

Our computational approach is based on the numerical formulation described in Ozcakir et al. (Reference Ozcakir, Tanveer, Hall and Overman2016), which uses a Galerkin truncation of Fourier modes in $\theta$ and a Chebyshev representation in $r$ with appropriately chosen velocity functions $\varPhi _j$, $\varPsi _j$ satisfying no-slip boundary conditions. Here, we give a brief discussion of the key elements that are added to the existing solver for the present problem. Within the core, the nonlinear equations (2.3) describing the roll–streak flow are solved numerically. It is convenient to decompose the velocity field $(u,v,w)$ as

(3.1)\begin{equation} (u,v,w)={\boldsymbol v}_B+{\boldsymbol V}(r,\theta), \end{equation}

where ${\boldsymbol v}_B(r,\theta )$ is a solenoidal base velocity field satisfying the boundary conditions (2.12) which includes Hagen–Poiseuille flow, ${\boldsymbol V}(r,\theta )$ is the perturbation velocity with no-slip boundary conditions. Using (3.1) in equations (2.3), one obtains

(3.2)\begin{equation} {\boldsymbol V} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol V}+ {\boldsymbol v}_B \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol V}+ {\boldsymbol V}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol v_B + \boldsymbol{\nabla} p-\nabla^2 {\boldsymbol V} = \nabla^2 {\boldsymbol v}_B - {\boldsymbol v}_B \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol v}_B, \end{equation}

where $p(r,\theta )$ is the perturbation pressure. This formulation allows us to find $N$-fold rotationally symmetric ($\mathbb {R}_{N}$) solutions ${\boldsymbol V}(r,\theta )$ using the truncated basis representations

(3.3)\begin{equation} \boldsymbol{V}\approx\sum[{u}_{jk} \varPhi_j ({ r};k N) \cos{k N \theta}, {v}_{jk} \varPhi_j ({ r};k N) \sin{k N \theta}, {w}_{jk} \varPsi_j ({r};k N) \cos{k N \theta}] \end{equation}

for the radial, azimuthal and axial perturbation velocity components, where $\varPhi _j$, $\varPsi _j$ are given in terms of Chebyshev polynomials $T_{j}$ as

(3.4a)\begin{equation} \varPhi_j (r; k) =\begin{cases} T_{2j+2} (r) - T_{2j} (r) & \text{if } k \text{ is odd},\\ T_{2j+3} (r) - T_{2j+1} (r) & \text{if } k \text{ is even}, \end{cases} \end{equation}

and

(3.4b)\begin{equation} \varPsi_j (r; k) =\begin{cases} T_{2j+3} (r) - T_{2j+1} (r) & \text{if } k \text{ is odd},\\ T_{2j+2} (r) - T_{2j} (r) & \text{if } k \text{ is even}. \end{cases} \end{equation}

It should be noted that this representation of $\varPhi _j$ and $\varPsi _j$ forms a complete radial basis for the boundary value problem. The sums are taken over $j\in [0, J]$ and $k\in [0, M]$ for large enough integer values of $J$ and $M$ to obtain three-digit relative accuracy for each velocity component in the solution. The cos/sin basis functions are chosen to fix the azimuthal phase of the flow structure. Note that the relative phase in $\theta$ between the radial and azimuthal velocity components is fixed by continuity. The relative phase of the streamwise velocity component is then determined through the ‘lift-up’ effect term $u \partial w/\partial r$ in (2.3c). Inserting this basis with given $N$, $\beta$ and $\varLambda$ into (3.2) produces a nonlinear algebraic system for $\{u_{jk},v_{jk},w_{jk} \}$, which is solved via Newton's method and the generalized minimal residual method.

For given $N,\beta$, we determine $\varLambda _0$ at which the flow becomes unstable using the linear stability results of § 2. A small amplitude neutral eigenfunction is then used as an initial guess in the Newton iteration procedure near this neutrally stable point $\varLambda =\varLambda _0$. At each Newton step, the linear ordinary differential equation (2.9) is solved using spectral methods for coefficients $\{p_k,1< k< M\}$, where $P(\theta )=\sum p_k \sin {k N \theta }$. The resulting $P(\theta )$ is then used to generate boundary conditions (2.12) which in turn are used to construct the base flow ${\boldsymbol v}_B(r,\theta )$. Once a converged solution is found, an arc-length continuation procedure with continuation parameter $\varLambda$ is used to determine other solutions. The corresponding eigenvalue problems were solved numerically using an Arnoldi method.