A linear system for pipe flow stability analysis allowing for boundary condition modifications
AA linear system for pipe flow stability analysis allowing for boundary conditionmodifications
M. Malik a , Martin Skote b, ∗ a Block 32 Ghim Moh Link, b Cranfield University, United Kingdom
Abstract
An accurate system to study the stability of pipe flow that ensures regularity is presented. The system produces a spectrum thatis as accurate as Meseguer & Trefethen (2000), while providing flexibility to amend the boundary conditions without a need tomodify the formulation. The accuracy is achieved by formulating the state variables to behave as analytic functions. We showthat the resulting system retains the regular singularity at the pipe centre with a multiplicity of poles such that the wall boundaryconditions are complemented with precisely the needed number of regularity conditions for obtaining unique solutions. In thecase of axisymmetric and axially constant perturbations the computed eigenvalues match, to double precision accuracy, the valuespredicted 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.
Keywords:
Pipe flow, Flow stability, Optimal patterns, Algebraic growth
1. Introduction
Flow through pipes is of a great importance due to vast appli-cations ranging from cooling systems to fluid transportations.As drag reduction in such systems will enhance human life bylowering energy consumption, understanding laminar to turbu-lent 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 in-duced research to tackle the understanding in the perspectivesof nonlinear dynamics, chaos and intermediate edge states [1–3], relatively recently.However, for the nonlinearity to set in, the disturbances needto be brought to a certain finite-scale size. Algebraic transientgrowth [4–8] 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 in-viscid interaction of the mean flow and perturbation throughthe lift-up effect for the modes with azimuthal modulations andthrough the Orr mechanism of enhancing the perturbations thatoppose the mean shear in the case of modes with streamwisemodulations. Inherently, since the laminar mean flow exhibitsa mean shear, the flow has a tendency to transfer energy into in-finitesimal disturbances, which can be inflow borne or originatefrom wall roughness and vibrations.One of the other several routes to turbulence is distortionsof the mean flow. Spatial analysis of this flow configuration ∗ Corresponding author
Email addresses: [email protected] (M. Malik), [email protected] (Martin Skote)©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/ predicted instability under suitable mean flow distortion [9].Furthermore, irregularities in the pipe geometry can serve asa cause for instability. For example, the flow has been foundunstable in slowly diverging pipes [10]. In addition, flowthrough sinusoidally corrugated pipe has exhibited linear insta-bility [11]. Wall roughness in the molecular scales has beenfound to be a destabilizing factor compared to smooth pipe,though instability has not been observed [12].Among the early investigations that predicted the linear sta-bility of flow through straight pipes in temporal setting [13–16],Burridge & Drazin [14] proposed an ansatz for their workingvariables for an asymptotic analysis, which was later adoptedfor numerical calculations [7, 17, 18].Priymak & Miyazaki [19] predicted the form of the r − , theradial coordinate, dependence of the velocity and pressure vari-ables through an ansatz, which is suitable if solutions of thelinear system for the perturbations that are analytic at the cen-treline are sought. It should be noted that this ansatz already in-corporates the necessary behaviour for velocity fields deducedby Khorrami et al. [20], thus leaving only regularity as meansfor further conditioning of the system. We would like to notethat in the light of this ansatz, the working variables of Bur-ridge & Drazin [14], namely, φ and Ω (defined in section 3.2),which are related to the radial velocity and radial vorticity havethe property, ( φ, Ω) ∼ ( r , const ) for n = 0 to the leading or-der close to the centreline, and ( φ, Ω) ∼ ( r | n | , r | n | ) for n (cid:54) = 0 ,where n is the azimuthal wavenumber. This suggests that forhigher values of | n | (e.g., | n | ≥ ), all the eigenfunctions willbe vanishingly small in a finite neighbourhood of the centreline.Consequently, there will be precision loss in determining theeigenfunctions. This complication arises because the eigenval-ues are multiplied by exceedingly low values in the discretized Preprint submitted to Elsevier August 27, 2019 a r X i v : . [ m a t h . NA ] A ug igensystem equations close to the the centreline.Furthermore, an increase of collocation points would onlyworsen the situation by enhancing the round-off errors since theeigenfunctions ( φ, Ω) T do not undergo steep variations closeto the centreline. These round-off errors at the end of the it-erative procedures for eigenvalues demonstrate themselves aspseudospectra. Since pseudospectra generally affects modesdecaying at higher rate than the least decaying ones, the op-timal patterns can get distorted at such high values of n .Priymak & Miyazaki [19] also developed a method for solv-ing the linearized Navier-Stokes equations cast as a system oftwo unknowns, the radial and azimuthal perturbation velocities,that utilizes the above findings, thus producing very accuratespectra. However, the linear operator itself needs to be deter-mined in a complex algorithmic way, unlike the case of Bur-ridge & Drazin [14]. In [19], the linear system is determinedthrough a sequence of steps numerically, where the pressureneeds to be found by inverting a Poisson operator, and the con-tinuity is being imposed by a correcting pressure which leads tocorrecting velocities, similar to the steps of the SIMPLE algo-rithm.Meseguer & Trefethen [21, 22] deployed the ansatz of Priy-mak and Miyazaki, and have produced very accurate resultswith a number of collocation points as low as . Their work isin Petrov-Galerkin formalism which uses solenoidal trial basisand test functions that satisfies the particular boundary condi-tion of no-slip. However, this method comes with a task ofchanging the trial and test bases if the boundary conditions areother than no-slip. In some occasions, the boundary conditionsare even time dependent (see for example, [23]). Finding ap-propriate solenoidal trial and test bases that satisfy arbitraryboundary conditions and the properties prescribed by the ansatzof Priymak & Miyazaki is not a simple task, in general.The above facts motivate us to formulate a method that isas accurate and efficient as Priymak & Miyazaki [19] andMeseguer & Trefethen [21], but as flexible and simple as a usualspectral collocation method with 2-tuple state variable similarto that of Ref. [7] or [14]. To this end, we use the ansatz ofPriymak & Miyazaki for velocity fields to determine a 2-tupleworking variables that do not go like a power law close to thecentreline. As the boundary conditions will be imposed on theunknowns unlike the Petrov-Galerkin method, the formulatedequations are applicable for a range of boundary conditions.The theoretical derivations are presented in section 2, whilethe numerical results are discussed in section 3. Finally, in sec-tion 4, we derive the Ellingsen and Palm solutions [24] thatare valid for inviscid axially constant modes that captures non-modal algebraic growth. Derivation of this is remarkably sim-pler in the working variables in this paper. Such solutions havebeen found useful in the case of plane-shear flows to computethe optimal patterns that demonstrate the lift-up effect leadingto streaks. We find that the pipe flow configuration retains thesefeatures. Although this is known through earlier viscous com-putations [7] and DNS [25], we arrive at these results in theinviscid limit.
2. Theoretical derivations
Let us consider the cylindrical coordinates, r = ( x, r, θ ) ,which are the axial, radial and azimuthal coordinates, respec-tively. We emphasize the unusual order of the coordinates usedhere. Let u ( r ) = ( u, v, w ) T be the velocities normalized withrespect to the centreline velocity. Let us use the decomposi-tion, u = U + ˜ u where the mean flow, U = ( U, , T and U = 1 − r , and the considered perturbation, ˜ u is such that (cid:107) ˜ u (cid:107) (cid:28) (cid:107) U (cid:107) , which are governed by, ˜ u t + ( U · ∇ ) ˜ u + ( ˜ u · ∇ ) U = − ∇ ˜ p + Re − ∇ ˜ u , (1)and by the continuity condition, ∇ · ˜ u = 0 . In Eq. (1), Re = U c R/ν , where U c is the dimensional centreline ve-locity, R is the pipe radius, and ν is the kinematic viscos-ity. The ˜ p in Eq. (1) is the perturbation pressure. Substituting ( ˜ u , ˜ p ) = ( u (cid:48) , p (cid:48) ) exp[i( αx + nθ − ωt )] into Eq. (1) to consider asingle Fourier mode, and upon taking the curl to remove the ap-pearance of pressure, we obtain the following three equations: i( αU − ω )( Du (cid:48) − i αv (cid:48) ) = − i αU r u (cid:48) − ( U rr + U r D ) v (cid:48) + Re − (cid:0) ∆( Du (cid:48) − i αv (cid:48) ) − nr − η (cid:48) (cid:1) (2) i( αU − ω )[ D ( rw (cid:48) ) − i nv (cid:48) ] = − i αrU r w (cid:48) + Re − (cid:0) ∆[ D ( rw (cid:48) ) − i nv (cid:48) ] (cid:1) (3) i( αU − ω ) η (cid:48) = − i nU r r − v (cid:48) + Re − (cid:0) ∆ η (cid:48) + 2 nr − ( αv (cid:48) + i Du (cid:48) ) (cid:1) (4)where ∆ = (cid:2) D + r − D − r − ( d + 1) (cid:3) , ∆ = (cid:2) D − r − D − r − ( d − (cid:3) , η (cid:48) = i( nr − u (cid:48) − αw (cid:48) ) and d = n + α r . As a matter of convention, we use suffix r to imply differentiation of mean flow variables with respect to r , and the operator D to represent the same for the fluctuationvariables. Using the above definition of η (cid:48) and the continuityequation, the variables u (cid:48) and w (cid:48) can be written as u (cid:48) =i rd − ( αv (cid:48) + αrDv (cid:48) − nη (cid:48) ) and (5) w (cid:48) =i d − ( nv (cid:48) + nrDv (cid:48) + αr η (cid:48) ) . (6)Priymak & Miyazaki [19] observed the following behaviourclose to the centreline. ( u (cid:48) , v (cid:48) , w (cid:48) ) = (cid:26) ( ψ , rφ, rψ ) for n = 0 , ( r (cid:96) +1 ψ , r (cid:96) φ, r (cid:96) ψ ) for n (cid:54) = 0 (7)where (cid:96) = | n | − and, ψ , ψ and φ are analytic functionshaving Taylor expansion around the centreline with vanishingcoefficients of odd powers of r . Substituting Eq. (7) into thedefinition of η (cid:48) , we get its behaviour. In summary, v (cid:48) and η (cid:48) have the following forms: ( v (cid:48) , η (cid:48) ) = (cid:26) ( rφ, r Ω) for n = 0 and ( r (cid:96) φ, r (cid:96) Ω) for n (cid:54) = 0 (8)2ith Ω an analytical function having similar Taylor expansionas φ . Substituting Eq. (8) into Eq. (5) and Eq. (6), we get u (cid:48) = (cid:26) i r (cid:96) +1 d [ α ( (cid:96) + 1) φ + αrDφ − n Ω] n (cid:54) = 0 i α (2 φ + rDφ ) n = 0 (9) w (cid:48) = (cid:26) i r (cid:96) d − [ n ( (cid:96) + 1) φ + nrDφ + αr Ω] n (cid:54) = 0 , i α − r Ω n = 0 . (10)First, let us consider the case of n (cid:54) = 0 . Upon substitutingEqs. (8)-(10) into Eqs. (2)-(4) we get ( ω − αU ) ( αL φ − nL Ω) = αL φ − αnL Ω+ i Re − ( αL φ − nL Ω) , (11) ( ω − αU ) ( nL φ + αL Ω) = nL φ − α L Ω+ i Re − ( nL φ + αL Ω) , and (12) ( ω − αU )Ω = nr − U r φ + i Re − (2 αnL φ + L Ω) (13)where the operators, L − are given in Appendix A. Eq. (11)and Eq. (12) are fourth order in φ and third order in Ω . As wehave three equations, Eqs. (11)-(13) for two unknowns, φ and Ω , we use Eq. (11) and Eq. (12) to reduce the order of Ω . Theresulting equation that is fourth order in φ and second order in Ω is ( ω − αU ) (cid:0) rdL φ + 2 αnr d − Ω (cid:1) = αr ( U r − rU rr ) φ + i Re − ( L φ + 4 αnL Ω) , (14)where L and L are operators given in Appendix A.Four more tasks are carried out before arriving at the requiredfinal system to work with: Firstly, the appearance of secondderivative of Ω in Eq. (14) is removed with the help of Eq. (13).Though this will not reduce the order of the system, it helps inderiving 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 thegoverning equation.Secondly, the even parity nature of φ and Ω , which have Tay-lor series expansions φ ( r ) = (cid:80) n a n r n and Ω( r ) = (cid:80) n b n r n is implemented by the transformation, y = r . Since this willtransform φ and Ω to general analytic functions, Chebyshevpolynomials of both odd and even orders can be used for thespectral expansion without any loss of efficiency. (For anothermethod, where only even Chebyshev polynomials are used, seePriymak & Miyazaki [19] and Meseguer & Trefethen [21]).Thirdly, similar operations performed for the case of { n =0 , α (cid:54) = 0 } by substituting φ and Ω from Eq. (8), give therequired system for axisymmetric disturbances. However, itshould 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 befully decoupled in this case with six as the order of the system.Finally, a system suitable for the case of { n = 0 , α = 0 } ,the axisymmetric and axially constant modes, is derived with y ( ≡ r ) as the independent variable. This special case requiresa different set of working variables, as v (cid:48) ( y ) = 0 (hence, φ =0 ) for all y due to continuity, and η (cid:48) ( y ) = 0 (hence, Ω = 0 )by definition. We choose the variables, ψ and ψ defined inEq. (7) for this purpose. The required system for this case canbe derived from the axial and azimuthal components of Eq. (1). In summary, the final system of equations are ( ω − αU ) (cid:0) [ − α g + g D + 4 yd D ] φ − αnd Ω (cid:1) = − αd n U y φ + i Re − (cid:0) [ α g − α g D + g D + 8 y ( g + 4 d ) D + 16 y d D ] φ − α n [ d + 2 dy D ]Ω (cid:1) , (15) ( ω − αU ) d Ω = 2 nd U y φ + i Re − (cid:0) αn { α g − d + ( (cid:96) + 3) d ] D − yd D } φ + {− α g +4 d [( (cid:96) + 1) d + n ] D + 4 yd D } Ω (cid:1) , (16)for the case of n (cid:54) = 0 , and ( ω − αU )(8 D + 4 y D − α ) φ = i Re − (cid:0) y D +96 y D + (96 − α y ) D − α D + α (cid:1) φ (17) ( ω − αU )Ω = i Re − (cid:0) y D + 8 D − α (cid:1) Ω (18)for the case of { n = 0 , α (cid:54) = 0 } , and − i ωψ =4 Re − (cid:0) y D + D (cid:1) ψ (19) − i ωψ =4 Re − (cid:0) y D + 2 D (cid:1) ψ (20)for the case of { n = 0 , α = 0 } , where y ≡ r , D = dd y andthe d and g i ( y ) ( i = 1 · · · ) are functions as given in Ap-pendix A. The order of the system is six for each of the casesof n (cid:54) = 0 and { n = 0 , α (cid:54) = 0 } , whereas it is four in the case of { n = 0 , α = 0 } . In the latter case, the order is only four owingto the fact that we did not have to take the curl of Eq. (1) toremove the appearance of ∇ p (cid:48) , as pressure is a constant in theaxial and azimuthal directions similar to the other perturbationquantities. All modes of this flow configuration, for all cases of α and n , are observed to be decaying. For the special case of { n = 0 , α = 0 } , a physical reason can be deduced since thetime evolution is dictated only by the viscous dissipation be-cause both the mean shear and the pressure gradient terms arezero. The systems given by Eqs. (15)-(16) and Eqs. (17)-(18) needto be solved with six conditions for a unique solution, and fourconditions for the system of Eqs. (19)-(20). In the case of n (cid:54) = 0 and { n = 0 , α (cid:54) = 0 } , three of the needed six conditions arestraight forward: they are the no-slip and no-penetration con-ditions at the pipe boundary. These conditions in terms of thevariables φ and Ω read as φ = D φ = Ω = 0 at y = 1 . (21)3he other three conditions should come from the regularityof solutions at y = 0 . Note that, although we are seek-ing analytic solutions by making the required transformations, { v (cid:48) , η (cid:48) } → { φ, Ω } with relations given in Eq. (8), the appear-ance of regular singularity in Eq. (15)-(18) implies that theystill can allow non-analytic solutions such as, for analogy, theBessel functions of second kind Y ν for positive integer ν in thecase of the Bessel equation. In a numerical procedure, the trans-formation in Eq. (8) and the transformation of the coordinate, y = r would only guarantee the accuracy by shifting the focusto the factors of v (cid:48) and η (cid:48) , namely, φ and Ω , whose r − variationsare unknown, and allow us to work with a reduced number ofChebyshev polynomials in the spectral expansion, while stillallowing non-analytic solutions. Therefore, we need regularityconditions to explicitly rule out non-analytic solutions in oursolution procedure.Since we seek φ and Ω as analytic functions, they have Taylorseries expansion around y = 0 . This warrants all orders ofthe derivatives to be of O (1) , which implies that all derivativeterms in Eq. (15)-(18) that are multiplied by y should vanish inthe limit y → . For analogy, such requirement is satisfied bythe the analytic J ν , the Bessel functions of first kind for positiveinteger ν , in the case of the Bessel equation. This instantly givesus two conditions of regularity for each of the cases, which are ( ω − αU ) (cid:0) [ − α g + g D ] φ − αnd Ω (cid:1) = − αn d U y φ + i Re − (cid:0) [ α g − α g D + g D ] φ − α nd Ω (cid:1) , (22) ( ω − αU ) d Ω = 2 nd U y φ + i Re − (cid:0) αn { α g − d + ( (cid:96) + 3) d ] D} φ + {− α g +4 d [( (cid:96) + 1) d + n ] D} Ω (cid:1) , (23)for the case of n (cid:54) = 0 , and ( ω − αU )(8 D − α ) φ = i Re − × (cid:0) D − α D + α (cid:1) φ (24) ( ω − αU )Ω = i Re − (cid:0) D − α (cid:1) Ω (25)for the case of { n = 0 , α (cid:54) = 0 } . We need one more conditionfor each of these cases of n (cid:54) = 0 and { n = 0 , α (cid:54) = 0 } . Note thatin Eq. (15) and (17) the term D φ is being multiplied by y .Requiring lim y → y D φ = 0 would only yield the constraintthat y D φ ∼ O (1) . Therefore, the other condition is the onethat would enforce lim y → y D φ = 0 , since φ is analytic. Thefollowing two facts should be noted in order to derive the re-quired condition: (1) As the Eq. (15) and (17) are both valid forthe whole domain, their derivatives with respect to y are alsovalid through out the flow field; (2) Due to regular singularity,Eq. (22) and (24) are at most second order in φ and zero-th orderin Ω , allowing them to be treated as boundary conditions. Fur-thermore, the derivatives of Eq. (15) and (17) with respect to y at the centreline can be used as boundary conditions since theyare of at most first order in Ω and third order in φ , which are less by one order comparing to the orders of these variables in thegoverning equations. Hence, the additional regularity conditionthat guarantees the differentiability of φ until fourth order foreach cases of n (cid:54) = 0 and { n = 0 , α (cid:54) = 0 } are the following: ( ω − α ) (cid:0) [ − α g + α g D + g D ] φ − αn [2 α + n D ]Ω (cid:1) = [ α g − αn (cid:96) D ] φ + 2 α n (1 − D )Ω+ i Re − (cid:0) [ α g + α g D − α g D + g D ] φ − α n [ α ( (cid:96) −
1) + n ( (cid:96) + 3) D ]Ω (cid:1) (26)for n (cid:54) = 0 , and ( ω − α )(12 D − α D ) φ = α ( α − D ) φ + i Re − (192 D − α D + α D ) φ (27)for the case of { n = 0 , α (cid:54) = 0 } , where g ··· are constantsdefined in Appendix A. In Eq. (26), we have substituted thevalues of mean flow variables for brevity.Note that the multiplicity of the pole in the governing equa-tion is such that it yields exactly three conditions to comple-ment the three boundary conditions at the wall. Had we insteadderived another regularity condition by differentiating the gov-erning equations once more and evaluated at the pipe centre,the condition would not have served as a boundary conditionsince it would have had the same orders of derivatives of theunknowns as they appear in the governing equations. Similarly,had the multiplicity been lower than in the present case, wewould not have been able to obtain a total of six conditions tosolve the sixth-order system.In the case of { n = 0 , α = 0 } , the boundary and regularityconditions are given by ψ = ψ = 0 at y = 1 , and , (28) (cid:0) ω − Re − D (cid:1) ψ = (cid:0) ω − Re − D (cid:1) ψ = 0 . (29)at y = 0 .The regularity condition for ψ in Eq. (29) and thegoverning equation Eq. (20) can also be obtained from the limit α → in Eq. (25) and Eq. (18), respectively, since ψ = i α − Ω for n = 0 . Therefore, one can expect that a part of the spectrumfor the case of { n = 0 , α = 0 } can be obtained from the systemfor the case of { n = 0 , α (cid:54) = 0 } in the limit of α → . How-ever, there are a new set of modes for this case that cannot beobtained from the system for finite α in the limit α → . Thesemodes originate from Eq. (19). Both sets of these modes aregoverned by Sturm and Liouville theory (see for example [26])since the underlying operators are self-adjoint. In the next sub-section, the solutions for this case { α = n = 0 } are derived in y − coordinate so that the characteristic relation are presented ina form that is less susceptible to numerical errors. This will fa-cilitate to validate the eigenvalues obtained from the numericalspectral collocation solutions of Eq. (19) and (20).It should be highlighted that the boundary conditions canbe changed depending on the flow configurations, without af-fecting the formulation, i.e., Eq. (15)-(20). This is the advan-tage of the present system compared with Meseguer and Tre-fethen [21], where one would be tasked with finding appro-priate basis and test functions that satisfies the boundary and4egularity conditions, besides satisfying the other condition ofcontinuity. The operators on the RHS in Eq. (19)- (20) are examplesof Stokes operators. The general Stokes operators are definedas Laplacians acted upon by Leray projectors which ensurescontinuity [27]. In the present case of { α = 0 , n = 0 } , thecontinuity is automatically satisfied without any constraint as v (cid:48) = 0 , and hence ∇ p (cid:48) = 0 . This implies that the Leray projec-tor is an identity operator in the present case. It is well knownthat the eigenvalues of a general Stokes operators are decayingand that the operator itself is self-adjoint for no-slip boundaryconditions. For this particular case of { α = 0 , n = 0 } , the so-lutions and characteristic relation are known to be Bessel func-tions of r (cid:112) |(cid:61){ ω }| Re and their roots at r = 1 [15, 28]. (Forsuch modes of plane Poiseuille flow, see [29]. For applicationof these Stokes modes, see [30].)Eq. (19)- (20) can be rewritten as λψ = (cid:0) y D + D (cid:1) ψ , (30) λψ = (cid:0) y D + 2 D (cid:1) ψ , (31)where λ = − i ωRe/ . The differential operators on the RHSare self-adjoint under the inner product defined as the cross-sectional area integral of velocity fields, i.e., ψ or √ yψ . TheSturm and Liouville theory suggests that the eigenvalues, λ j ofEq. (30) and λ k of Eq. (31) with j, k ∈ , , , · · · are real andthe eigenfunctions, ψ ,j and ψ ,j are orthogonal in the sense, (cid:82) ψ ∗ ,j ψ ,k d y = C δ jk (cid:82) ψ ∗ ,j ψ ,k y d y = C δ jk (cid:41) (32)where δ jk is the Kronecker delta and, C and C are constants.As we are interested in solutions of Eq. (30)-(31) that areanalytic at the centreline, we can bypass the method of Frobe-nius and resort to a simple power series method. Let ψ = (cid:80) k ≥ a k y k and ψ = (cid:80) k ≥ b k y k , where the constants, a k and b k , are to be determined. Substituting these into Eq. (30)-(31), and matching the coefficients of like powers of y , we getthe eigenfunctions, ψ ,l and ψ ,l as, ψ ,l ( y ) = a (cid:88) k ≥ (cid:0) λ ,l y (cid:1) k ( k !) − and (33) ψ ,l ( y ) = b (cid:88) k ≥ (cid:0) λ ,l y (cid:1) k [ k !( k + 1)!] − , (34)where λ ,l with l = 1 , , , · · · are the real solutions of thecharacteristic equation, (cid:88) k ≥ λ k ( k !) − = 0 , (35)and λ ,l are the real solutions of (cid:88) k ≥ λ k [ k !( k + 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 J ( (cid:112) |(cid:61){ ω }| Re ) = 0 and J ( (cid:112) |(cid:61){ ω }| Re ) = 0 [15, 28],however the above expressions are helpful to obtain the eigen-values accurately as the roots of polynomials that approximatethe LHS of Eqs. (35) and (36) by cutting off the summation atsome k = k max . The constants, a and b in Eqs. (33) and (34),respectively, can be of any value as long as the linearity ofthe perturbations with respect to the mean flow is respected.Without loss of generality they can be the normalization con-stants, a = 1 / (cid:112) ( ψ , ψ ) and b = 1 / (cid:112) (cid:104) ψ , ψ (cid:105) , where theinner-products are defined as ( ψ , ψ ) = (cid:82) ψ ∗ ( y ) ψ ( y ) d y and (cid:104) ψ , ψ (cid:105) = (cid:82) yψ ∗ ( y ) ψ ( y ) d y . It should be mentionedthat (cid:82) d y is in fact an area element, (cid:82) r d r . Under suchnormalizations, the constants C and C of Eq. (32) become, C = C = 1 . The convergence of the series in Eq. (33)and (34) is obvious by comparison tests with the series for exp[ λ ,l y ] and exp[ λ ,l y ] , respectively. The implication ofthe orthogonality is that superposition of such modes cannothave transient growth. A velocity state can be represented us-ing the notation introduced at the beginning of this section as u (cid:48) kl ( y ) = { ψ ,k ( y ) , , √ yψ ,l ( y ) } T . Then the superpositionvelocity is given by, ˜ u ( t, y ) = exp(i[ αx + nθ ]) (cid:88) k,l a kl exp ( Λ kl t ) u (cid:48) kl , (37)where Λ kl = diag { λ k , , λ l } and a kl = diag { a k , , a l } arediagonal matrices for each chosen k and l , where in turn λ k and λ l are the real solutions of Eq. (35) and (36), respectively, and a k and a l are some constants.Let us define the energy as, E ( t ) = 2 − (cid:82) | ˜ u | dy . E ( t ) = 2 − (cid:88) k,l (cid:88) i,j a H ij a kl × exp ([ Λ kl + Λ ij ] t ) (cid:90) u (cid:48) ij H u (cid:48) kl d y, (38)where the superscript, ‘H’ denotes Hermitian transpose. Usingthe orthogonality relations given by Eq. (32) and normalizationsof ψ and ψ , we get, E ( t ) = 12 (cid:32)(cid:88) k | a k | exp (2 λ k t ) + (cid:88) l | a l | exp (2 λ l t ) (cid:33) . (39)As said in the subsection 2.1, in this present case of { α =0 , n = 0 } , the evolution of perturbations are determined bya dissipative phenomena. Therefore, λ k and λ l are all negative.This shows that the energy, E ( t ) in Eq. (39), is a monotonouslydecaying function of time. Since the interaction with the meanflow is null for these modes, they should be supplied with en-ergy by other means of receptivity such as vibration of the wallor by inflow-borne disturbances. The role of acoustic distur-bances may also be ruled out since extremely longwaves implyextremely low (near-zero) infrasonic frequencies of the source(since, frequency ∼ /wavelength).5 . Numerical Results The system of Eqs. (15)-(18) can be written as, Re − Aq = ω Bq , (40)where q = { φ, Ω } T and, A and B are matrices of operators.The elements of A and B can be figured out from those equa-tions (Eq. (15)-(18)) for each case of n (cid:54) = 0 and n = 0 . Fordiscretization, we use spectral collocation method in transformspace (see, for example, Ref. [31] or Appendix A.5 of Ref. [32]and S. C. Reddy’s codes therein). The generalized eigensystem,Eq. (40) is discretized at the points, y j = 2 − (1 + ξ j ) , if α ≤ and Re ≤ , (41a)else y j = (exp( a ) − − (exp[ a (1 + ξ j ) / − , (41b)where the Gauss Lobatto points, ξ j = cos( πj/N ) , j =0 · · · N and N is the highest order of Chebyshev polynomials.Eq. (41b) represents a stretching function that maps ξ on y suchthat the grid points cluster towards y = 0 for the parameter a > . For n ≤ , a = 2 gives best result, whereas a = 3 works well for n > .The stretching is needed in our present formulation for largevalues of α and Re , due to the following. As the algebraic r − variations of the velocity fields in a power-law fashion (i.e., y − variations) has been factored out from the unknowns by thePriymak and Miyazaki ansatz given in Eq. (7), the functions φ and Ω do not have any other constraints apart from those of sat-isfying regularity. This allows Ω and φ to have steep variationsat the pipe centre, given that there is a regular singularity at thatlocation. The freedom of Ω and φ and their derivatives from therequirement to vanish at the pipe centre, coupled with their van-ishing at the wall and the presence of regular singularity there,causes a boundary layer behaviour at pipe centre. Hence, thestretching introduced in Eq. (41b) is required in our formalism.However, it is precisely this boundary layer behaviour at thepipe centre that allows for avoiding the pseudospectra, as willbe shown further down.Note that the requirement for stretching does not arise in theformulation of Burridge & Drazin, where the functions and oneof their derivatives vanish at both boundaries, namely, wall andpipe centre. Therefore, upon solving their system, the varia-tions of the solutions are spread over the entire domain, avoid-ing steep variations. Such clustering is neither needed in theformulation of Meseguer & Trefethen since the computation isperformed with r as independent variable, which would allowthe collocation to resolve the near-centre of the pipe.The no-slip boundary conditions are implemented through aset of row and column operations of reducing the order of thesystem as described through commented lines in the code pro-vided in Appendix B. We would like to note in passing that thealternate method of spurious modes technique, which is a con-cise and clever way of implementing such homogeneous bound-ary conditions described in Ref. [32], renders our results differ-ent at the order of − . The method of elimination that we used gives N − number of modes in the cases of n (cid:54) = 0 and { n = 0 , α (cid:54) = 0 } , and N modes in the case of { n = 0 , α = 0 } .Care is taken to ward off round-off errors by evaluating thefunctions g i , i = 1 · · · , which are polynomials in α y , in anested manner.The eigensystems were solved using QZ algorithm imple-mented in the tool, EIG of Matlab (version R2016a Update 7),and the least decaying modes are evaluated through the Arnoldimethod implemented in the tool,
EIGS of Matlab which allowsspecification of the error tolerance.
The obtained eigenvalues are compared with that ofMeseguer & Trefethen [21], Schmid & Henningson [32] andPriymak & Miyazaki [19] with identical parameters as in thosereferences, and shown in Table 1. It can be noted that, on theaverage, the accuracy matches with that of Meseguer & Tre-fethen [21]. We also found that all 41 eigenvalues listed in [21]were matching at similar accuracy with the present computa-tion. The figures for Schmid & Henningson reported in theTable 1 are after converting the complex phase-speeds reportedin [32] into complex ω ’s by the relation ω = αc . The improvedaccuracy in the present case is pronounced in comparison withSchmid & Henningson, which is due to the deployment of theansatz of Priymak & Miyazaki [19] given in Eq. (8). Schmid& Henningson’s [32] results were obtained by solving a systemequivalent to that of Burridge & Drazin [14].As shown in Table 1, a matching accuracy is reached at col-location points as low as N + 1 = 48 , which corresponds to thesystem size of N − . Increasing the value of N fur-ther did not increase the accuracy of the least decaying mode inthe present 64-bit calculations, although the modes with higherdecay rates improved in convergence.The attaining of the maximum possible accuracy of the leastdecaying mode for the case of α = n = 0 shown in Table 1,is due to that the system is normal, hence the effectiveness ofeigenvalue iterative algorithms is maximized. In this case, theArnoldi method reduces to the Lanczos method [33].The eigenvalue of the least decaying mode of the case withparameters, α = n = 20 and Re = 4000 matches well withthat of Priymak & Miyazaki [19]. It should be noted that thepresent computations have been performed in 64-bit calcula-tions. From the number of significant places in the data of Priy-mak & Miyazaki shown in Table 1, one can note that it is aresult from 80-bit computation.Figure 1 shows representative spectra for various combina-tions of parameters. The sub-figures show the overall convergedspectrum. As can be observed from the caption of Fig. 1, an in-crease in the number of collocation points is needed when thereis an increase in α , shown in Fig. 1(c), or Re , Fig. 1(b), but sel-dom for an increase in n , shown in Fig. 1(e) and (f). However,it should be observed that there is a surfacing of mild distortionin the spectrum for n = 10 with α = 1 and Re = 3000 asshown in Fig. 1(f). This can be explained using Eq. (4). It iswell-known that the non-normality is caused by the convectiveterms that interact with mean shear. These terms are enhanced6 able 1: Least decaying eigenvalues of converged accuracy. ‘Size’ refers to the order of eigensystemParameters Present Method Meseguer & Trefethen(2000) α n Re Size ω r ω i Size ω r ω i . − . . − . . − . . − . . − . . − . . − . . − . − . − . − . − . − . − . − . − . . − . . − . Parameters Present Method Schmid & Henningson(2001) α n Re
Size ω r ω i Size ω r ω i . − . − . − . . . − . − . − . .
25 2 2000 77 0 . − . − . − . − . − − . Parameters Present Method Priymak & Miyazaki(1998) α n Re
Size ω r ω i Size ω r ω i
20 20 4000 201 1 . − . . · · · − . · · ·· · · · · · -1-0.50 c i Re =3000 α =1 , n =1 Re =10000 α =1 , n =1 Re =3000 α =3 , n =1 c r -1-0.50 c i Re =3000 α =1 , n =0 c r Re =3000 α =1 , n =5 c r Re =3000 α =1 , n =10 (a) (b) (c)(d) (e) (f) Figure 1: Spectra for parameters shown in each panels. Other parameters:(a) N = 47 ; (b) N = 68 ; (c) N = 69 ; (d) N = 53 ; (e) N = 53 ; (f) N = 53 ; for large values of n by the first term on the RHS of Eq. (4)giving rise to such distortions. The spectrum in Fig. 1(b) hassame parameters as a sub-figure of Fig. 1 of Meseguer & Tre-fethen [22] and can be compared as a validation. For the caseof { n = 0 , α = 1 } shown in Fig. 1(d), the two wall-modebranches appears to have merged at the accuracy of the fig-ure. As the azimuthal wavenumber is increased the numberof modes in the centre branch decrease as shown in Fig. 1(e)and (f).Table 2 compares the spectrum obtained for the case of { α = 0 , n = 0 } from the system of Eqs. (19)-(20) with thoseobtained from the characteristic relations given by Eqs. (35)-(36). These characteristic relations have been approximated bypolynomials, as explained in the caption of Table 2, by settinga cut-off value of 90 for the running index k in the summation.In the table, we have shown only the converged decimal fig-ures that does not change when N is changed from to , orwhen max { k } is changed from 89 to 90. As can be noted fromthis table, the results obtained by computation from Eqs. (19)- Table 2: Comparison of first seven least decaying eigenvalues, ω i ’s for the caseof { α = 0 , n = 0 } with Re = 3000 and N = 47 obtained from system ofEq. (19)-(20) versus those obtained from ω i = λRe/ where the λ ’s are thereal solutions of the characteristic Eq. (35)-(36) approximated with the setting max { k } = 90 . The reported figures are of the converged accuracy with respectto N and max { k } . ω i ’s from Eq. (19) λRe/ with λ ’s from Eq. (35) − · − · − · − · − · − · − · − · − · − · − · − · − · − · ω i ’s from Eq. (20) λRe/ with λ ’s from Eq. (36) − · − · − · − · − · − · − · − · − · − · − · − · − · − · (20) have higher converging precision in comparison with thosegiven by the characteristic relations Eqs. (35)-(36). The charac-teristic relations contains a quadratic term of factorials of k inthe denominator, which causes loss of precision for large valuesof k . The set value of max { k } = 90 is the largest that we couldafford in the present 64-bit (double precision) calculations. Es-pecially, the highly decaying modes are prone to the precisionloss since | λ | are large for these modes.Finally, in Fig. 2 we compare the result from the present sys-tem for large wavenumbers with that from solving the systemof Burridge & Drazin [14]. The system of Burridge and Drazinis given by ( ω − αU ) T φ = − αdr − D (cid:0) rd − U r (cid:1) φ + i Re − (cid:0) T φ − αnT Ω (cid:1) and (42)7 c r -1-0.8-0.6-0.4-0.20 c i Present c r -1-0.8-0.6-0.4-0.20 From thesystem of BD -1 0 1 φ / || φ || ∞ r Present -1 0 1 φ / || φ || ∞ From thesystemof BD (a) (b)(c) (d)
Figure 2: Comparison of spectra and eigenfunctions from the present formula-tion with that computed from the system of Burridge & Drazin, Eq. (42)-(43).Parameters: N = 200 , α = 10 , n = 7 , Re = 2000 . (a) and (b) are spectra;(c) and (d) are the ∞− norm normalized real parts of φ and φ , respectively,close to the centreline of 25 modes selected from the region shown by rectan-gles in (a) and (b). In (c) φ has been plotted with respect to r for the purpose ofcomparison with (d). ( ω − αU )Ω = − n ( rd ) − U r φ + i Re − (cid:0) αnd − T φ + S Ω (cid:1) , (43)where D = dd r , φ = − i rv (cid:48) and Ω = ( αrw (cid:48) − nu (cid:48) ) /d , and theoperators T and S are defined as T [ · ] = dr − D ( rd − D [ · ]) − dr − [ · ] and S [ · ] = ( rd ) − D ( rdD [ · ]) − dr − [ · ] . The boundaryconditions are given in Schmid & Henningson [32]. (However,there is a trivial typographical sign error in the first term ofEq. (3.41) in that reference. The system of Burridge & Drazinas it appear in Ref. [14], has a typographical error in the defini-tion of operator, T .)As can be noted from Fig. 2(a) and (b), the computation fromthe system, Eq. (42)-(43), suffers from pseudospectra, which isabsent in the present formulation. This can be understood asfollows. In Fig. 2(a) and (b) the blue rectangle shows a re-gion of the spectrum containing 25 modes that undergoes dis-tortion when computation is performed from the system of Bur-ridge and Drazin. Fig. 2(c) and (d) shows φ ’s of the modes thatare present within those rectangles. As mentioned earlier, thecomponents of eigenfunctions of the present formulation, φ forexample, undergoes steep variations close to the centreline forlarge α . Such variations have been suppressed in the case of thesystem of Eq. (42)-(43), as φ = − i r n φ , which shows that Bur-ridge and Drazin’s unknown φ is a multiplication of the presentunknown φ by a function that is vanishingly small for a rangeof the domain. When n is large, φ is closer to zero for almosta fifth of the pipe radius as can be seen from Fig. 2(d) for themodes of distorted regions of the spectrum shown in Fig. 2(b). This implies that the details of the factor of φ in the expressionfor φ (= − i r n φ ) , where the former, i.e., φ is the distinguishingfeature among different eigenfunctions, would be poorly repre-sented even at double precision. For evidence that the distortionof the spectrum in Fig. 2(b) originiates in round-off errors, onecan note that this pattern is similar to Fig. 2(b) of Ref. [7], whichshows a spectrum in 16 bit, although for a set of small values ofparameters, i.e., α , n and Re .Let us assume that we cast Eqs. (42)-(42) in the form of Re − Aq = ω Bq where q = ( φ, Ω) T , similar to the presentformulation as we saw in Eq. (40). Let us compare the con-dition numbers, cond ( L − ωI ) and cond ( L − ωI ) where L = B − A /Re and L = B − A /Re . We found thatcond ( L − ωI ) ∼ O (10 ) and cond ( L − ωI ) ∼ O (10 ) , for ω (cid:48) s chosen from their respective rectangles shown in Fig. 2(a)and (b). This shows that the present eigensystem is much closerto singularity for these modes than that of Burridge and Drazin.In other words, the contours of resolvent norms in the complexplane will be of smaller value for the same radial distance fromthese eigenvalues in the present formulation in comparison withthe system of Burridge and Drazin.
4. Inviscid Algebraic Instability of Axially constant Modes
In this section we deduce the present flow configuration’sanalogue of the findings of Ellingsen and Palm [24] for planeshear flows for the streamwise constant modes. The aim is tofind the initial φ and Ω , i.e. at time t = 0 , that maximizes theenergy amplification, and the output pattern at a later time t .Ellingsen and Palm solution concerns inviscid solutions with α = 0 in a non-modal setting and results in a initial valueproblem. (For such initial value problem in viscous scenario,see Ref. [34] for series solutions, and Ref. [7] for optimal pat-terns). For axially constant inviscid modes that has a generalnon-modal time dependence, Eq. (15) and (16) become ∂∂t (cid:2) ( (cid:96) + 2) D φ + y D φ (cid:3) =0 (44) ∂ Ω ∂t = − nU y φ. (45)The Eq. (44) implies conservation of kinetic energy in the ra-dial and azimuthal directions. As will be shown later in thissection, the sum of energies in these directions is proportionalto ( (cid:96) + 2) D φ + y D φ (apart from a factor of − n − y (cid:96) +1 ).Eq. (44) governs the rate at which energy is pumped into theradial vorticity by the mean shear. Since the azimuthal velocityis conserved, this implies the rate at which the axial componentof the kinetic energy is enhanced, giving rise to streaks.Eq. (44) can be readily integrated to yield ( (cid:96) + 2) D φ + y D φ = C ( y ) , where the function, C ( y ) = [( (cid:96) + 2) D φ + y D φ ] | t =0 is equivalent to the specification of an initial valueof the kinetic energy in the non-streamwise directions. Suchspecification of C ( y ) at t = 0 is equivalent to the specificationof φ , the initial value of φ itself. (In fact, for a specified C ( y ) ,the φ ( y ) that is finite at centreline and satisfying φ (1) = 0 is8iven by φ ( y ) = (cid:82) [ U ( y − y ) − y − ( (cid:96) +2) (cid:82) y C (˜ y )˜ y (cid:96) +1 d˜ y d y where U ( y ) is the unit step function). Therefore, φ ( t, y ) = φ ( y ) (46) Ω( t, y ) =Ω ( y ) − nU y tφ ( y ) , (47)where Ω ( y ) = Ω( t = 0 , y ) . Eq. (47) is obtained from in-tegration of Eq. (45) and describes algebraic growth since thesecond term is secular, which implies the collapse of linear per-turbation theory after certain initial time unless viscosity actsand kills it. In terms of velocities these equations translate into: v (cid:48) ( t, y ) = v (cid:48) (0 , y ) , u (cid:48) ( t, y ) = u (cid:48) (0 , y ) − t √ yU y v (cid:48) (0 , y ) and w (cid:48) ( t, y ) = i n − [ v (cid:48) (0 , y ) + 2 y D v (cid:48) (0 , y )] .To find the optimal patterns in this inviscid limit, we fol-low the method of Ref. [35] for compressible flow with nec-essary modifications for the present incompressible configura-tions (see also Ref [36]). As we are working with a 2-tuplevariable, the constraint of continuity is already taken into con-sideration. The energy, E ( t ) ≡ − (cid:82) | u (cid:48) | dy can be writtenwith the help of Eq. (8)-(10) as E ( t ) = 12 n (cid:90) (cid:2) y (cid:96) +1 | Ω | + 2 n y (cid:96) | φ | + 4 y (cid:96) +2 |D φ | +2( (cid:96) + 1) y (cid:96) +1 ( φ ∗ D φ + φ D φ ∗ ) (cid:3) d y, (48)which can be written after integration by parts as, E ( t ) = (2 n ) − (cid:90) y y (cid:96) +1 q H A H MA q d y, (49)where q = { φ , Ω } T , the positive definite × matrix M = diag {− (cid:96) +2) D + y D ] , } . The operator A is the propagatorof q and can be obtained by casting Eq. (46) and Eq. (47) inmatrix form. The elements of A are A = A = 1 , A = 0 ,and A = − nU y t . Let g ( t ) be the propagator of the initialenergy, i.e., E ( t ) = g ( t ) E (0) . (50)Note that E (0) = (2 n ) − (cid:82) y y (cid:96) +1 q H M q d y . Let us define, G ( t ) = max q g ( t ) . This G ( t ) is found as follows. Taking thefunctional derivative of Eq. (50), we get δδ q H E ( t ) = (cid:18) δδ q H g ( t ) (cid:19) E (0) + g ( t ) δδ q H E (0) . (51)Setting δg/δ q H = 0 for maximization we get, A H MA q = g ( t ) M q , (52)which is a differential equation that can also be written as (cid:2) ( nU y t ) − ( (cid:96) + 2) D − y D (cid:3) φ + i nU y t Ω = − g ( t ) (cid:2) ( (cid:96) + 2) D + y D (cid:3) φ (53) − nU y tφ + Ω = g ( t )Ω (54)This is an eigenvalue problem with the Lagrange multiplier, g ( t ) , as the eigenvalue for a chosen t . Since the operatorsare Hermitian, g ( t ) ’s are all real, and as M is positive definite Streamlines at t = 0 Axial velocity at t = 50 (a) -20-1001020 (b)(c) -50050 (d)(e) -80-60-40-20020406080 (f) Figure 3: The inviscid initial (left) and the instantaneous velocity patterns at alater time (right) for axially constant modes that are optimized for maximumenergy growth: (a) and (b) n = 1 ; (c) and (d) n = 4 ; and, (e) and (f) n = 10 .The sub-figures on the left column are the 2D streamlines in the cross sectionof the pipe at time, t = 0 . The sub-figures on the right column are the contoursof the axial velocity at time, t = 50 with levels indicated by colour-bar. by definition, g ( t ) ’s are all positive. The boundary conditionsare no-penetration at the wall and regularity at the centreline,i.e., φ (1) = 0 and g ( t ) − (cid:96) + 2) D + ( nU y t ) ] y =0 φ +i nU y t Ω (0) = 0 , respectively. As Ω appears algebraicallyin Eq. (54), it does not require boundary conditions. In fact,we can explicitly solve for Ω in terms of φ and expressEq. (52) as an eigenvalue problem with only φ as the un-known. However, this will lead to that the eigenvalue param-eter, g ( t ) , appears quadratically and will warrant a compan-ion matrix method to solve the equation numerically, which isequivalent to the problem as it appears currently in Eq. (52).The maximum amplification, G ( t ) as defined above, is givenby G ( t ) = max j g j ( t ) where g j ( t ) ’s for a chosen t are the setof eigenvalues of Eq. (52). However, the variation G ( t ) ∼ t can be anticipated from Eq. (47). In addition, such a variationof G ( t ) with respect to t has been observed in wall-boundedplane shear flows [37]. Here, we focus is on the optimal q and the corresponding q ( t ) as obtained from Eq. (52). To solvethe system given by Eq. (52), we adopt the the spectral methoddescribed in previous section.Figs. 3(a), (c) and (e) show the streamlines of initial veloc-ity pattern in the cross sectional plane for n = 1 , and , re-spectively, that are prone to the maximum energy amplification.These patterns represent counter rotating vortices. The expla-9ation for the appearance of these patterns is similar to the caseof plane shear flows, such as boundary layer, Couette and planePoiseuille flows. Since counter rotating vortices contain move-ment of fluid particles in the radial direction, they can transferenergy from mean shear into axial velocity. This phenomenonis widely known as lift-up effect [32] in these plane shear flows.In the context of this pipe flow, the term, up refers to the radialdirection.Such patterns have been observed in this flow configura-tion through DNS of full nonlinear viscous equations whenthe mean flow was perturbed initially with axially constantmodes [25]. Fig. 3(a) matches with the optimal pattern fromthe viscous flow computation at Re = 3000 shown in Ref. [22].Since this inviscid analysis captures those features, one canconclude that the inviscid lift-up effect is dominant at initialtimes before viscosity break this vortices into small scales.Figs. 3(b), (d) and (f) show the contours of axial componentof perturbation velocity at a later time chosen as t = 50 in thenon-dimensional scale, where the initial state q has been nor-malized such that E (0) = 1 . These patterns describe near-wallstreaks when superimposed with the mean flow. These streaksare known to be the characteristic feature of bypass transition .(For a matching of these optimal perturbations with experimen-tal data in the case of boundary layers, see [38]). As the valueof n is increased, these streaks are pushed closer to the wall, aresult tied to the following two facts: (1) The streamwise veloc-ity’s growth is proportional to the radial velocity as can be seenfrom the relation, u (cid:48) ( t, y ) = u (cid:48) (0 , y ) − t √ yU y v (cid:48) (0 , y ) ; and,(2) the radial velocity culminates at the periphery of the cross-section when n is increased as can be to noted from v (cid:48) = r l φ .It should be noted that the G ( t ) at a chosen time would in-crease with an increase in n as evident from Eq. (47). This alsoshows up on the range of the colour bars in Fig. 3 as we movedown from panels (b) to (f). However, as can be seen from theoptimal patterns (in the inviscid limit), the size of these vor-tices reduces when n is increased, making them more prone tothe action of viscosity to break them further and eventually be-ing killed. This phenomenon can be observed from Fig.4(a) ofRef. [7] and Fig.4 of Ref. [39].
5. Conclusion
The present formulation, which uses the ansatz of Priymak& Miyazaki [19] results in accurate spectra even for reason-ably large values of the parameters, namely, the wavenumbersand Reynolds number. The obtained accuracy is on par withthat of Meseguer & Trefethen [21]. However, as the boundaryconditions of the wall are directly imposed on the unknownsused, which is in contrast to that of Meseguer & Trefethen [21],where they are imposed on the basis functions, the present for-mulation enables the derived system to be applied on a range offlow configurations with various boundary conditions.The multiplicity of the poles of regular singularity yielded thenecessary number of regularity conditions in order to comple-ment the boundary conditions at wall, providing 6 conditionsto solve the 6th order system for the cases other than that of { α = 0 , n = 0 } . In the case of { α = 0 , n = 0 } , the working variables had tobe changed from a representative of normal vorticity and nor-mal velocity to that of streamwise and azimuthal velocities. Theobtained self-adjoint system’s analytical characteristic relationpredicted eigenvalues which were used to validate the numer-ical results for this case. The orthogonality showed that thesuperposition of these modes are always decaying.As a demonstration of the performance of this formulation,the results were contrasted against those from the system ofBurridge and Drazin [14]. The performance of the current for-mulation is found to be highlighted at large values of the param-eters, and found to postpone the appearance of pseudo-spectrato the regimes of further higher values of these parameters. Thiswas tracked to the fact that the numerical representation of dis-tinctive features of different eigenfunctions are enhanced in thepresent formulation, and that it undergoes precision loss in thecase of Burridge and Drazin due to round-off errors.Finally, the inviscid algebraic (nonmodal) growth of the per-turbation for streamwise constant modes were studied. The re-sult was the pipe flow’s counterpart of Ellingsen and Palm so-lutions for plane shear flows. The present working variablesfor the unknowns and the independent variable, y , resulted inthe solutions which were almost identical in structure to thatof plane shear-flows, hence retaining the well known character-istics such as having the counter rotating vortices and streaksas optimal patterns. These streaks were found to be closer tothe wall when the azimuthal wavenumber is larger. For the ap-pearance of such streaks, it is not necessary that the infinitesi-mal perturbations acquired through receptivity should have theFourier components with n very large. The nonlinearity canplay the role even at very early stages to give birth to suchmodes through convolution of modes with lower n , and sub-sequently amplified by the linear mechanism of lift-up . Appendix A. Operators and functions
To facilitate locating the functions and operatorsbelow, their names have been given in boldface. L = ( f + f D + rd − D ) , L = ( f + d − D ) , L = ( f + f D + α rU d − D ) , L = d − ( f + U D ) , L = ( f + f D + f D + f D + rd − D ) , L = ( f + f D + f D + d − D ) , L = f + f D + ( rd ) − D , L = f + d − D , L = f + U r ( rd ) − D , L = d − ( U r + U D ) , L = f + f D + f D + f D + ( rd ) − D , L = f + f D + f D + d − D , L = f − f D − d − D , L = f + f D + D , L = f + f D + f D + f D + r D , L = f + f D + r d − D f = [( (cid:96) + 1) d − d ] / ( d r ) , f = [ d + ( (cid:96) + 2) d ] /d , f = d / ( d r ) , f = [ U r ( α r − n (cid:96) ) − rdU rr ] / ( dr ) , f = − n U r / ( rd ) , f = ( U d + rdU r ) / ( rd ) , f = [( (cid:96) +1)( d + dd − ( d +1) d d )+( d +1 − (cid:96) ) d ] / ( d r ) , f = [ d + (3 (cid:96) + 7) dd + (2 (cid:96) + 3 − d ) d d − (cid:96) + 1) d − ( (cid:96) +2) d ] / ( d r ) , f = [3 d +(3 (cid:96) +11) dd +( (cid:96) +2 − d ) d ] / ( rd ) , f = [3 d + ( (cid:96) + 5) d ] /d , f = [ d + dd − ( d + 1) d d + d ] / ( d r ) , f = [3 d + 2 dd − ( d + 1) d ] / ( d r ) , f = [3 d + d ] / ( d r ) , f = [( (cid:96) + 1)( d + d ) − d ] / ( d r ) , f =[ d + ( (cid:96) + 3) d ] / ( d r ) , f = ( d + 3 d ) / ( rd ) , f = ( (cid:96) + 1) U r / ( r d ) , f = ( (cid:96) + 1)( r d ) − [ d +2 dd − ( d + 1) d d + (1 − d ) d ] + r − ( d − (cid:96) + 2 (cid:96) − , f = [ d + (3 (cid:96) + 8) dd + (4 (cid:96) + 7 − d ) d d − ( (cid:96) +1) d − (3 (cid:96) + 2) d ] / ( rd ) , f = [3 d + (3 (cid:96) + 13) dd +(2 (cid:96) + 5) d − d ] / ( rd ) , f = [3 d + ( (cid:96) + 6) d ] / ( rd ) , f = [ d + 8 dd + (13 − d ) d d + 3(1 − d ) d ] / ( d r ) , f =(3 d + 15 dd + 13 d − d ) / ( d r ) , f = (3 d + 8 d ) / ( d r ) , f = [ d − ( (cid:96) + 1) d ] / ( rd ) , f = [ d + ( (cid:96) + 2) d ] / ( d r ) , f = [( (cid:96) − d − d + 2 n d ] / ( rd ) , f =[(2 (cid:96) + 1) d + 2 n ] / ( rd ) , f = [( (cid:96) + 1)( d − d d − d d + dd d )+( d − (cid:96) ) d +(2 (cid:96)n − d ) d ] / ( d r ) , f = [ d +2(2 (cid:96) +3) dd − d d − (cid:96) +1) d +( (cid:96) +2) d d − (2 (cid:96) +3) dd d ] / ( d r ) , f = [3 d + (5 (cid:96) + 9) dd − d − ( (cid:96) + 2) dd − d ] /d , f = 2 r [ d + ( (cid:96) + 2) d ] /d , f = ( d + dd − d ) /d , f = r [2 d + d ] /d d = n − α r , d = ( (cid:96) + 1) d − α r , d = ( (cid:96)d − α r ) d + 2( (cid:96) − α r d , d = ( (cid:96) − dd − α r d + 2 α r d [( (cid:96) − d + ( (cid:96) − (cid:96) + 2) d − (cid:96) − α r ] , d = (cid:96)d − α r , d = ( (cid:96) − dd − α r d +2( (cid:96) − α r d , d = ( (cid:96) − dd − α r d + 2 α r d [( (cid:96) − d + ( (cid:96) − (cid:96) +1) d − (cid:96) − α r ] g = n + 2( (cid:96) + 1) n + n [3 n + 4( (cid:96) + 1)] α y + [3 n +2( (cid:96) + 1)] α y + α y , g = 4( (cid:96) + 2) n + 4 n (3 (cid:96) + 5) α y +4 n (3 (cid:96) + 4) α y + 4( (cid:96) + 1) α y , g = n ( (cid:96) + 8 (cid:96) + 26 (cid:96) +32 (cid:96) + 13) + (3 (cid:96) + 20 (cid:96) + 50 (cid:96) + 36 (cid:96) + 3) α y + (3 (cid:96) + 10 (cid:96) +7) α y + α y , g = n ( (cid:96) +6 (cid:96) +11 (cid:96) +6)+ n (3 (cid:96) +15 (cid:96) +21 (cid:96) + 5) α y + (3 (cid:96) + 12 (cid:96) + 13 (cid:96) + 4) α y + ( (cid:96) + 1) α y , g = 8[ n (2 (cid:96) +10 (cid:96) +12)+ n (5 (cid:96) +22 (cid:96) +21) α y + n (3 (cid:96) +12 (cid:96) + 9) α y − n α y − α y ] , g = n + 2( (cid:96) + 1) + α y , g = n ( (cid:96) + 4 (cid:96) + 7) + 2( (cid:96) + 3 (cid:96) + 2) α y + α y , g = n [3 n +4( (cid:96) +1)] , g = n [2(5 (cid:96) +9) − n ] , g = 4 n ( (cid:96) +3) , g = n [ n + 2 (cid:96) + 18] , g = 3 (cid:96) + 20 (cid:96) + 50 (cid:96) + 36 (cid:96) + 3 , g = n ( (cid:96) − (cid:96) − (cid:96) − (cid:96) − , g = 8 n ( (cid:96) + (cid:96) − (cid:96) − , g = 16 n ( (cid:96) + 7 (cid:96) + 12) Appendix B. MATLAB Code %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% LINEAR STABILITY CODE FOR PIPE FLOW %% %% M Malik & Martin Skote, %% A linear system for pipe flow stability analysis allowing for %% boundary condition modifications, Computers and Fluids (2019), %% https://doi.org/10.1016/j.compfluid.2019.104267 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%{N − Highest order of Chebyshev polynomials used. N=47 is best foral=0 to 1 and Re < 3000.al − Axial wavenumber alpha as in exp[i(alpha x + n theta) − i omega t]n − Azimuthal wavenumber as in exp[i(alpha x + n theta) − i omega t]Re − Reynolds numbers_bool − Stretching enabler. If s_bool=1, clusters grids towards y=0, elsegrids from − − Stretching parameter. a=2 gives best result for al>3 or Re>6000for small n. a=3 is best for n>6.xi − Gauss Lobatto grid points. ( − − r^2 (square of radial coordinate) (0<=y<=1)D0 − Cheby. poly. matrix (maps from coeff. to Gauss Lobatto(GL) pts) D1 − D4 − Differentiation matrices. Maps Cheb coeffs. on 0<=y_j<=1U,Uy − Mean velocity,U(y), and its derivative with respect to yA,B − Operators in the equation: A q = e Re B qe − Complx frequencies omega's as in exp[i(alpha x + n theta − omega t)]q − Eigenfunctions. To note: phi(y) or psi_1(y) = D0*q(1:N+1,:);Omega(y) or psi_2(y) = D0*q(N+2:2N+2,:) (For definition of phi,Omega, psi_1 and psi_2, see the paper by the authors)%}close all; clear variables% −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− % Control variables are N, al, n, Re, s_bool and a only.N = 200; al = 10; n =7; Re = 2000; s_bool = 1; a=3.0;% −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−− l=abs(n) −
1; N1=N+1; j=(0:N)'; xi=cos(pi*j/N); D1xi=zeros(N1); D2xi = D1xi;D3xi = D1xi; D4xi = D1xi; ons=ones(1,N1); D0=cos((j*(pi*j'))/N);D1xi(:,2:3) = [D0(:,1) 4*D0(:,2)]; D2xi(:,3) = 4*D0(:,1);for ord=3:N%D^kT_n=2nD^{k − − − − − − − − − − − − − − − − − − − − − (2/a)*xiyf.^2; xi_yyy=(4/a)*xiyf.^3; xi_yyyy=(12/a)*xiyf.^4;D1=(xi_y*ons).*D1xi; D2=((xi_y.^2)*ons).*D2xi + (xi_yy*ons).*D1xi;D3=((xi_y.^3)*ons).*D3xi+3*((xi_y.*xi_yy)*ons).*D2xi +(xi_yyy*ons).*D1xi;D4=((xi_y.^4)*ons).*D4xi+6*((xi_y.^2.*xi_yy)*ons).*D3xi + ...((3*xi_yy.^2 + 4*xi_y.*xi_yyy)*ons).*D2xi + (xi_yyyy*ons).*D1xi;elsey=(xi+1)/2; D1 = 2*D1xi; D2 = 4*D2xi; D3 = 8*D3xi; D4 = 16*D4xi;endU = 1 − y; Uy = − ones(N1,1); Uyy = zeros(N1,1);ar2 =al^2*y; d =ar2+n^2; d2 =(l+1)*d − − al*al*D0;A22 = Re*al*(U*ons).*D0 +1i*B11; A11 = Re*al*(U*ons).*B11 − Re*4*al ...*((y.*Uyy)*ons).*D0 +1i*( ((4*y).^2*ons).*D4 + 96*(y*ons).*D3 ...+ ((96 − − − g15 are as defined in above paper)g2 =4*( (l+2)*n^6 + ar2.*((3*l+5)*n^4+ar2.*((3*l+4)*n^2+ar2*(l+1))));g3 =(l^4+8*l^3+26*l^2+32*l+13)*n^2+ar2.*(3*l^4+20*l^3+50*l^2+36*l+3 ...+ ar2.*(3*l^2+10*l+7+ ar2)); g4=(l^3+6*l^2+11*l+6)*n^4+ar2.*( ...(3*l^3+15*l^2+21*l+5)*n^2+ar2.*(3*l^3+12*l^2+13*l+4 + ar2*(l+1)));g5 =8*((2*l^2+10*l+12)*n^6+ar2.*((5*l^2+22*l+21)*n^4+ar2.*( ...(3*l^2+12*l+9)*n^2+ar2.*( − n^2 − ar2)))); g6 =d + 2*(l+1);g7=n^2*(l^2+4*l+7)+ar2.*(2*(l^2+3*l+2)+ar2); g8=n^2*(3*n^2+4*(l+1));g9=n^4*(2*(5*l+9) − n^2); g10=4*n^6*(l+3); g11=n^4*(n^2+2*l+18);g12=3*l^4+20*l^3+50*l^2+36*l+3; g13=n^2*(l^4 − − − − − − − al^2*(g1*ons).*D0 +(g2*ons).*D1 + 4*((y.*d.^3)*ons).*D2;B12 = − − Re*4*al*((y.*d.^3.*Uyy)*ons).*D0 ... − Re*8*al*n^2*((d.^2.*Uy)*ons).*D0 +1i*( ((al^4*g3)*ons).*D0 ... − − − − − − al^2*(g7*ons).*D0 ...+4*((d.*((l+1)*d+n^2))*ons).*D1 + 4*((y.*d.^2)*ons).*D2);end% For n~=0, the system of phi and Omega are coupled.if n~=0, % Enforcing no − slip and regularity Conditions:A = [[D0(1,:);D1(1,:);A11(3:N1,:)] [zeros(2,N1); A12(3:N1,:)]; ...[zeros(1,N1); A21(2:N1,:)] [D0(1,:); A22(2:N1,:)]];B = [[zeros(2,N1); B11(3:N1,:)] [zeros(2,N1); B12(3:N1,:)]; ...[zeros(1,N1); B21(2:N1,:)] [zeros(1,N1); B22(2:N1,:)]];Breg_phi= − al^4*g8*D0(N1,:) +al^2*g9*D1(N1,:)+g10*D2(N1,:);Breg_omg= − − − Re*4*al*l*n^6*D1(N1,:) ...+1i*( al^6*g12*D0(N1,:)+al^4*g13*D1(N1,:) − al^2*g14*D2(N1,:) ...+g15*D3(N1,:)); Areg_omg=Re*al*Breg_omg+Re*2*al^2*n^5*D0(N1,:) ...+1i*( − − − − slip conditions at top rows of the system [Aq =e Re Bq]:tmp1=A(3:N1,:); A(3,:) =A(N1+1,:); A=[A(1:3,:); tmp1; A(N1+2:end,:)];tmp1=B(3:N1,:); B(3,:) =B(N1+1,:); B=[B(1:3,:); tmp1; B(N1+2:end,:)];%Gather chosen 3 unknowns of q to be eliminated at first 3 locs. of q:tmp1 = A(:,3:N1); A(:,3) =A(:,N1+1); A = [A(:,1:3) tmp1 A(:,N1+2:end)];tmp1 = B(:,3:N1); B(:,3) =B(:,N1+1); B = [B(:,1:3) tmp1 B(:,N1+2:end)];%Remove singularities due to no − slip in [Aq = e Re Bq] by elimination:for kk = 1:3, vec = (A(kk,:)/A(kk,kk)); jj = kk+1:2*N1;A(jj,:) = A(jj,:) − A(jj,kk)*vec; B(jj,:) = B(jj,:) − B(jj,kk)*vec;end % Gauss elimination ended; Now solve the reduced system:L = (B(4:end,4:end)\A(4:end,4:end))/Re; [q,e] = eig(L); e =diag(e);ok=((imag(e)> − − − − imag(e));e=e(I);q=q(:,I);tmp1=q(1:N − − A(1:3,1:3)\ (A(1:3,4:end)*q);%Elim'ted coefs.q(1:N1,:) = [q_temp(1:2,:); tmp1];q(1+N1:2*N1,:) = [q_temp(3,:); tmp2];elseif ((n==0) && (al~=0)) %Enforce Conditions for this and the other case:A11=[D0(1,:);D1(1,:);A11(3:N1,:)]; B11 = [zeros(2,N1); B11(3:N1,:)];A22 = [D0(1,:); A22(2:N1,:)]; B22 = [zeros(1,N1); B22(2:N1,:)];A11(N,:)=al*Re*(12*D2(N1,:) − (al^2+8)*D1(N1,:)+al^2*D0(N1,:))...+1i*(192*D3(N1,:) − − al^2*D1(N1,:);else A11=[D0(1,:);A11(2:N1,:)]; B11 = [zeros(1,N1); B11(2:N1,:)];A22 = [D0(1,:); A22(2:N1,:)]; B22 = [zeros(1,N1); B22(2:N1,:)];end% For n==0, the systems are decoupled. So solve for q1 = [phi or psi1] and% q2 = [Omega or psi2] individually for accuracy and speed.if n==0, M =2*(al~=0)+1*(al==0);%M= 2 or 1; M is No. of zero − rows in B_11 or kk = 1:M, vec = (A11(kk,:)/A11(kk,kk));%Gauss eliminationjj = kk+1:N1; A11(jj,:) = A11(jj,:) − A11(jj,kk)*vec;B11(jj,:) = B11(jj,:) − B11(jj,kk)*vec;endL1 = (B11(M+1:end,M+1:end)\A11(M+1:end,M+1:end))/Re; [q1,e1] = eig(L1);e1 =diag(e1);ok =((imag(e1)> − − − − imag(e1));e1=e1(I);q1=q1(:,I);q1=[ − A11(1:M,1:M)\ (A11(1:M,M+1:end)*q1); q1];vec = (A22(1,:)/A22(1,1)); A22(2:N1,:) = A22(2:N1,:) − A22(2:N1,1)*vec;B22(2:N1,:) = B22(2:N1,:) − B22(2:N1,1)*vec;L2 = (B22(2:end,2:end)\A22(2:end,2:end))/Re; [q2,e2] = eig(L2);e2 =diag(e2); ok =((imag(e2)> − − − imag(e2));e2=e2(I);q2=q2(:,I);q2=[ − A22(1,1)\ (A22(1,2:end)*q2); q2];e = [e1; e2]; [eimag,I]=sort( − imag(e));e=e(I);q = [[q1; zeros(N1,size(q1,2))] [zeros(N1,size(q2,2)); q2]]; q=q(:,I);end References [1] T. Itano, S. Toh, The dynamics of bursting process in wall turbulence,Journal of the Physical Society of Japan 70 (3) (2001) 703–716 (2001).[2] T. M. Schneider, B. Eckhardt, J. A. Yorke, Turbulence transition and theedge of chaos in pipe flow, Physical review letters 99 (3) (2007) 034502(2007).[3] N. B. Budanur, B. Hof, Complexity of the laminar-turbulent boundary inpipe flow, Physical Review Fluids 3 (5) (2018) 054401 (2018).[4] L. Böberg, U. Brosa, Onset of turbulence in a pipe, Zeitschrift für Natur-forschung A 43 (8-9) (1988) 697–726 (1988).[5] K. M. Butler, B. F. Farrell, Three-dimensional optimal perturbations inviscous shear flow, Physics of Fluids A: Fluid Dynamics 4 (8) (1992)1637–1650 (1992).[6] L. N. Trefethen, A. E. Trefethen, S. C. Reddy, T. A. Driscoll, Hydrody-namic stability without eigenvalues, Science 261 (5121) (1993) 578–584(1993).[7] P. J. Schmid, D. S. Henningson, Optimal energy density growth in hagen–poiseuille flow, Journal of Fluid Mechanics 277 (1994) 197–225 (1994).[8] P. J. Schmid, Nonmodal stability theory, Annu. Rev. Fluid Mech. 39(2007) 129–162 (2007).[9] M. Gavarini, A. Bottaro, F. Nieuwstadt, The initial stage of transition inpipe flow: role of optimal base-flow distortions, Journal of Fluid Mechan-ics 517 (2004) 131–165 (2004).[10] K. C. Sahu, R. Govindarajan, Stability of flow through a slowly divergingpipe, Journal of Fluid Mechanics 531 (2005) 325–334 (2005).[11] S. Loh, H. Blackburn, Stability of steady flow through an axially corru-gated pipe, Physics of Fluids 23 (11) (2011) 111703 (2011).[12] V. Pruša, On the influence of boundary condition on stability of hagen–poiseuille flow, Computers & Mathematics with Applications 57 (5)(2009) 763–771 (2009).[13] M. Lessen, S. G. Sadler, T.-Y. Liu, Stability of pipe poiseuille flow, ThePhysics of Fluids 11 (7) (1968) 1404–1409 (1968).[14] D. Burridge, P. Drazin, Comments on â ˘AIJstability of pipe poiseuilleflowâ ˘A˙I, The Physics of Fluids 12 (1) (1969) 264–265 (1969).[15] H. Salwen, C. E. Grosch, The stability of poiseuille flow in a pipe ofcircular cross-section, Journal of Fluid Mechanics 54 (1) (1972) 93–112(1972).[16] H. Salwen, F. W. Cotton, C. E. Grosch, Linear stability of poiseuille flowin a circular pipe, Journal of Fluid Mechanics 98 (2) (1980) 273–284(1980).[17] P. Oâ ˘A ´Zsullivan, K. Breuer, Transient growth in circular pipe flow. I. lin-ear disturbances, Physics of fluids 6 (11) (1994) 3643–3651 (1994).[18] C. Thomas, A. Bassom, P. Blennerhassett, The linear stability of oscillat-ing pipe flow, Physics of Fluids 24 (1) (2012) 014106 (2012).[19] V. Priymak, T. Miyazaki, Accurate navier–stokes investigation of tran-sitional and turbulent flows in a circular pipe, Journal of ComputationalPhysics 142 (2) (1998) 370–411 (1998).[20] M. R. Khorrami, M. R. Malik, R. L. Ash, Application of spectral colloca-tion techniques to the stability of swirling flows, Journal of ComputationalPhysics 81 (1) (1989) 206–229 (1989).[21] A. Meseguer, L. Trefethen, A spectral petrov-galerkin formulation forpipe flow I: linear stability and transient growth, Tech. rep., Oxford Uni-versity Computing Laboratory (2000).[22] A. Meseguer, L. N. Trefethen, Linearized pipe flow to reynolds number , Journal of Computational Physics 186 (1) (2003) 178–197 (2003).[23] M. Malik, M. Skote, R. Bouffanais, Growth mechanisms of perturbations in boundary layers over a compliant wall, Physical Review Fluids 3 (1)(2018) 013903 (2018).[24] T. Ellingsen, E. Palm, Stability of linear flow, The Physics of Fluids 18 (4)(1975) 487–488 (1975).[25] O. Y. Zikanov, On the instability of pipe poiseuille flow, Physics of Fluids8 (11) (1996) 2923–2932 (1996).[26] W. E. Boyce, R. C. DiPrima, Elementary differential equations andboundary value problems., John Wiley & Sons, 2009 (2009).[27] C. Foias, O. Manley, R. Rosa, R. Temam, Navier-Stokes equations andturbulence, Vol. 83, Cambridge University Press, 2001 (2001).[28] B. Rummler, The eigenfunctions of the stokes operator in specialdomains. I, ZAMM-Journal of Applied Mathematics and Mechanic-s/Zeitschrift für Angewandte Mathematik und Mechanik 77 (8) (1997)619–627 (1997).[29] B. Rummler, The eigenfunctions of the stokes operator in specialdomains. II, ZAMM-Journal of Applied Mathematics and Mechanic-s/Zeitschrift für Angewandte Mathematik und Mechanik 77 (9) (1997)669–675 (1997).[30] P. F. Batcho, G. E. Karniadakis, Generalized stokes eigenfunctions: anew trial basis for the solution of incompressible navier-stokes equations,Journal of Computational Physics 115 (1) (1994) 121–146 (1994).[31] L. N. Trefethen, Spectral methods in MATLAB, Vol. 10, Siam, 2000(2000).[32] P. J. Schmid, D. S. Henningson, Stability and transition in shear flows,Vol. 142, Springer Verlag, 2001 (2001).[33] G. H. Golub, C. F. Van Loan, Matrix computations, 3rd Edition, JHUpress, 1996 (1996).[34] L. Bergström, Initial algebraic growth of small angular dependent dis-turbances in pipe poiseuille flow, Studies in Applied Mathematics 87 (1)(1992) 61–79 (1992).[35] A. Hanifi, D. S. Henningson, The compressible inviscid algebraic insta-bility for streamwise independent disturbances, Physics of Fluids 10 (8)(1998) 1784–1786 (1998).[36] M. Malik, J. Dey, M. Alam, Linear stability, transient energy growth,and the role of viscosity stratification in compressible plane couette flow,Physical Review E 77 (3) (2008) 036322 (2008).[37] L. H. Gustavsson, Energy growth of three-dimensional disturbances inplane poiseuille flow, Journal of Fluid Mechanics 224 (1991) 241–260(1991).[38] P. Luchini, Reynolds-number-independent instability of the boundarylayer over a flat surface: optimal perturbations, Journal of Fluid Mechan-ics 404 (2000) 289–309 (2000).[39] L. Bergström, Optimal growth of small disturbances in pipe poiseuilleflow, Physics of Fluids A: Fluid Dynamics 5 (11) (1993) 2710–2720(1993)., Journal of Computational Physics 186 (1) (2003) 178–197 (2003).[23] M. Malik, M. Skote, R. Bouffanais, Growth mechanisms of perturbations in boundary layers over a compliant wall, Physical Review Fluids 3 (1)(2018) 013903 (2018).[24] T. Ellingsen, E. Palm, Stability of linear flow, The Physics of Fluids 18 (4)(1975) 487–488 (1975).[25] O. Y. Zikanov, On the instability of pipe poiseuille flow, Physics of Fluids8 (11) (1996) 2923–2932 (1996).[26] W. E. Boyce, R. C. DiPrima, Elementary differential equations andboundary value problems., John Wiley & Sons, 2009 (2009).[27] C. Foias, O. Manley, R. Rosa, R. Temam, Navier-Stokes equations andturbulence, Vol. 83, Cambridge University Press, 2001 (2001).[28] B. Rummler, The eigenfunctions of the stokes operator in specialdomains. I, ZAMM-Journal of Applied Mathematics and Mechanic-s/Zeitschrift für Angewandte Mathematik und Mechanik 77 (8) (1997)619–627 (1997).[29] B. Rummler, The eigenfunctions of the stokes operator in specialdomains. II, ZAMM-Journal of Applied Mathematics and Mechanic-s/Zeitschrift für Angewandte Mathematik und Mechanik 77 (9) (1997)669–675 (1997).[30] P. F. Batcho, G. E. Karniadakis, Generalized stokes eigenfunctions: anew trial basis for the solution of incompressible navier-stokes equations,Journal of Computational Physics 115 (1) (1994) 121–146 (1994).[31] L. N. Trefethen, Spectral methods in MATLAB, Vol. 10, Siam, 2000(2000).[32] P. J. Schmid, D. S. Henningson, Stability and transition in shear flows,Vol. 142, Springer Verlag, 2001 (2001).[33] G. H. Golub, C. F. Van Loan, Matrix computations, 3rd Edition, JHUpress, 1996 (1996).[34] L. Bergström, Initial algebraic growth of small angular dependent dis-turbances in pipe poiseuille flow, Studies in Applied Mathematics 87 (1)(1992) 61–79 (1992).[35] A. Hanifi, D. S. Henningson, The compressible inviscid algebraic insta-bility for streamwise independent disturbances, Physics of Fluids 10 (8)(1998) 1784–1786 (1998).[36] M. Malik, J. Dey, M. Alam, Linear stability, transient energy growth,and the role of viscosity stratification in compressible plane couette flow,Physical Review E 77 (3) (2008) 036322 (2008).[37] L. H. Gustavsson, Energy growth of three-dimensional disturbances inplane poiseuille flow, Journal of Fluid Mechanics 224 (1991) 241–260(1991).[38] P. Luchini, Reynolds-number-independent instability of the boundarylayer over a flat surface: optimal perturbations, Journal of Fluid Mechan-ics 404 (2000) 289–309 (2000).[39] L. Bergström, Optimal growth of small disturbances in pipe poiseuilleflow, Physics of Fluids A: Fluid Dynamics 5 (11) (1993) 2710–2720(1993).