# A linear system for pipe flow stability analysis allowing for boundary condition modifications

An accurate system to study the stability of pipe flow that ensures regularity is presented. The system produces a spectrum that is as accurate as Meseguer & Trefethen (2000), while providing flexibility to amend the boundary conditions without a need to modify the formulation. The accuracy is achieved by formulating the state variables to behave as analytic functions. We show that the resulting system retains the regular singularity at the pipe centre with a multiplicity of poles such that the wall boundary conditions are complemented with precisely the needed number of regularity conditions for obtaining unique solutions. In the case of axisymmetric and axially constant perturbations the computed eigenvalues match, to double precision accuracy, the values predicted by the analytical characteristic relations. The derived system is used to obtain the optimal inviscid disturbance pattern, which is found to hold similar structure as in plane shear flows.

## Authors

• 2 publications
• 1 publication
01/07/2022

### Wavenumber-explicit hp-FEM analysis for Maxwell's equations with impedance boundary conditions

The time-harmonic Maxwell equations at high wavenumber k in domains with...
10/15/2019

### Numerical multiscale methods and effective boundary conditions

We develop numerical multiscale methods for viscous boundary layer flow....
06/08/2019

### Modified symmetry technique for mitigation of flow leak near corners for compressible inviscid fluid flow

Using the standard symmetry technique for applying boundary conditions f...
04/27/2021

### An optimal control approach to determine resistance-type boundary conditions from in-vivo data for cardiovascular simulations

The choice of appropriate boundary conditions is a fundamental step in c...
11/01/2019

### Cut Bogner-Fox-Schmit Elements for Plates

We present and analyze a method for thin plates based on cut Bogner-Fox-...
07/02/2020

### Weak Boundary Condition Enforcement for Linear Kirchhoff-Love Shells: Formulation, Error Analysis, and Verification

Stable and accurate modeling of thin shells requires proper enforcement ...
06/21/2021

### The implementation of a broad class of boundary conditions for non-linear hyperbolic systems

We propose methods that augment existing numerical schemes for the simul...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Flow through pipes is of a great importance due to vast applications ranging from cooling systems to fluid transportations©2019. This manuscript version is made available under the CC-BY-NC-ND 4.0 license http://creativecommons.org/licenses/by-nc-nd/4.0/. As drag reduction in such systems will enhance human life by lowering energy consumption, understanding laminar to turbulent transition has been the subject of research since Reynolds’ experiments. Despite that the flow exhibits instability in reality, it is linearly stable mathematically. This phenomenon has induced research to tackle the understanding in the perspectives of nonlinear dynamics, chaos and intermediate edge states itano2001dynamics; schneider2007turbulence; budanur2018complexity, relatively recently.

However, for the nonlinearity to set in, the disturbances need to be brought to a certain finite-scale size. Algebraic transient growth boberg1988onset; butler1992three; trefethen1993hydrodynamic; schmid1994optimal; schmid2007nonmodal caused by non-normality of the linearized Navier-Stokes operator can mathematically explain this initial growth. In physical terms, the transient growth is described by the inviscid interaction of the mean flow and perturbation through the lift-up effect for the modes with azimuthal modulations and through the Orr mechanism of enhancing the perturbations that oppose the mean shear in the case of modes with streamwise modulations. Inherently, since the laminar mean flow exhibits a mean shear, the flow has a tendency to transfer energy into infinitesimal disturbances, which can be inflow borne or originate from wall roughness and vibrations.

One of the other several routes to turbulence is distortions of the mean flow. Spatial analysis of this flow configuration predicted instability under suitable mean flow distortion gavarini2004initial. Furthermore, irregularities in the pipe geometry can serve as a cause for instability. For example, the flow has been found unstable in slowly diverging pipes sahu2005stability. In addition, flow through sinusoidally corrugated pipe has exhibited linear instability loh2011stability. Wall roughness in the molecular scales has been found to be a destabilizing factor compared to smooth pipe, though instability has not been observed pruuvsa2009influence.

Among the early investigations that predicted the linear stability of flow through straight pipes in temporal setting lessen1968stability; burridge1969comments; salwen1972stability; salwen1980linear, Burridge & Drazin burridge1969comments

proposed an ansatz for their working variables for an asymptotic analysis, which was later adopted for numerical calculations

schmid1994optimal; o1994transient; thomas2012linear.

Priymak & Miyazaki priymak1998accurate predicted the form of the , the radial coordinate, dependence of the velocity and pressure variables through an ansatz, which is suitable if solutions of the linear system for the perturbations that are analytic at the centreline are sought. It should be noted that this ansatz already incorporates the necessary behaviour for velocity fields deduced by Khorrami et al. khorrami1989application, thus leaving only regularity as means for further conditioning of the system. We would like to note that in the light of this ansatz, the working variables of Burridge & Drazin burridge1969comments, namely, and (defined in section 3.2), which are related to the radial velocity and radial vorticity have the property, for to the leading order close to the centreline, and for , where is the azimuthal wavenumber. This suggests that for higher values of (e.g.,

), all the eigenfunctions will be vanishingly small in a finite neighbourhood of the centreline. Consequently, there will be precision loss in determining the eigenfunctions. This complication arises because the eigenvalues are multiplied by exceedingly low values in the discretized eigensystem equations close to the the centreline.

Furthermore, an increase of collocation points would only worsen the situation by enhancing the round-off errors since the eigenfunctions do not undergo steep variations close to the centreline. These round-off errors at the end of the iterative procedures for eigenvalues demonstrate themselves as pseudospectra. Since pseudospectra generally affects modes decaying at higher rate than the least decaying ones, the optimal patterns can get distorted at such high values of .

Priymak & Miyazaki priymak1998accurate also developed a method for solving the linearized Navier-Stokes equations cast as a system of two unknowns, the radial and azimuthal perturbation velocities, that utilizes the above findings, thus producing very accurate spectra. However, the linear operator itself needs to be determined in a complex algorithmic way, unlike the case of Burridge & Drazin burridge1969comments. In priymak1998accurate, the linear system is determined through a sequence of steps numerically, where the pressure needs to be found by inverting a Poisson operator, and the continuity is being imposed by a correcting pressure which leads to correcting velocities, similar to the steps of the SIMPLE algorithm.

Meseguer & Trefethen meseguer2000spectral; meseguer2003linearized deployed the ansatz of Priymak and Miyazaki, and have produced very accurate results with a number of collocation points as low as . Their work is in Petrov-Galerkin formalism which uses solenoidal trial basis and test functions that satisfies the particular boundary condition of no-slip. However, this method comes with a task of changing the trial and test bases if the boundary conditions are other than no-slip. In some occasions, the boundary conditions are even time dependent (see for example, malik2018growth). Finding appropriate solenoidal trial and test bases that satisfy arbitrary boundary conditions and the properties prescribed by the ansatz of Priymak & Miyazaki is not a simple task, in general.

The above facts motivate us to formulate a method that is as accurate and efficient as Priymak & Miyazaki priymak1998accurate and Meseguer & Trefethen meseguer2000spectral, but as flexible and simple as a usual spectral collocation method with 2-tuple state variable similar to that of Ref. schmid1994optimal or burridge1969comments. To this end, we use the ansatz of Priymak & Miyazaki for velocity fields to determine a 2-tuple working variables that do not go like a power law close to the centreline. As the boundary conditions will be imposed on the unknowns unlike the Petrov-Galerkin method, the formulated equations are applicable for a range of boundary conditions.

The theoretical derivations are presented in section 2, while the numerical results are discussed in section 3. Finally, in section 4, we derive the Ellingsen and Palm solutions ellingsen1975stability that are valid for inviscid axially constant modes that captures nonmodal algebraic growth. Derivation of this is remarkably simpler in the working variables in this paper. Such solutions have been found useful in the case of plane-shear flows to compute the optimal patterns that demonstrate the lift-up effect leading to streaks. We find that the pipe flow configuration retains these features. Although this is known through earlier viscous computations schmid1994optimal and DNS zikanov1996instability, we arrive at these results in the inviscid limit.

## 2 Theoretical derivations

Let us consider the cylindrical coordinates, , which are the axial, radial and azimuthal coordinates, respectively. We emphasize the unusual order of the coordinates used here. Let be the velocities normalized with respect to the centreline velocity. Let us use the decomposition, where the mean flow, and , and the considered perturbation, is such that , which are governed by,

 ~ut+(U⋅∇)~u+(~u⋅∇)U=−∇~p+Re−1∇2~u, (1)

and by the continuity condition, . In Eq. (1), , where is the dimensional centreline velocity, is the pipe radius, and is the kinematic viscosity. The in Eq. (1) is the perturbation pressure. Substituting into Eq. (1) to consider a single Fourier mode, and upon taking the curl to remove the appearance of pressure, we obtain the following three equations:

 i(αU−ω)(Du′−iαv′)=−iαUru′−(Urr+UrD)v′+Re−1(Δ(Du′−iαv′)−2inr−2η′) (2)
 i(αU−ω)[D(rw′)−inv′]=−iαrUrw′+Re−1(¯¯¯¯¯Δ[D(rw′)−inv′]) (3)
 i(αU−ω)η′=−inUrr−1v′+Re−1(Δη′+2nr−2(αv′+iDu′)) (4)

where , , and . As a matter of convention, we use suffix to imply differentiation of mean flow variables with respect to , and the operator to represent the same for the fluctuation variables. Using the above definition of and the continuity equation, the variables and can be written as

 u′= ird−1(αv′+αrDv′−nη′)and (5) w′= id−1(nv′+nrDv′+αr2η′). (6)

Priymak & Miyazaki priymak1998accurate observed the following behaviour close to the centreline.

 (7)

where and, , and

are analytic functions having Taylor expansion around the centreline with vanishing coefficients of odd powers of

. Substituting Eq. (7) into the definition of , we get its behaviour. In summary, and have the following forms:

 (v′,η′)={(rϕ,rΩ)% forn=0and(rℓϕ,rℓΩ)forn≠0 (8)

with an analytical function having similar Taylor expansion as . Substituting Eq. (8) into Eq. (5) and Eq. (6), we get

 u′= ⎧⎨⎩irℓ+1d[α(ℓ+1)ϕ+αrDϕ−nΩ]n≠0iα(2ϕ+rDϕ)n=0 (9) w′= {irℓd−1[n(ℓ+1)ϕ+nrDϕ+αr2Ω]n≠0,iα−1rΩn=0. (10)

First, let us consider the case of . Upon substituting Eqs. (8)-(10) into Eqs. (2)-(4) we get

 (ω−αU)(αL1ϕ−nL2Ω)=αL3ϕ−αnL4Ω+iRe−1(αL5ϕ−nL6Ω), (11)
 (ω−αU)(nL7ϕ+αL8Ω)=nL9ϕ−α2L10Ω+iRe−1(nL11ϕ+αL12Ω),and (12)
 (ω−αU)Ω=nr−1Urϕ+iRe−1(2αnL13ϕ+L14Ω) (13)

where the operators, are given in A. Eq. (11) and Eq. (12) are fourth order in and third order in . As we have three equations, Eqs. (11)-(13) for two unknowns, and , we use Eq. (11) and Eq. (12) to reduce the order of . The resulting equation that is fourth order in and second order in is

 (ω−αU)(rdL1ϕ+2αnr2d−1Ω)=αr(Ur−rUrr)ϕ+iRe−1(L15ϕ+4αnL16Ω), (14)

where and are operators given in A.

Four more tasks are carried out before arriving at the required final system to work with: Firstly, the appearance of second derivative of in Eq. (14) is removed with the help of Eq. (13). Though this will not reduce the order of the system, it helps in deriving simpler conditions of regularity which are shown later. The outcome is that the order of the derivatives of in the yet-to-be-derived regularity condition is less than the order of the governing equation.

Secondly, the even parity nature of and , which have Taylor series expansions and is implemented by the transformation, . Since this will transform and to general analytic functions, Chebyshev polynomials of both odd and even orders can be used for the spectral expansion without any loss of efficiency. (For another method, where only even Chebyshev polynomials are used, see Priymak & Miyazaki priymak1998accurate and Meseguer & Trefethen meseguer2000spectral).

Thirdly, similar operations performed for the case of by substituting and from Eq. (8), give the required system for axisymmetric disturbances. However, it should be noted that this system is derived without using Eq. (3) as there will not be a need to reduce the order of in the system. The obtained system from Eq. (2) and Eq. (4) will already be fully decoupled in this case with six as the order of the system.

Finally, a system suitable for the case of , the axisymmetric and axially constant modes, is derived with as the independent variable. This special case requires a different set of working variables, as (hence, ) for all due to continuity, and (hence, ) by definition. We choose the variables, and defined in Eq. (7) for this purpose. The required system for this case can be derived from the axial and azimuthal components of Eq. (1).

### 2.1 Linear System

In summary, the final system of equations are

 (ω−αU)([−α2g1+g2D+4yd3D2]ϕ−2αnd2Ω)=−8αd2n2Uyϕ+iRe−1([α4g3−8α2g4D+g5D2+8y(g2+4d3)D3+16y2d3D4]ϕ−8α3n[d2+2dyD]Ω), (15)
 (ω−αU)d2Ω=2nd2Uyϕ+iRe−1(2αn{α2g6−2[d2+(ℓ+3)d]D−4ydD2}ϕ+{−α2g7+4d[(ℓ+1)d+n2]D+4yd2D2}Ω), (16)

for the case of , and

 (ω−αU)(8D+4yD2−α2)ϕ=iRe−1(16y2D4+96yD3+(96−8α2y)D2−16α2D+α4)ϕ (17)
 (ω−αU)Ω=iRe−1(4yD2+8D−α2)Ω (18)

for the case of , and

 −iωψ1= 4Re−1(yD2+D)ψ1 (19) −iωψ2= 4Re−1(yD2+2D)ψ2 (20)

for the case of , where , and the and () are functions as given in A. The order of the system is six for each of the cases of and , whereas it is four in the case of . In the latter case, the order is only four owing to the fact that we did not have to take the curl of Eq. (1) to remove the appearance of , as pressure is a constant in the axial and azimuthal directions similar to the other perturbation quantities. All modes of this flow configuration, for all cases of and , are observed to be decaying. For the special case of , a physical reason can be deduced since the time evolution is dictated only by the viscous dissipation because both the mean shear and the pressure gradient terms are zero.

### 2.2 Boundary and Regularity Conditions

The systems given by Eqs. (15)-(16) and Eqs. (17)-(18) need to be solved with six conditions for a unique solution, and four conditions for the system of Eqs. (19)-(20). In the case of and , three of the needed six conditions are straight forward: they are the no-slip and no-penetration conditions at the pipe boundary. These conditions in terms of the variables and read as

 ϕ=Dϕ=Ω=0aty=1. (21)

The other three conditions should come from the regularity of solutions at . Note that, although we are seeking analytic solutions by making the required transformations, with relations given in Eq. (8), the appearance of regular singularity in Eq. (15)-(18) implies that they still can allow non-analytic solutions such as, for analogy, the Bessel functions of second kind for positive integer in the case of the Bessel equation. In a numerical procedure, the transformation in Eq. (8) and the transformation of the coordinate, would only guarantee the accuracy by shifting the focus to the factors of and , namely, and , whose variations are unknown, and allow us to work with a reduced number of Chebyshev polynomials in the spectral expansion, while still allowing non-analytic solutions. Therefore, we need regularity conditions to explicitly rule out non-analytic solutions in our solution procedure.

Since we seek and as analytic functions, they have Taylor series expansion around . This warrants all orders of the derivatives to be of , which implies that all derivative terms in Eq. (15)-(18) that are multiplied by should vanish in the limit . For analogy, such requirement is satisfied by the the analytic , the Bessel functions of first kind for positive integer , in the case of the Bessel equation. This instantly gives us two conditions of regularity for each of the cases, which are

 (ω−αU)([−α2g1+g2D]ϕ−2αnd2Ω)=−8αn2d2Uyϕ+iRe−1([α4g3−8α2g4D+g5D2]ϕ−8α3nd2Ω), (22)
 (ω−αU)d2Ω=2nd2Uyϕ+iRe−1(2αn{α2g6−2[d2+(ℓ+3)d]D}ϕ+{−α2g7+4d[(ℓ+1)d+n2]D}Ω), (23)

for the case of , and

 (ω−αU)(8D−α2)ϕ=iRe−1×(96D2−16α2D+α4)ϕ (24)
 (ω−αU)Ω=iRe−1(8D−α2)Ω (25)

for the case of . We need one more condition for each of these cases of and . Note that in Eq. (15) and (17) the term is being multiplied by . Requiring would only yield the constraint that . Therefore, the other condition is the one that would enforce , since is analytic. The following two facts should be noted in order to derive the required condition: (1) As the Eq. (15) and (17) are both valid for the whole domain, their derivatives with respect to are also valid through out the flow field; (2) Due to regular singularity, Eq. (22) and (24) are at most second order in and zero-th order in , allowing them to be treated as boundary conditions. Furthermore, the derivatives of Eq. (15) and (17) with respect to at the centreline can be used as boundary conditions since they are of at most first order in and third order in , which are less by one order comparing to the orders of these variables in the governing equations. Hence, the additional regularity condition that guarantees the differentiability of until fourth order for each cases of and are the following:

 (ω−α)([−α4g8+α2g9D+g10D2]ϕ−2αn3[2α2+n2D]Ω)=[α3g11−4αn6ℓD]ϕ+2α2n5(1−D)Ω+iRe−1([α6g12+α4g13D−α2g14D2+g15D3]ϕ−8α3n[α2(ℓ−1)+n2(ℓ+3)D]Ω) (26)

for , and

 (ω−α)(12D2−α2D)ϕ=α(α2−8D)ϕ+iRe−1(192D3−24α2D2+α4D)ϕ (27)

for the case of , where are constants defined in A. In Eq. (26), we have substituted the values of mean flow variables for brevity.

Note that the multiplicity of the pole in the governing equation is such that it yields exactly three conditions to complement the three boundary conditions at the wall. Had we instead derived another regularity condition by differentiating the governing equations once more and evaluated at the pipe centre, the condition would not have served as a boundary condition since it would have had the same orders of derivatives of the unknowns as they appear in the governing equations. Similarly, had the multiplicity been lower than in the present case, we would not have been able to obtain a total of six conditions to solve the sixth-order system.

In the case of , the boundary and regularity conditions are given by

 ψ1=ψ2=0aty=1,and, (28)
 (ω−4iRe−1D)ψ1=(ω−8iRe−1D)ψ2=0. (29)

at .The regularity condition for in Eq. (29) and the governing equation Eq. (20) can also be obtained from the limit in Eq. (25) and Eq. (18), respectively, since for . Therefore, one can expect that a part of the spectrum for the case of can be obtained from the system for the case of in the limit of . However, there are a new set of modes for this case that cannot be obtained from the system for finite in the limit . These modes originate from Eq. (19). Both sets of these modes are governed by Sturm and Liouville theory (see for example boyce2009elementary) since the underlying operators are self-adjoint. In the next subsection, the solutions for this case are derived in coordinate so that the characteristic relation are presented in a form that is less susceptible to numerical errors. This will facilitate to validate the eigenvalues obtained from the numerical spectral collocation solutions of Eq. (19) and (20).

It should be highlighted that the boundary conditions can be changed depending on the flow configurations, without affecting the formulation, i.e., Eq. (15)-(20). This is the advantage of the present system compared with Meseguer and Trefethen meseguer2000spectral, where one would be tasked with finding appropriate basis and test functions that satisfies the boundary and regularity conditions, besides satisfying the other condition of continuity.

### 2.3 Stokes Modes with no Transient Growth

The operators on the RHS in Eq. (19)- (20) are examples of Stokes operators. The general Stokes operators are defined as Laplacians acted upon by Leray projectors which ensures continuity foias2001navier. In the present case of , the continuity is automatically satisfied without any constraint as , and hence . This implies that the Leray projector is an identity operator in the present case. It is well known that the eigenvalues of a general Stokes operators are decaying and that the operator itself is self-adjoint for no-slip boundary conditions. For this particular case of , the solutions and characteristic relation are known to be Bessel functions of and their roots at  salwen1972stability; rummler1997eigenfunctions. (For such modes of plane Poiseuille flow, see rummler1997eigenfunctionsii. For application of these Stokes modes, see batcho1994generalized.)

Eq. (19)- (20) can be rewritten as

 λψ1= (yD2+D)ψ1, (30) λψ2= (yD2+2D)ψ2, (31)

where . The differential operators on the RHS are self-adjoint under the inner product defined as the cross-sectional area integral of velocity fields, i.e., or . The Sturm and Liouville theory suggests that the eigenvalues, of Eq. (30) and of Eq. (31) with are real and the eigenfunctions, and are orthogonal in the sense,

 (32)

where is the Kronecker delta and, and are constants.

As we are interested in solutions of Eq. (30)-(31) that are analytic at the centreline, we can bypass the method of Frobenius and resort to a simple power series method. Let and , where the constants, and , are to be determined. Substituting these into Eq. (30)-(31), and matching the coefficients of like powers of , we get the eigenfunctions, and as,

 ψ1,l(y)= ¯¯¯a0(1+∑k≥1(λ1,ly)k(k!)−2)and (33) ψ2,l(y)= ¯¯b0(1+∑k≥1(λ2,ly)k[k!(k+1)!]−1), (34)

where with are the real solutions of the characteristic equation,

 1+∑k≥1λk(k!)−2=0, (35)

and are the real solutions of

 1+∑k≥1λk[k!(k+1)!]−1=0, (36)

which arise due to the no-slip conditions given by Eq. (28). The Eqs. (35) and (36) can also be written as the roots of and  salwen1972stability; rummler1997eigenfunctions, however the above expressions are helpful to obtain the eigenvalues accurately as the roots of polynomials that approximate the LHS of Eqs. (35) and (36) by cutting off the summation at some . The constants, and in Eqs. (33) and (34), respectively, can be of any value as long as the linearity of the perturbations with respect to the mean flow is respected. Without loss of generality they can be the normalization constants, and , where the inner-products are defined as and . It should be mentioned that is in fact an area element, . Under such normalizations, the constants and of Eq. (32) become, . The convergence of the series in Eq. (33) and (34) is obvious by comparison tests with the series for and , respectively. The implication of the orthogonality is that superposition of such modes cannot have transient growth. A velocity state can be represented using the notation introduced at the beginning of this section as . Then the superposition velocity is given by,

 ~u(t,y)=exp(i[αx+nθ])∑k,laklexp(Λklt)u′kl, (37)

where and are diagonal matrices for each chosen and , where in turn and are the real solutions of Eq. (35) and (36), respectively, and and are some constants.

Let us define the energy as, .

 E(t)=2−1∑k,l∑i,jaHijakl×exp([Λkl+Λij]t)∫10u′ijHu′kldy, (38)

where the superscript, ‘H’ denotes Hermitian transpose. Using the orthogonality relations given by Eq. (32) and normalizations of and , we get,

 E(t)=12(∑k|ak|2exp(2λkt)+∑l|al|2exp(2λlt)). (39)

As said in the subsection 2.1, in this present case of , the evolution of perturbations are determined by a dissipative phenomena. Therefore, and are all negative. This shows that the energy, in Eq. (39), is a monotonously decaying function of time. Since the interaction with the mean flow is null for these modes, they should be supplied with energy by other means of receptivity such as vibration of the wall or by inflow-borne disturbances. The role of acoustic disturbances may also be ruled out since extremely longwaves imply extremely low (near-zero) infrasonic frequencies of the source (since, frequency /wavelength).

## 3 Numerical Results

### 3.1 Numerical method

The system of Eqs. (15)-(18) can be written as,

 Re−1Aq=ωBq, (40)

where and, and are matrices of operators. The elements of and can be figured out from those equations (Eq. (15)-(18)) for each case of and . For discretization, we use spectral collocation method in transform space (see, for example, Ref. trefethen2000spectral or Appendix A.5 of Ref. schmid2001stability and S. C. Reddy’s codes therein). The generalized eigensystem, Eq. (40) is discretized at the points,

 yj=2−1(1+ξj),if α≤3 and Re≤6000, (41a) else yj=(exp(a)−1)−1(exp[a(1+ξj)/2]−1), (41b)

where the Gauss Lobatto points, and is the highest order of Chebyshev polynomials. Eq. (41b) represents a stretching function that maps on such that the grid points cluster towards for the parameter . For , gives best result, whereas works well for .

The stretching is needed in our present formulation for large values of and , due to the following. As the algebraic variations of the velocity fields in a power-law fashion (i.e., variations) has been factored out from the unknowns by the Priymak and Miyazaki ansatz given in Eq. (7), the functions and do not have any other constraints apart from those of satisfying regularity. This allows and to have steep variations at the pipe centre, given that there is a regular singularity at that location. The freedom of and and their derivatives from the requirement to vanish at the pipe centre, coupled with their vanishing at the wall and the presence of regular singularity there, causes a boundary layer behaviour at pipe centre. Hence, the stretching introduced in Eq. (41b) is required in our formalism. However, it is precisely this boundary layer behaviour at the pipe centre that allows for avoiding the pseudospectra, as will be shown further down.

Note that the requirement for stretching does not arise in the formulation of Burridge & Drazin, where the functions and one of their derivatives vanish at both boundaries, namely, wall and pipe centre. Therefore, upon solving their system, the variations of the solutions are spread over the entire domain, avoiding steep variations. Such clustering is neither needed in the formulation of Meseguer & Trefethen since the computation is performed with as independent variable, which would allow the collocation to resolve the near-centre of the pipe.

The no-slip boundary conditions are implemented through a set of row and column operations of reducing the order of the system as described through commented lines in the code provided in B. We would like to note in passing that the alternate method of spurious modes technique, which is a concise and clever way of implementing such homogeneous boundary conditions described in Ref. schmid2001stability, renders our results different at the order of . The method of elimination that we used gives number of modes in the cases of and , and modes in the case of . Care is taken to ward off round-off errors by evaluating the functions , , which are polynomials in , in a nested manner.

The eigensystems were solved using QZ algorithm implemented in the tool, "EIG" of Matlab (version R2016a Update 7), and the least decaying modes are evaluated through the Arnoldi method implemented in the tool, "EIGS" of Matlab which allows specification of the error tolerance.

### 3.2 Validations and comparisons

The obtained eigenvalues are compared with that of Meseguer & Trefethen meseguer2000spectral, Schmid & Henningson schmid2001stability and Priymak & Miyazaki priymak1998accurate with identical parameters as in those references, and shown in Table 1. It can be noted that, on the average, the accuracy matches with that of Meseguer & Trefethen meseguer2000spectral.

We also found that all 41 eigenvalues listed in meseguer2000spectral were matching at similar accuracy with the present computation. The figures for Schmid & Henningson reported in the Table 1 are after converting the complex phase-speeds reported in schmid2001stability into complex ’s by the relation . The improved accuracy in the present case is pronounced in comparison with Schmid & Henningson, which is due to the deployment of the ansatz of Priymak & Miyazaki priymak1998accurate given in Eq. (8). Schmid & Henningson’s schmid2001stability results were obtained by solving a system equivalent to that of Burridge & Drazin burridge1969comments.

As shown in Table 1, a matching accuracy is reached at collocation points as low as , which corresponds to the system size of . Increasing the value of further did not increase the accuracy of the least decaying mode in the present 64-bit calculations, although the modes with higher decay rates improved in convergence.

The attaining of the maximum possible accuracy of the least decaying mode for the case of shown in Table 1, is due to that the system is normal, hence the effectiveness of eigenvalue iterative algorithms is maximized. In this case, the Arnoldi method reduces to the Lanczos method golub2012matrix.

The eigenvalue of the least decaying mode of the case with parameters, and matches well with that of Priymak & Miyazaki priymak1998accurate. It should be noted that the present computations have been performed in 64-bit calculations. From the number of significant places in the data of Priymak & Miyazaki shown in Table 1, one can note that it is a result from 80-bit computation.

Figure 1 shows representative spectra for various combinations of parameters.

The sub-figures show the overall converged spectrum. As can be observed from the caption of Fig. 1, an increase in the number of collocation points is needed when there is an increase in , shown in Fig. 1(c), or , Fig. 1(b), but seldom for an increase in , shown in Fig. 1(e) and (f). However, it should be observed that there is a surfacing of mild distortion in the spectrum for with and as shown in Fig. 1(f). This can be explained using Eq. (4). It is well-known that the non-normality is caused by the convective terms that interact with mean shear. These terms are enhanced for large values of by the first term on the RHS of Eq. (4) giving rise to such distortions. The spectrum in Fig. 1(b) has same parameters as a sub-figure of Fig. 1 of Meseguer & Trefethen meseguer2003linearized and can be compared as a validation. For the case of shown in Fig. 1(d), the two wall-mode branches appears to have merged at the accuracy of the figure. As the azimuthal wavenumber is increased the number of modes in the centre branch decrease as shown in Fig. 1(e) and (f).

Table 2 compares the spectrum obtained for the case of from the system of Eqs. (19)-(20) with those obtained from the characteristic relations given by Eqs. (35)-(36). These characteristic relations have been approximated by polynomials, as explained in the caption of Table 2, by setting a cut-off value of 90 for the running index in the summation.

In the table, we have shown only the converged decimal figures that does not change when is changed from to , or when is changed from 89 to 90. As can be noted from this table, the results obtained by computation from Eqs. (19)-(20) have higher converging precision in comparison with those given by the characteristic relations Eqs. (35)-(36). The characteristic relations contains a quadratic term of factorials of in the denominator, which causes loss of precision for large values of . The set value of is the largest that we could afford in the present 64-bit (double precision) calculations. Especially, the highly decaying modes are prone to the precision loss since are large for these modes.

Finally, in Fig. 2 we compare the result from the present system for large wavenumbers with that from solving the system of Burridge & Drazin burridge1969comments. The system of Burridge and Drazin is given by

 (ω−αU)T¯¯¯ϕ= −αdr−1D(rd−1Ur)¯¯¯ϕ +iRe−1(T2¯¯¯ϕ−2αnT¯¯¯¯Ω)and (42) (ω−αU)¯¯¯¯Ω= −n(rd)−1Ur¯¯¯ϕ +iRe−1(2αnd−2T¯¯¯ϕ+S¯¯¯¯Ω), (43)

where , and , and the operators and are defined as and . The boundary conditions are given in Schmid & Henningson schmid2001stability. (However, there is a trivial typographical sign error in the first term of Eq. (3.41) in that reference. The system of Burridge & Drazin as it appear in Ref. burridge1969comments, has a typographical error in the definition of operator, .)

As can be noted from Fig. 2(a) and (b), the computation from the system, Eq. (42)-(43), suffers from pseudospectra, which is absent in the present formulation. This can be understood as follows. In Fig. 2(a) and (b) the blue rectangle shows a region of the spectrum containing 25 modes that undergoes distortion when computation is performed from the system of Burridge and Drazin. Fig. 2(c) and (d) shows ’s of the modes that are present within those rectangles. As mentioned earlier, the components of eigenfunctions of the present formulation, for example, undergoes steep variations close to the centreline for large . Such variations have been suppressed in the case of the system of Eq. (42)-(43), as , which shows that Burridge and Drazin’s unknown is a multiplication of the present unknown by a function that is vanishingly small for a range of the domain. When is large, is closer to zero for almost a fifth of the pipe radius as can be seen from Fig. 2(d) for the modes of distorted regions of the spectrum shown in Fig. 2(b). This implies that the details of the factor of in the expression for , where the former, i.e., is the distinguishing feature among different eigenfunctions, would be poorly represented even at double precision. For evidence that the distortion of the spectrum in Fig. 2(b) originiates in round-off errors, one can note that this pattern is similar to Fig. 2(b) of Ref. schmid1994optimal, which shows a spectrum in 16 bit, although for a set of small values of parameters, i.e., , and .

Let us assume that we cast Eqs. (42)-(42) in the form of where , similar to the present formulation as we saw in Eq. (40). Let us compare the condition numbers, and where and . We found that and , for chosen from their respective rectangles shown in Fig. 2(a) and (b). This shows that the present eigensystem is much closer to singularity for these modes than that of Burridge and Drazin. In other words, the contours of resolvent norms in the complex plane will be of smaller value for the same radial distance from these eigenvalues in the present formulation in comparison with the system of Burridge and Drazin.

## 4 Inviscid Algebraic Instability of Axially constant Modes

In this section we deduce the present flow configuration’s analogue of the findings of Ellingsen and Palm ellingsen1975stability for plane shear flows for the streamwise constant modes. The aim is to find the initial and , i.e. at time , that maximizes the energy amplification, and the output pattern at a later time . Ellingsen and Palm solution concerns inviscid solutions with in a non-modal setting and results in a initial value problem. (For such initial value problem in viscous scenario, see Ref. bergstrom1992initial for series solutions, and Ref. schmid1994optimal for optimal patterns). For axially constant inviscid modes that has a general non-modal time dependence, Eq. (15) and (16) become

 ∂∂t[(ℓ+2)Dϕ+yD2ϕ]= 0 (44) ∂Ω∂t= −2inUyϕ. (45)

The Eq. (44) implies conservation of kinetic energy in the radial and azimuthal directions. As will be shown later in this section, the sum of energies in these directions is proportional to (apart from a factor of ). Eq. (44) governs the rate at which energy is pumped into the radial vorticity by the mean shear. Since the azimuthal velocity is conserved, this implies the rate at which the axial component of the kinetic energy is enhanced, giving rise to streaks.

Eq. (44) can be readily integrated to yield , where the function, is equivalent to the specification of an initial value of the kinetic energy in the non-streamwise directions. Such specification of at is equivalent to the specification of , the initial value of itself. (In fact, for a specified , the that is finite at centreline and satisfying is given by