Hocine Bentrad

Ahmed Esmael

Chérif Nouar
email: cherif.nouar@univ-lorraine.fr

Alain Lefevre

Noureddine Ait-Messaoudene

Energy growth in Hagen-Poiseuille flow of Herschel-Bulkley fluid

Keywords: Herschel-Bulkley fluid, cylindrical pipe, transient-growth, energy method

Linear stability of Poiseuille flow of Herschel-Bulkley fluid in a cylindrical pipe is studied using modal and non-modal approaches. The first part of the present study thus deals with the classical normal mode approach in which the resulting eigenvalue problem is solved using a Chebyshev collocation method. Within the considered range of parameters, the modal-linear theory predicts that perturbations are dumped exponentially. In the second part, the effect of the rheological behaviour of the fluid on the pseudospectra and the most amplified perturbations is investigated. At very low Herschel-Bulkley number (Hb << 1 ), the optimal perturbation consists of almost streamwise vortices, and the amplification of the kinetic energy is provided by the lift-up mechanism. In contrast, for sufficiently large values of Hb, the optimal perturbation is axisymmetric and the growth of the kinetic energy is provided by the Orr-mechanism. For intermediate values of Hb, the optimal perturbation is oblique. The amplification of such perturbation is due to a synergy between Orr and lift-up mechanisms. In the last part of the study, the maximal value of the Reynolds number, Re cE , below which the perturbation energy decreases monotonically with time is computed for a large range of Hb. Asymptotic behaviors of Re cE for Hb << 1 and Hb >> 1 are established. The influence of the terms arising from the viscosity perturbation is highlighted throughout this study.

Introduction

The first step in the study of the transition to turbulence consists in considering the linearized equations and the starting point is the standard normal mode approach. Existence of unstable modes means that there is an exponential growth of infinitesimal perturbation and the eigenfunction represents the pattern that appear in early stage of transition. The predictions of the normal mode analysis match the laboratory experiments for some flows, notably Taylor Couette flow and Rayleigh-Bénard convection. For other flows, such as, open parallel shear flows, the normal mode approach fails to match the experiments. The failure of the normal mode linear stability analysis to describe the early stage of transition is attributed in part to the non-normal nature of the linearized equations, i.e. their eigenfunctions are not orthogonal. The non-normal nature manifests itself by a transient growth of perturbations and a large receptivity to ambient disturbances. The amplification of the kinetic energy of the perturbation may reach a significant amplitude that can trigger nonlinear mechanisms before its long time decay due to viscous effects [START_REF] Grossmann | The onset of shear flow turbulence[END_REF]. The physical mechanism behind the transient growth is related to the inviscid vortex tilting mechanism in the presence of base flow shear. The Orr-mechanism [START_REF] Orr | The stability or instability of the steady motions of a perfect liquid and of a viscous liquid. part I: A perfect liquid[END_REF] and the lift-up mechanism [START_REF] Ellingsen | Stability of linear flow[END_REF], [START_REF] Landhall | A note on an algebraic instability of inviscid parallel shear flows[END_REF], [START_REF] Buttler | Three optimal perturbations in viscous shear flow[END_REF] are two such commonly identified disturbance growth phenomena in a shear flow. The Orrmechanism involves disturbance field that consists of spanwise-uniform vortices initially tilted against the direction of the base flow. The disturbance is then sheared and tilted. Its kinetic energy is amplified by the base shear via the working of the Reynolds stress. In the lift-up mechanism, streamwise uniform vortices superimposed on a parallel shear flow lift up fluid at low velocity from the wall and push down towards the wall fluid with high velocity. Hence, elongated streaks with high and low velocity are generated. For arbitrary 3D perturbations, both Orr and lift-up mechanisms are operating, but in a hierarchical manner. The Orr mechanism generates a transient normal-wall velocity, which in turn creates by the lift-up mechanism streamwise streaks that decay slowly by viscosity [START_REF] Farrell | Optimal excitation of three-dimensional perturbations in viscous constant shear flows[END_REF]. The optimal perturbation, i.e. the initial perturbation that undergoes the largest amplification of the kinetic energy, is a fundamental concept of the linear stability theory [START_REF] Schmid | Stability and transition in shear flows[END_REF]. In simple shear flows, such as Couette, Poiseuille or Hagen-Poiseuille flows, the optimal perturbation takes the form of streamwise or quasi streamwise vortices which evolve, via lift-up effect, into spanwiseperiodic streaks (regions of positive and negative streamwise velocity) elongated along the main flow direction. A review on the lift-up process in shear flows can be found in [START_REF] Brandt | The lift-up effect: The linear mechanism behind transition and turbulence in shear flows[END_REF] Comparatively to the Newtonian case, very few studies have been devoted to the transient growth in non-Newtonian fluid flows. This is perhaps not sur-prising, given the inherent complexities involved. Non-Newtonian fluids can be divided in three broad groups: (i) time-independent fluids, for which the viscosity is solely dependent on the instantaneous second invariant of the strain-rate tensor; (ii) viscoelastic fluids, for which stresses depend on the flow history and (iii) time-dependent fluids, for which the strain-rate is a function of both the magnitude and the duration of the applied stress.

Transient growth analysis of inertialess Couette and Poiseuille flows of Oldroyd-B models was studied by Jovanovic and Kumar [START_REF] Jovanovic | Transient growth without inertia[END_REF], [START_REF] Jovanovic | Nonmodal amplification of stochastic disturbances in strongly elastic channel flows[END_REF]. It is shown that the presence of elasticity can produce significant energy amplification. The most amplified disturbances tend to be elongated in the streamwise direction. The mechanism responsible for the energy amplification is related to the stretching of the polymer molecules by the background shear which results in a lift-up of the disturbances. In the case of inertia dominated channel flow of viscoelastic fluids modeled by Oldroyd-B and FENE-P closures, Zhang et al. [START_REF] Zhang | Linear stability analysis of channel flow of viscoelastic Oldroyd-B and FENE-P fluids[END_REF] showed that the ratio between the polymer relaxation time and the characteristic time of instability plays a key role on the energy amplification.

Studies on non-modal transient growth in viscosity-stratified fluid flows were recently reviewed by Govindarajan and Sahu [START_REF] Govindarajan | Instabilities in viscosity-stratified flows[END_REF]. In the present study, we focus on purely viscous (inelastic) shear-thinning fluid flows. The mean features of such fluid flows are: (i) an increase of the wall shear-rate with increasing shear-thinning effects, (ii) a viscosity stratification and (iii) an anisotropy of the shear-stress tensor perturbation which arises from the viscosity perturbation. Chikkadi et al. [START_REF] Chikkadi | Preventing transition to turbulence: A viscosity stratification does not always help[END_REF] considered the case of plane Poiseuille flow of a Carreau fluid, without taking into account for the anisotropic nature of the shear-stress disturbance tensor. Rather, unexpectedly, the authors found that the transient growth is relatively unaffected by the viscosity gradient. This problem has been revisited by Nouar et al. [START_REF] Nouar | Delaying transition to turbulence in channel flow: revisiting the stability of shear-thinning fluids[END_REF] , and obtained a substantial increase of the transient growth when the viscosity perturbation is taken into account. The optimal perturbation is analogous to that of Newtonian fluids, i.e. longitudinal streamwise vortices which transform into streaks by the lift-up mechanism. Similar results are also obtained in the case of pipe flow of shearthinning fluids [START_REF] Liu | Nonmodal stability in hagen poiseuille flow of a shearthinning fluid[END_REF]. In the case of the plane Couette flow, the shear-thinning effects reduce only to the anisotropy of the shear-stress tensor perturbation. A significant increase of the transient growth is observed [START_REF] Liu | Non-modal instability in plane Couette flow of a power-law fluid[END_REF]. However, the optimal perturbation remains similar to that for a Newtonian fluid.

Studies dealing with the transient growth in yield stress shear-thinning fluids are very limited. Non-modal stability in plane Poiseuille flow of a yield stress fluid was studied by Nouar et al. [START_REF] Nouar | Modal and non-modal linear stability of the plane Bingham-Poiseuille flow[END_REF]. The rheological behavior of the fluid is described by the Bingham model. Comparatively to pure shear-thinning fluids a fundamental additional parameter intervenes in the problem: the Bingham or the generalized Bingham number, B, which is the ratio of the yield to viscous stresses. This additional parameter will modify the viscosity stratification, the wall shear rate, the anisotropy of the perturbation stress tensor and the geometry of the yielded zone, where the exchange between the base flow and the perturbation occurs. Nouar et al. [START_REF] Nouar | Modal and non-modal linear stability of the plane Bingham-Poiseuille flow[END_REF] found that a large amplification of the kinetic energy of the perturbation can be obtained for such fluid flows. At small B, the optimal perturbation is in the form of streamwise uniform vortices, whereas for large B, the optimal transient growth occurs for an oblique wave. Similar results were also obtained by Liu & Liu [START_REF] Liu | Non-modal stability in Hagen-Poiseuille flow of a Bingham fluid[END_REF] in the case of Hagen-Poiseuille flow of a Bingham fluid. However, in these studies the physical mechanisms associated with the obliquity of the optimal perturbation are not provided. The influence of the viscosity stratification and that of the anisotropy of the perturbation shear-stress tensor is not clarified.

The objective of the present work is to provide a comprehensive description of the temporal linear stability analysis of Hagen-Poiseuille flow of a shearthinning yield stress fluid. The rheological behavior of the fluid is described by the Herschel-Bulkley which is used generally in the rheological analysis of yield stress fluids. The goal is to clarify the influence of the viscosity stratification and that of the anisotropy of the stress-tensor perturbation. The structure of the paper is as follows. The governing equations and the base flows details are presented in Sec. 2 and 3 respectively. The linear stability problem is formulated in Sec. 4. The concept of tangent viscosity is introduced. The results of the modal approach and the behavior of the least stable mode are discussed in Sec. 5. The eigenvalue spectra are quite similar to those obtained in Newtonian fluids. Nevertheless, we show that for axisymmetric disturbances, the influence of the yield stress on eigenvalue spectra is very weak. An interpretation based on the tangent viscosity is given. The nonnormality of the linear operators is characterized by the ǫ-pseudospectra and the numerical range in Sec. 6. The results of the energy growth, the optimal perturbations structures and the energy growth scaling are also provided in Sec. 6. The influence of tangent viscosity on the transient growth mechanism is highlighted. The energy stability analysis is presented in Sec. 7. Finally, the summary and conclusions are provided in Sec. 8.

Governing equations

The flow of an incompressible shear-thinning fluid with a yield stress τ0 in a circular duct of radius R is considered. A constant pressure-gradient ∂ P ∂ ẑ is imposed in the axial direction e z . The governing equations in dimensionless form are:

∇ • U = 0, (1) ∂U ∂t + U • ∇U = -∇P + ∇ • τ , ( 2
)
where U is the velocity, P the pressure and τ the deviatoric of the extra stress tensor. The velocity vector U is of the form U = Ue r + V e θ + W e z , where U, V, W are the velocity components and e r , e θ , e z are unit vectors in the radial r, circumferential θ and axial z directions respectively. The above equations are non-dimensionalized using the radius of the pipe R as the length scale, maximum velocity Ŵ0 of the basic flow as velocity scale and ρ Ŵ 2 0 as stress and pressure scale. The dimensional quantities are denoted with a hat symbol. The rheological behavior of the fluid is assumed to be described by the Herschel-Bulkley model, which is more realistic than Bingham law. According to [START_REF] Balmforth | Yielding to stress: Recent developments in viscoplsatic fluid mechanics[END_REF], the material structure that resists deformation and leads to the yield stress is typically not completely destroyed at τ = τ0 . The structure persists post-yield and renders the viscosity shear-rate dependent. For several viscoplastic fluids, the corresponding flow curve is well fitted by the Herschel-Bulkley model. Using the Von-Mises criterion, the scaled constitutive equations are

τ = 1 Re µ γ ⇐⇒ τ II > Hb Re , (3) γII
= 0 ⇐⇒ τ II ≤ Hb Re , (4)
where γII and τ II are the second invariant of the strain rate tensor γ and of the deviatoric of the extras stress tensor τ respectively, given by

γII = 1 2 γij γij 1/2 , τ II = 1 2 τij τ ij 1/2 (5)
and γij = U i,j + U j,i . The dimensionless effective viscosity µ is defined by

µ = Hb γII + ( γII ) n-1 . (6)
The Herschel-Bulkley Hb and the Reynolds Re numbers derived from the non-dimensional governing equations, are defined using a generalized viscosity μgen = K Ŵ0 / R n-1 . They are given by

Hb = τ0 K[ Ŵ0 / R] n , Re = ρ Rn Ŵ 2-n 0 K , ( 7
)
where τ0 is the yield stress, K the consistency and n the flow behavior index.

0 = - dP dz + 1 r d dr (rτ rz ) , (8)
where τ rz,b is given by:

τ rz,b = 1 Re   sgn dW b dr Hb + dW b dr n-1 dW b dr   ⇔ τ rz,b > Hb Re , (9)
dW b dr = 0 ⇔ τ rz,b ≤ Hb Re . (10)
Here, sgn mean sign of the argument. Integration of (8) combined with [START_REF] Jovanovic | Transient growth without inertia[END_REF][START_REF] Jovanovic | Nonmodal amplification of stochastic disturbances in strongly elastic channel flows[END_REF] and the non-slip boundary condition at the wall gives:

P b (z) = -2 Hb Re r 0 z + const, (11)
W b =      1 ; 0 ≤ r ≤ r 0 1 -r-r 0 1-r 0 n+1 n ; r 0 ≤ r ≤ 1, (12)
where r 0 is the dimensionless radius of the plug zone. The basic flow (Fig. 1) is characterized by a central region of radius r 0 moving as a rigid solid, i.e. a plug zone which is surrounded by a yielded region where the viscosity µ b varies nonlinearly with the shear rate: The viscosity stratification is more pronounced as r 0 is higher or n is lower, as it is shown in Fig. 2. The wall shear-rate, dW b /dr at r = 1, increases as the width of the yielded zone is reduced or as the shear-thinning index decreases.

µ b = Hb |dW b /dr| + |dW b /dr| n-1 . (13)
It is worthy to note that the yield surface r = r 0 is not a material surface. The variation of r 0 as function of Hb and n is determined on one hand from the integration of ( 8) and on the other hand from the constitutive equation:

|τ w | = Hb r 0 Re and |τ w | = 1 Re Hb + n + 1 n n 1 (1 -r 0 ) 1+n (14)
From [START_REF] Nouar | Delaying transition to turbulence in channel flow: revisiting the stability of shear-thinning fluids[END_REF] it can be shown that the radius of the yield surface is solution of the following equation

Hb (1 -r 0 ) n+1 - n + 1 n n r 0 = 0 , (15)
For low and large Hb, the following asymptotic behaviors are derived:

r 0 ≈ n n + 1 n Hb -(n + 1) n n + 1 2n Hb 2 as Hb → 0, ( 16
)
r 0 ≈ 1 - n + 1 n n n+1
1 Hb fluid with respect to infinitesimal perturbations. This influence arises from four different effects: (i) variation of the yielded zone, where the exchange of energy between the perturbation and the base flow occurs; (ii) variation of the axial velocity gradient dW b /dr; (iii) viscosity stratification and (iv) nonlinear variation of the viscosity with the shear rate. The contribution of each of these effects on the flow stability will be highlighted whenever possible.

µ b /µ b (r = 1) (1) (2) (3)
µ b /µ b (r = 1) (4) (3) (2) (1)

Linear stability analysis

An infinitesimal perturbation (ǫu ′ , ǫp ′ ), with ǫ ≪ 1, is imposed on the base flow (U b , P b ). The perturbed flow is given by:

(U b + ǫu ′ , P b + ǫ p ′ ) = (ǫ u ′ , ǫ v ′ , W b + ǫ w ′ , P b + ǫp ′ ) . (18)
Wherever the yield stress is exceeded, i.e. τ > Hb/Re, the effective viscosity of the perturbed flow is expanded about the base flow as:

µ (U b + ǫ u ′ ) = µ b + ǫ µ ′ + ..., (19)
where

µ ′ = ∂µ ∂ γij b γij (u ′ ) . ( 20
)
The deviatoric stresses in the disturbed flow can also be written as:

τ ij (U b + ǫ u ′ ) = τ ij (U b ) + ǫτ ′ ij + ..., (21)
with

τ ′ ij = µ b γij (u ′ ) + µ ′ γrz (U b ) (22)
For one dimensional shear flow, with W b (r) in the streamwise direction z, the components of the deviatoric stress perturbation read:

τ ′ ij = 1 Re µ b γij (u ′ ) if ij = rz, zr, (23)
τ ′ rz = 1 Re µ b + γrz (U b ) ∂µ ∂ γrz b γrz (u ′ ) = 1 Re µ t γrz (u ′ ) . (24)
In Eq. ( 24), µ t is the tangent viscosity [START_REF] Nouar | Delaying transition to turbulence in channel flow: revisiting the stability of shear-thinning fluids[END_REF], defined by:

µ t = (∂τ rz /∂ γrz ) b = n |DW b | n-1 . ( 25
)
It is worthy to note that: (i) µ t is independent of Hb, i.e. the disturbance "will not feel" the yield stress, and (ii) the deviatoric stress associated to the perturbation is anisotropic. The second invariant of the deviatoric of the stress tensor is linearly disturbed:

|τ (U b + ǫu ′ ) -τ (U b + ǫu)| = O (ǫ)
. Therefore, the yield surface position r y is linearly disturbed from its position r 0 :

r y = r 0 + ǫ h (θ, z, t) . ( 26
)
The disturbance is assumed periodic in the z and θ-directions:

(u ′ , v ′ , w ′ , p ′ , h ′ ) = [u (r, t) , v (r, t) , w (r, t) , p (r, t) , h (t)] exp i (αz + mθ), (27)
where α and m are the axial and azimuthal wave numbers respectively. After some algebra, it can be shown that the linearization of continuity and momentum equations around the base flow leads to the following initial value problem:

D (ru) + i [m v + α r w] = 0, ( 28
)
∂u ∂t = -iαW b u -Dp + 1 Re µ b ∆u - 2im r 2 v - u r 2 + 2 Re Dµ Du + iα Re (µ t -µ b ) (Dw + iαu) , (29)
∂v ∂t = -iαW b v - imp r + 1 Re µ b ∆v + 2im r 2 u - v r 2 + 1 Re Dµ b Dv + i m r u - v r , (30)
∂w ∂t = -iαW b w -DW b u -iα p + 1 Re µ b ∆w + 1 Re Dµ b (Dw + iαu) + 1 Re 1 r D [r (µ t -µ b ) (Dw + i α u)] , (31)
with ∆ ≡ D 2 + 1 r D - m 2 r 2 -α 2 and D ≡ d dr .
The boundary conditions for the perturbation velocity are obtained from the non-slip and non-penetration at the wall (r = 1):

u (1) = v (1) = w (1) = 0. ( 32
)
The continuity of stress at the yield surface requires

γij (U b + ǫu ′ ) = 0 at r = r y (33)
Expanding and linearizing about r = r 0 give:

u (r 0 ) = v (r 0 ) = w (r 0 ) = 0 , (34)
Du (r 0 ) = Dv (r 0 ) = 0; Dw (r 0 ) =        0 if n < 1 h 2 (1 -r 0 ) 2 if n = 1 . (35)
We end up with nine boundary conditions: three at the wall and six at the yield surface. The system seems overdetermined. Actually, the conditions Du(r 0 ) = Dv(r 0 ) = 0 and Dw(r 0 ) = 0 when n < 1 are needed to overcome the singularity in the viscosity at r = r 0 . When n = 1, the boundary condition ( 35) is a condition for h not for w. The system of equations ( 28) -( 31) may be expressed in terms of u and v if α = 0 or in terms of u and w if m = 0. One has to note that at the leading order, the unyielded zone is not disturbed and behaves as a rigid solid from kinematic point of view.

Long-time behaviour of the disturbance: eigenvalue problem

When the long time behavior is sought, the disturbance is assumed to behave exponentially as:

Ψ(r, t) = ψ(r) exp (-i C t), (36)
where Ψ stands for (u, v) T or (u, w) T and C = C r + iC i is the complex frequency. The phase speed of the perturbation is given by C r /α and the growth rate by C i . The initial value problem is transformed into the following generalized eigenvalue with C as eigenvalue

L uv (u, v) T = CM uv (u, v) T or L uw (u, w) T = CM uw (u, w) T (37)
depending on whether the (u, v) (if α = 0) or (u, w) (if m = 0) formulation is used. Note that the eigenvalue problem ( 37) is only defined on the yielded zone r ∈ [r 0 , 1]. With increasing the Herschel-Bulkley number, the width (1 -r 0 ) of the sheared fluid zone decreases. To take into account this geometrical effect, the following reduced parameters are introduced:

r = r(1 -r 0 ) + r 0 , z = z(1 -r 0 ), t = t(1 -r 0 ), (38) α
= α (1 -r 0 ) , C = C (1 -r 0 ) , ( 39
) Re = Re (1 -r 0 ) n , Hb = Hb (1 -r 0 ) n . ( 40
)
In terms of r, the basic velocity profile

W b = 1 - r 1+n n ; r ∈ [0, 1] (41)
is artificially independent of the Herschel-Bulkley number. The differential eigenvalue problem [START_REF] Esmael | Transitional flow of a nonnewtonian fluid in a pipe: Experimental evidence of weak turbulence induced by shear-thinning behavior[END_REF], written in terms of tilde variables, with the boundary conditions, is discretized using the Chebyshev spectral collocation method [START_REF] Schmid | Stability and transition in shear flows[END_REF] at (N + 1) Gauss-Lobatto collocation points. The resulting generalized eigenvalue problem is solved using the QZ-algorithm available in Matlab software package. Spectra with increasing collocation points were compared to determine the adequate number of Chebyshev polynomials (N + 1). It was found that with N = 120, the first 15 eigenvalues taken by increasing Ci were resolved accurately within four digits (invariant four digits with increasing N) for almost all the situations considered in this paper.

Case of a one-dimensional perturbation

In the case of a one-dimensional perturbation, i.e. α = m = 0, the linear inertial terms which ensure the exchange of energy between the base flow and the perturbation vanish. The initial value problem ( 28)-( 31) with the associated boundary conditions ( 32)-( 35) reduces to

u = 0, p(r) = constant, ( 42
) ∂v ∂ t = 1 Re 1 η 2 D η 3 μb D v η , (43) ∂w
∂ t = 1 Re 1 η D η μt Dw , (44)
where

η = r + r 0 1 -r 0 , D = d dr , μb = Hb | DW b | + | DW b | n-1 , μt = n D|W b | n-1 . ( 45
)
The boundary conditions are:

v(r = 1) = w(r = 1) = v(r = 0) = w(r = 0) = 0; (46)
Multiplying ( 30) by v and ( 31) by w and then integrating from 0 to 1, the following expressions are obtained: Looking for a solution of ( 43), [START_REF] Pralits | Weakly nonlinear optimal perturbations[END_REF] as exp -i C t (v, ŵ), it can be shown straightforwardly that the eigenvalues are purely imaginary. Fig. 3(a) is an example of spectrum obtained at Re = 5000 for Herschel-Bulkley fluid with r 0 = 0.5 and n = 0.5. The maximum growth rate Ci,max behaves as Re -1 (Fig.

1 2 ∂ ∂ t v 2 r = - 1 Re μb η D v η 2 r (47)

3(b))

. For a fixed r 0 , Ci,max decreases with increasing shear-thinning effects, whereas for fixed n, Ci,max increases with increasing r 0 . Opposite effects can be found by using a Reynolds number defined with the wall shear-viscosity.

The selection of the viscosity scale may be considered a matter of choice, however the conclusion that one reaches by comparing shear-thinning fluids among themselves and against Newtonian fluid can be radically different from one choice to another.

Case of axisymmetric perturbation: m = 0

The (u, v) formulation is used. The differential eigenvalue problem [START_REF] Esmael | Transitional flow of a nonnewtonian fluid in a pipe: Experimental evidence of weak turbulence induced by shear-thinning behavior[END_REF] reduces to two decoupled differential equations: the Orr-Sommerfeld and Squire equations.

-The Orr-Sommerfeld modes are solution of the equation

L os u = -i C ReLu, (49)
with the boundary conditions:

u(1) = Du(1) = u(0) = Du(0) = 0, (50)
where, L ≡ D D+ -α2 , D+ = D+1/η and L os is the Orr-Sommerfeld operator:

L os ≡ -iα ReW b L -iα Re DW b η -D2 W b + μb L 2 + D μb 2 DL + 1 η L + D2 μb L + 2 α2 + L + 2 α2 (μ t -μb ) L + 2 α2 . ( 51
)
The eigenvalues spectrum of C arising from the discretized operator -(i/Re)L -1 L os exhibits a similar three-branches as for the case of the Newtonian plane Poiseuille flow.

In order to highlight separately the influence of r 0 , i.e, the modification of the flow geometry on the eigenvalue spectra, we have fixed n = 1 and we have canceled artificially Hb terms in μb expression. The result is illustrated in Fig. 4(a). For a fixed Re, the maximum growth rate increases with decreasing the width of the yielded zone. This result may be related qualitatively with that of Newtonian fluid flow stability in an annular space. In such case, Mott and Joseph [START_REF] Mott | Stability of parallel flow between concentric cylinders[END_REF] have shown that the critical Reynolds number decreases as the ratio of the inner cylinder to that of the outer one increases and approaches 1. In Fig. 4(b), the eigenvalue spectrum obtained for Re = 5×10 3 , α = 1, m = 0, n = 0.5, r 0 = 0.8 and Hb = 6.9282 is compared with that obtained by neglecting the Hb terms in the Orr-Sommerfeld equation. As it can be observed, there is practically no effect of Hb terms on wall modes. For an axisymmetric perturbation, the wall modes are more sensitive to μt rather than μb . Indeed, in Eq. ( 51), μt terms involve a fourth derivative D4 v, while μb involve a second derivative D2 v. Hence, the wall modes do not "fill" the yield stress.

The influence of the shear-thinning index on the wall modes is illustrated by

-Ci Re -1/3 (1) (4)
Re 10 For large values of Re, the scaling with Re is the same as for a Newtonian fluid. Similarly, for interface and mean modes, the scaling with Re is the same as for a Newtonian fluid.

-The Squire modes are solution of the equation

αW b v + i μb Re D D+ -α2 + i Dµ b Re D - 1 η v = Cv. ( 52
) with v(0) = v(1) = 0. (53)
Multiplying (52) by v * , where the star designates the complex conjugate, and integrating between the yield surface and the wall gives:

Ci |v| 2 r = - μb Re   α2 |v| 2 + Dv - v η 2   r , (54)
It is clear that the Squire modes are always damped.

5.3

Case of streamwise homogeneous perturbation: α = 0 and m = 0

The (u, w) formulation is used. On setting α = 0, the eigenvalue problem (37) reduces to

i C D+ u 2 + 1 η D+ u - m 2 η 2 u = - μb Re D+ u 2 + 1 η D+ u - m 2 η 2 u 2 -D μb Re 1 η D D+ u + 2 D3 u + 4 D2 η u - 2m 2 η 2 Du + m 2 η 3 u - D2 μb Re D D+ u - m 2 η 2 u ( 55
) -i Cw = -DW ℓ u + μb Re 1 η D η Dw - m 2 η 2 w + D μb Re Dw + 1 η D η (μ t -μb ) Dw (56)
The set of eigenvalues of the system (55) and ( 56) can be divided into two eigenmodes classes. The first one corresponds to the set of the eigenmodes of (55) with a particular solution of (56). The second class represents the set of eigenmodes of (56) with u = 0. It can be shown straightforwardly that the eigenmodes are always damped. They are pure imaginary, i.e. the eigenmodes do not propagate and decrease monotonically with time. An example of the eigenvalues spectrum for α = 0, m = 1, Re = 5000, n = 0.5 and r 0 = 0.5 is shown in Fig. 6(a). The eigenvalues cluster to the origin with increasing Re. For the streamwise case m = 1 and α = 0, the maximum growth rate approaches zero as Re → ∞ almost exactly as Re -1 . Figure 6 index. The numerical results indicate that Ci,max decreases with decreasing n.

Remark: Comparison with the situation of annular Poiseuille flow with sliding inner cylinder In the linear theory, the unyielded zone behaves as a rigid solid moving axially with a constant velocity. The Hagen-Poiseuille flow of a yield stress fluid can be viewed a priori as a combined axial Poiseuille Couette flow of a shear-thinning fluid in a cylindrical annulus. To our knowledge, this problem has not been considered before in the literature. For a Newtonian fluid, Sadeghi and Higgins [START_REF] Sadeghi | Stability of sliding Couette-Poiseuille flow in an annulus subject to axisymmetric and asymmetric disturbances[END_REF] studied the stability of sliding Couette-Poiseuille flow in an annulus to both axisymmetric and asymmetric perturbations. Computations were performed for a radius ratio R 1 /R 2 = 0.5 and 0 ≤ m ≤ 3. Preziosi and Rosso [START_REF] Preziosi | Stability of a viscous liquid between sliding pipes[END_REF] studied the linear stability of a Newtonian fluid between sliding pipes. The linearized disturbance equations were solved numerically for a radius ratio R 1 /R 2 ≥ 0.1, α ≤ 10 and m ≤ 5 and a Reynolds number, based on the axial velocity of the moving cylinder, less than 10 4 . For this range of parameters, the authors did not observe any instability. Gitler [START_REF] Ph | Stability of axial Poiseuille-Couette flow between concentric cylinders[END_REF], using a long-wave version of the axisymmetric Orr-Sommerfeld equation (proposed initially by Cowley and Smith [START_REF] Cowley | On the stability of Poiseuille-Couette flow: a bifurcation from infinity[END_REF]), showed that for a radius ratio R 1 /R 2 < 0.1415, the Couette flow is linearly unstable from a finite Reynolds number. Based on Gitler's results, one might expect that Hagen-Poiseuille flow of a yield stress fluid would be unstable at a finite critical Reynolds number at least for r 0 ≤ 0.14. A long wave-approximation of the Orr-Sommerfeld equation for Herschel-Bulkley was established (see Appendix A) and solved numerically for two values of r 0 : 0.1 and 0.05. We have not found any instability. This is probably due to the fact that the velocity profile can not be written as the sum of Couette and Poiseuille flow of a yield stress fluid. In this situation, either the (u, v) or the (u, w) formulation can be used. The shape of the eigenvalue spectrum is similar to that of a Newtonian fluid, except that for a yield stress fluids, there are two separate vertical branches, one of them is associated to the Squire modes described in §5.2 as shown in figure 7. As an indication, we have represented the eigenfunctions u for the three selected eigenmodes. As expected, for the wall eigenmode, variations of the axial velocity occur mainly near the wall, while for the interfacial eigenmode, variations are observed near the yield surface (interface). The mean mode is mainly characterized by oscillations in the associated eigenfunctions.

Finally for the range of the rheological parameters considered in this study, all the eigenmodes lie in the stable half of the complex plane. Thus, it is conjectured that Poiseuille flow of Herschel-Bulkley fluid in a cylindrical pipe is linearly stable with respect to infinitesimal perturbations.

6 Short time behavior: pseudospectra, transient growth and optimal perturbation

Nonnormality: pseudospectra and numerical range

The eigenvalues describe the time asymptotic behavior of disturbances. Thus for Hagen-Poiseuile flow of yield stress fluids, a perturbation introduced at t = 0 decays to zero for large times. The stability analysis based on eigenvalues is not sufficient to describe the temporal behavior of the disturbance at all times because of the non-normality of the linear operator. This mathematical property means that there is a potential for extraction of energy from the base flow by a subspace of perturbation leading to transient growth, despite the absence of an exponential instability.

The non-normality of the linear stability operator L ≡ M -1 uw(uv) L uw(uv) , Eq. ( 37), is characterized by using ǫ-pseudospectrum and numerical range tools (Trefethen et al. [START_REF] Trefethen | Hydrodynamic stability without eigenvalues[END_REF], Reddy et al. [START_REF] Reddy | Pseudospectra of the Orr-Sommerfeld operator[END_REF], Schmid & Henningson [START_REF] Schmid | Optimal energy density growth in Hagen-Poiseuille flow[END_REF]). The computation needs to define an inner product and a norm. For a velocity disturbance vector-function q = (u, v) T , the scalar product based on the energy density is defined as

(q 1 , q 2 ) E = 1 r 0 q H
2 Qq 1 rdr with (57)

q H 2 Qq 1 = u 1 u * 2 + 1 α 2 D + u 1 D + u * 2 + 1 + m 2 α 2 r v 1 v * 2 + im α 2 r v 1 D + u * 2 -D + u 1 v * 2 ,
where the superscript 'H' means transpose conjugate. The associated energy norm is given by

||q|| E = (q, q) E = 1 r 0 q H Qqrdr. (58)
Q is a positive definite matrix. It can be decomposed into Q = F H F using a Cholesky decomposition. The energy norm of a perturbation is equivalent to the standard (Eucledian)L2-norm of the vector F q.

Let F be the discrete representation of L. The ǫ-pseudospectrum of F is defined as the set of complex numbers z for which

|| (zI -F ) -1 || E ≥ ǫ -1 .
It is usually displayed graphically with contours of the norm of the resolvent (zI -F ) -1 for various values of ǫ. Pseudospectra can also be defined as the set of complex numbers z which are eigenvalues of F + ∆F for some perturbation matrix ∆F with ||∆F || E ≤ ǫ. The more non-normal the linear operator L, the greater the potential for a disturbance operator ∆F to affect the eigenvalues.

A third equivalent definition of the ǫ-pseudospectrum, closer to computation, involves the singular value decomposition. It is the set of complex numbers z for which σ min zI -F F F -1 ≤ ǫ, where σ min is the smallest singular value. The numerical range of F is the set of complex numbers (F q, q) E / (q, q) E . For a normal operator, the numerical range is the convex hull of the spectrum. Therefore, for a stable normal operator, the numerical range is always in the stable half plane. However, it can extend significantly to even protrude into the unstable half plane for stable non-normal operators. Its maximum protrusion determines the maximum energy growth rate at t = 0 + ([28], [START_REF] Gallino | The stability of a rising droplet: an inertialess non-modal growth mechanism[END_REF])

Axisymmetric perturbation

In figure 9, spectra, pseudospectra and numerical ranges of the linear stability operator L os (51) are shown for axisymmetric perturbations with α = 1, n = 0.5, Re = 5000 and two values of the plug radius r 0 = 0.01 and 0.8. It is clear that the Orr-Sommerfeld operator is non-normal. At Re = 5000, perturbations with norm ǫ > 0.01 is necessary to protrude the pseudospectra in the unstable half plane. As expected for axisymmetric perturbations, the influence of r 0 or Hb on the pseudospectra is very weak. Indeed, the perturbation is mainly sensitive to the tangent viscosity which is independent of the yield stress.

streamwise homogeneous perturbation

Spectra, pseudospectra and numerical ranges of the linear stability operator (55), (56) for streamwise homogeneous perturbation (longitudinal rolls) are depicted in Figure 10 at Re = 5000, α = 0, m = 1, n = 0.5 and two values of r 0 : 0.01 ( Hb = 0.017) and 0.8 ( Hb = 6.93). On contrast with the former case, the increase of Hb reduces significantly the extension of the pesudospectra in the unstable half plane. For instance, at r 0 = 0.01, a perturbation of norm ǫ ≈ 5.2 × 10 -6 is necessary to reach the unstable half plane, while for r 0 = 0.8, a much more intense perturbation ǫ ≈ 1.0 × 10 -3 is needed.

Oblique perturbation

For an oblique perturbation, the pseudospectra and numerical ranges are qualitatively similar to those obtained for an axisymmetric perturbation. However, the influence of Hb is much more significant. Increasing Hb reduces the maximum protrusion of the numerical range and the extension of the pseudospectra in the unstable half plane.

Norm of the minimal perturbation

As for Newtonian fluids, the norm of the minimal destabilizing perturbation, ǫ min behaves as ( Re -1/2 ) for large Re. Note that for a given Re, ǫ min decreases with decreasing n or r 0 .

Transient growth and optimal disturbances

The transient evolution of perturbations in the linear regime is determined following the methodology described by Schmid and Henningson [START_REF] Schmid | Optimal energy density growth in Hagen-Poiseuille flow[END_REF]. Let g(t) the ratio between the energy norm ||q(t)|| E of the perturbation at time t and its initial norm

||q 0 || E g(t) = ||q(t)|| E ||q 0 || E (59)
For a given Fourier mode, the amplification of the energy at time t maximized over all possible non-zero initial conditions is denoted by

G(t, α, m) = sup q 0 =0 g(t). ( 60
)
The maximum transient energy growth possible over all times is

G max (α, m) = sup t G(t, α, m). ( 61
)
The maximum of G max for all the pairs (α, m) is

G opt = sup α,m G(α, m), (62)
which is reached by the optimal perturbation at time t opt . Unlike the modal approach, here the growth of disturbances occurs over relatively short initial time and is related to an inviscid mechanism. The viscosity acts to moderate the amplification.

Transient growth for a given Fourier mode

-Case of oblique and streamwise homogeneous perturbations Figure 12 shows the kinetic energy amplification G( t) as a function of time at Re = 3000 and three values of the dimensionless yield stress: Hb = 0.02 (r 0 = 0.01), Hb = 0.22 (r 0 = 0.1), Hb = 0.5 (r 0 = 0.2). The flow behavior index n is fixed to n = 1. Two different initial conditions are considered: longitudinal rolls with α = 0, m = 1 (Fig. 12a) and oblique perturbation α = 1, m = 1 (Fig. 12b). With increasing Hb, the viscosity increases thus reducing the viscous diffusion time and therefore t max and G max . At very low Hb, the transient growth remains significantly lower than that obtained for a Newtonian fluid, particularly for the case of α = 0, m = 1. The unmatched transient growth at the limit Hb = 0 is explained as follows. In the framework of the linear stability analysis, it is assumed that the perturbation is infinitesimal with respect to all scales of the base flow. If ǫ is a characteristic scale of the perturbation size, therefore ǫ = o(r 0 ). In other words, the yield surface is linearly disturbed and the plug zone remains intact ( [START_REF] Nouar | Modal and non-modal linear stability of the plane Bingham-Poiseuille flow[END_REF], [START_REF] Metivier | Linear stability involving the Bingham model when the yield stress approaches zero[END_REF]). Hence, when the Hb terms are canceled artificially in the perturbation equations, one recovers the linear stability problem of Couette-Poiseuille flow between two coaxial cylinder where the outer is fixed and the inner is moving with a constant axial velocity W = 1.

-Case of axisymmetric perturbations For axisymmetric perturbations, the influence of r 0 on the transient growth is relatively weak. Indeed, for an axisymmetric perturbation, spectra, pseudospectra and numerical range do not substantially depend on r 0 .

Characteristics of the optimal perturbation and optimal transient growth

We have computed G max at Re = 5000 and for the set of rheological parameters n = 0.3, 0.5, 0.7 and 1, with 0.001 ≤ r 0 ≤ 0.96 (0.0017 ≤ Hb ≤ 48.0). The characteristics of the optimal perturbations, i.e. αopt and m opt corresponding to max α,m (G max ) are displayed in the Figure 13 as a function of the plug zone radius. Three regions can be distinguished depending mainly on r 0 : (i) The first region holds for low values of r 0 , say r 0 0.33, which corresponds to Hb ≤ 0.74 when n = 0.5. In this region, α opt = 0 and the optimal azimuthal wavenumber m opt increases with increasing shear-thinning effects. The optimal perturbation consists of pairs of counter-rotating streamwise vortices, In Figure 14, we have represented the structure of the optimal perturbation at t = 0 and t = topt for r 0 = 0.01 and n = 0.5. These streamwise counter-rotating vortices allow the transfer of energy to the streamwise-velocity component by the lift-up mechanism creating high (+) and low (-) streamwise streaks displayed in Fig. 14(b). Note that the location of the maximum streamwise velocity component approaches the wall with increasing Hb or r 0 . At t = 0, almost all the energy is in the azimuthal (61.52%) and radial (38.46%) components and only a negligible part is in the streamwise component. At optimal time, t = t opt , the kinetic energy is merely concentrated in the axial component. (ii) The second region holds for 0.35 r 0 0.93. In this region, αopt increases gradually until α ≈ 0.5 at r 0 ≈ 0.6, then strongly until a maximum at r 0 ≈ 0.93. With increasing Hb or shear-thinning effects, the width of the zone where the exchange of energy between the base flow and the perturbation occurs, is reduced leading to an increase of m opt . The optimal perturbation is oblique and its wave vector has an angle which increases with increasing r 0 . The time evolution of the axial velocity component, in the (r, z) and (r, θ) planes is shown in Fig. 15. The optimal perturbation is initially oriented in the direction opposite to the mean shear, then aligned with it at the optimal time. As indicated by Farrell and Ioannou [START_REF] Farrell | Optimal excitation of three-dimensional perturbations in viscous constant shear flows[END_REF], the kinetic energy of the perturbation is amplified thanks to Orr-and lift-up mechanisms. An alternative explanation of the transient growth was proposed by Vitoshkin et al. [START_REF] Vitoshkin | On the role of vortex stretchning in energy optimal growth of three-dimensional perturbations on plane parallel shear flows[END_REF]. These authors pointed out the role played by the interplay between the spanwise (here azimuthal) vorticity q and the divergence plane d (in the shear-plane, rz). This is highlighted in the analysis of the energy growth via the Reynolds stress production term. For this, a 2D Helmoltz decomposition is performed:

u = u d + u q with u d = ∇ (ϕ) and u q = -rot (ψe θ ) , (63)
where ϕ and ψ are scalar functions such that the planar divergence d = ∂u/∂r+u/r+∂w/∂z and the azimuthal vorticity q = ∂v/∂r+v/r-(1/r)∂u/∂θ satisfy

d = ∆ H (ϕ) and q = ∆ H (ψ) (64)
where ∆ H ≡ ∂ 2 /∂r 2 + (1/r)∂/∂r + ∂ 2 /∂z 2 is the 2D Laplacian in the shear plane (r, z). The energy growth via the Reynolds stress is

+ + - - (a) (b)
∂ E ∂t = -ℜ γb rz uw * = -ℜ γb rz (u d + u q ) w * d + w * q (65) = -ℜ γb rz u d w * d -ℜ γb rz u q w * q -ℜ γb rz u d w * q -ℜ γb rz u q w * d ,
with E = |u| 2 and (•) = 1 r 0 (•)rdr. Time evolution of the different terms is depicted in Fig. 16(a) for n = 0.5, r 0 = 0.5, Re = 5000. The first term on the right-hand side of Eq. (65) represents the contribution of the planar divergence to the energy growth. This contribution is positive, but relatively small. The second term, i.e. the rotational term -γb rz u q w * q , is the contribution of the 2D Orr-mechanism. Initially, this term is positive and participates to the increase of E. From t = t opt /2, the vortices are tilted with the shear and -γb rz u q w * q is negative: the Orr-mechanism contributes to the decay of the energy. The third term, -γb rz u d w * q is the first mixed q -d term. It is positive and contributes strongly to the energy growth and is able to overwhelm the large negative contribution of the Orr-mechanism, when the vortices align with the base shear. The fourth term, -γb rz u q w * d is the second mixed q -d term. It is small and negative. With increasing r 0 , the contribution of the Orr-mechanism in the amplification of the perturbation kinetic energy becomes more significant as illustrated in Fig. 16(b) and 16(c).

(iii) For very large plug radius, 0.93 ≤ r 0 < 1, the optimal perturbation is axisymmetric. The transient growth of the perturbation energy arises from the Orr-mechanism. It is associated with the tilting of the disturbance into the direction of the mean shear. This tilting phenomenon is illustrated by Fig. 17, where contours of axial velocity of the optimal perturbation are displayed at different times. In contrast with the previous situation (ii), the 2D optimal perturbation is oriented transversally to the main shear at the optimal time.

γb rz u d w * q γb rz u q w * q -γb rz uw *

(1) (2) -γrz u q w * q ;

-γb rz u d w * d (4) (3) (c) (d)
(3) -γrz u d w * q ; (4) -γrz u q w * d

Discussion

• If we set n = 1, and we cancel artificially Hb terms in the linear stability operators, we recover the Couette-Poiseuille flow of a Newtonian fluid. The optimal perturbation consists of longitudinal counter rotating vortices pairs (α = 0) for any r 0 < 1. The lift-up mechanism is the only operating mechanism in the energy amplification. The azimuthal wavenumber increases significantly with increasing r 0 . A slight increase of G opt with increasing r 0 is observed. • The obliquity of the optimal perturbation which may be characterized by the ratio arctang(α/m) is a consequence of the anisotropy terms, i.e. terms arising from the viscosity perturbation. Indeed, if we cancel artificially µ t -µ b terms as well as their derivatives, the obliquity is strongly reduced. For instance at r 0 = 0.9, the ratio α/m is much small (10 times smaller) when the anisotropy terms are neglected. Note that the energy amplification is also significantly reduced and the lift-up is the main mechanism operating. The contribution of the different terms of Eq. (65) in the transient growth at r 0 = 0.9 and n = 0.5 is represented in Fig. 16(d). Comparatively to Fig. 16(c), the contribution of the Orr mechanism is limited.

Scaling laws for G opt and t opt

For given r 0 and n, the dependence of G opt and t opt on Re is studied. It is found that G opt increases with Re. The scaling with Re 2 , Fig 18(a) is recovered for Hershel-Bulkley fluids in the first region where α = 0. It also applies in the second region until r 0 0.6, where the lift-up mechanism is the dominant mechanism in the transient growth (Fig. 16(a)). Similarly, the scaling of t opt with Re is satisfied, Fig. 18(b). Analysis of the numerical results show that they can be fitted as

G opt Re 2 = A n exp (B n r 0 ) and topt Re = A ′ n exp (B ′ n r 0 ) , (66)
where the coefficients A n , A ′ n and the exponents B n , depend on n. They are given in Table 1 7 Energy stability, Euler-Lagrange equations

In this section, the maximum Reynolds number, Re cE , below which, the perturbation kinetic energy decreases monotonically with time is determined. In the yielded zone, Ω = [0; 2 π]×]r 0 ; 1] × [0; 2π/α], the mean kinetic energy, E, of a three dimensional perturbation is defined by:

wDW b = -Dξ + 2Λ µ b △u - 2im r 2 v - u r 2 + 2 dµ dr du dr +2Λ [iα (µ t -µ b ) (Dw + iα)] , (74) 0
= - im r ξ + 2Λ µ b △v + 2im r 2 u - v r 2 + dµ b dr Dv + im r u - v r , (75)
uDW b = -i α ξ + 2Λ µ b △ w + dµ dr (Dw + iαu) +2Λ 1 r D [r (µ t -µ b ) (Dw + iαu)] . (76)
with

Du + u r + i m r v + α w = 0. (77)
Equations ( 74)-( 76) define an eigenvalue problem, where Λ is the eigenvalue to compute and

Re 1 (α, m, n, Hb) = max Λ. ( 78
)
The eigenvalue problem (74)-( 76) is solved numerically using the procedure described in §5 with mapping [r 0 , 1] into [0, 1] and the transformations ( 38) - [START_REF] Mack | The initial stage of transition in cylindrical pipe flow: role of optimal base-flow distorsion[END_REF]. The results of the numerical computation are presented in terms of marginal stability curves, i.e. curves in the α, Re 1 plane, that separate regions of initial energy density growth from regions of initial energy decay. Influence of the plug radius In order to highlight the influence of the unyielded zone radius, r 0 , on the global stability conditions, the Hb terms are canceled artificially and n is fixed to 1. Figure 19 shows in this particular case, i.e. Couette-Poiseuille flow of a Newtonian fluid, that the most dangerous perturbation is a streamwise vortices (α = 0) with an azimuthal wavenumber m increasing with r 0 . The maximum Reynolds number, Re CE = min α,m

Re 1 ( α, m), that ensures no energy growth, depends weakly on r 0 . When r 0 → 1, Re cE tends towards the value obtained by Nouar et al. [START_REF] Nouar | Modal and non-modal linear stability of the plane Bingham-Poiseuille flow[END_REF] for a Bingham plane Poiseuille flow, i.e. Re cE = 79.82. When r 0 → 0, Re cE tends to Re cE = 91.08. This value is different from that obtained by Joseph and Carmi [START_REF] Joseph | Stability of Poiseuille flow in pipes, annuli and channels[END_REF] for Hagen-Poiseuille flow of a Newtonian fluid, i.e. Re cE = 81.49. This singularity is a consequence of the boundary conditions at the yield surface. For sufficiently high value of Hb, the critical mode is axisymmetric, as shown in figure [START_REF] Balmforth | Yielding to stress: Recent developments in viscoplsatic fluid mechanics[END_REF], where variations of the critical axial α and azimuthal m wavenumbers as function of Hb are represented. It is worth noting that for the ax- isymmetric perturbation, the viscous dissipation involves the tangent viscosity µ t < µ b . The variation of the critical Reynolds number, Re CE is given by the curve (1) in figure 23. The conditional stability derived by Nouar and Frigaard [START_REF] Nouar | Nonlinear stability of Poiseuille flow of a Bingham fluid: theoretical results and comparison with phenomenological criteria[END_REF] for Bingham fluid is also represented by curve [START_REF] Orr | The stability or instability of the steady motions of a perfect liquid and of a viscous liquid. part I: A perfect liquid[END_REF]. In this analysis, the effect of the yield stress is limited to the modification of the unyielded zone radius. For large Hb ( Hb ≥ 30), Re cE varies as Hb

Conclusion

Linear stability analysis and receptivity of Hagen-Poiseuille flow of Herschel-Bulkley fluid are investigated using modal and non modal approaches. The base flow is mainly characterized by (i) a central plug zone of radius r 0 moving with a constant velocity and (ii) a nonlinear variation of the viscosity between the wall and the yield surface. Two dimensionless parameters govern the problem: the shear-thinning index n and the modified Herschel-Bulkley number Hb or equivalently n and r 0 . In the linear stability analysis, taking into account the viscosity perturbation leads to an anisotropy of the stress-tensor perturbation. The component τ ′ rz = τ ′ zr involves the tangent viscosity µ t , which is independent of the yield stress, while the other components τ ′ ij , ij = rz, zr involve the effective viscosity µ b > µ t . The results of the modal approach show that: (i) In the case of an axisymmetric perturbation, the influence of Hb (dimensionless yield stress) is reduced mainly to the modification of the width of the yielded zone; (ii) In the case of an oblique perturbation, the Hb terms contribute to stabilize the flow through an increase of the viscous dissipation; (iii) The Hagen-Poiseuille flow of Herschel-Bulkley fluid is asymptotically stable to infinitesimal disturbances.

Because of the non-normality of the linear stability operators (with respect to a scalar product based on the energy norm), a transient algebraic amplification of the kinetic energy of the perturbation is observed. The influence of the rheological parameters on the transient growth and the characteristics of the optimal perturbation is studied. It is shown that: (i) with increasing the dimensionless yield stress Hb, the amplification of the kinetic energy G max is reduced. However, for an axisymmetric perturbation, Hb has practically no effect on G max . (ii) For Hb << 1, the optimal perturbation is in the form of longitudinal rolls. The amplification of the kinetic energy is governed by the lift-up mechanism. The scaling laws in Re 2 for G opt and Re for t opt are recovered. (iii) With increasing Hb, shear-thinning effects become more significant, the optimal perturbation is oblique in the plane (r, z) with an axial wavenumber α = O(1). The optimal azimuthal wavenumber increases with decreasing the width of the yielded zone. The optimal perturbation evolves by the effects of Orr and Lift-up mechanisms. (iv) For sufficiently large value of Hb, the optimal perturbation is axisymmetric and the transient growth occurs only by the Orr-mechanism.

Concerning the energy stability analysis and the determination of the maximum Reynolds number Re CE below which the kinetic energy of the perturbation decreases uniformly with time, the numerical results show that for Hb << 1, with 0.3 ≤ n ≤ 1, the most dangerous perturbation is in the form of longitudinal rolls with Re cE ≈ 91.1 -26(1 -n). For large Hb, the most dangerous perturbation is axisymmetric and Re cE ≈ 218n 0.46 Hb 0.5 , with α ∝ Hb -0.5 . If the viscosity perturbation is not taken into account, Re cE ∝ Hb.

The different scaling laws, given in sections §6- §7, involve the dimensionless parameters Hb and Re defined with a generalized viscosity, μgen ( §2) derived from the dimensionless form of the momentum equations. However, the structure of: (i) the base flow, the less stable mode and (iii) the optimal perturbation, suggests that wall shear viscosity is a more appropriate choice of viscosity scale. For an axisymmetric perturbation, tangent wall shear viscosity is even more relevant. Results obtained in terms of Re and can be expressed in terms of wall Reynolds number, Re w , using the relation Re w = Re (1 -r 0 ) n [n/(n + 1)] n-1 . Using Re w rather than Re does not globally change the conclusions given above. Nevertheless, the comparison of yield-stress shear-thinning fluids among themselves does require to indicate the scale viscosity adopted.

This work puts just a brick to a more comprehensive building to identify possible paths of transition. Experimental studies of transition to turbulence in a pipe for yield stress fluids were performed by Peixinho et al. [START_REF] Peixinho | Transitional and turbulent flow of yield stress fluid in a pipe[END_REF], Esmael et al. [START_REF] Esmael | Transitional flow of a nonnewtonian fluid in a pipe: Experimental evidence of weak turbulence induced by shear-thinning behavior[END_REF] and Guzel et al. [START_REF] Guzel | Observation of laminarturbulent transition of a yield stress fluid in Hagen-Poiseuille flow[END_REF]. A weakly turbulent flow is observed, where the time-averaged axial velocity profiles exhibit an asymmetry. The mechanism associated to this nonlinear asymmetry is not yet clarified. The present work can be pursued in two directions. The first one, also based on linear stability equations, starts from the realization that the base flow around which the linearization is performed is just an idealization, whereas in reality small defects inevitably occur. This approach proposed by Bottaro et al. [START_REF] Bottaro | The effect of base flow variation on flow stability[END_REF] has shown some success in capturing features of transition in pipe flow (Gavarinni et al. [START_REF] Mack | The initial stage of transition in cylindrical pipe flow: role of optimal base-flow distorsion[END_REF], Ben-Dov & Cohen [START_REF] Bendov | Critical Reynolds number for a natural transition to turbulence in pipe flows[END_REF], [START_REF] Bendov | Instability of optimal non-axisymmetric base flow deviations in pipe poiseuille flow[END_REF]). It was applied for the first time to the channel flow of Bingham flow by Nouar et al. [START_REF] Nouar | Stability of the flow of a Bingham fluid in a channel: eigenvalue sensitivity, minimal defects and scaling laws of transition[END_REF], where scaling laws for transition were proposed. The second direction, is to determine a weakly nonlinear optimal perturbation. This approach, proposed recently by Pralits et al. [START_REF] Pralits | Weakly nonlinear optimal perturbations[END_REF] allows to identify initial states which relaminarize from those which grow without bound.

A Long wave approximation

Following the asymptotic method introduced by Cowley and Smith [START_REF] Cowley | On the stability of Poiseuille-Couette flow: a bifurcation from infinity[END_REF] for a Couette-Poiseuille flow, the following "long-wave" eigenvalue problem is derived for an axisymmetric perturbation:

3

 3 Base flow The base flow is a one directional shear flow, U = U b = W b (r)e z driven by a constant pressure gradient. The subscript b means base flow. The only non zero elements of the strain rate tensor are off-diagonal γrz,b = γzr,b , so that the deviatoric stress tensor elements are all zero except for τ rz,b = τ zr,b . The momentum equations reduce to:

Fig. 1 .

 1 Fig. 1. shematic representation of the base flow.

  base flow depends on two dimensionless parameters, n and Hb or n and r 0 . The objective of the present work is to examine the influence of n and Hb on the stability of the Hagen-Poiseuille flow of a Herschel-Bulkley r 0 )/(1 -r 0 )

  r 0 )/(1 -r 0 )

Fig. 2 .

 2 Fig. 2. Viscosity profiles between the yield surface and the wall. (a) Influence of the radius of the plug zone on the viscosity stratification for n = 0.5: (1) r 0 = 0.01, Hb = 0.0176; (2) r 0 = 0.1, Hb = 0.203; (3) r 0 = 0.5, Hb = 2.45; (4) r 0 = 0.8, Hb = 15.49. (b) Influence of the shear-thinning index with r 0 = 0.1: (1) n = 1; (2) n = 0.7; (3) n = 0.5; (4) n = 0.3.

  η dr. It is clear that the kinetic energy of the perturbation decreases uniformly

Fig. 3 .

 3 Fig. 3. (a) Eigenvalue spectrum for a pipe flow of Herschel-Bulkley fluid at Re = 5 × 10 3 , n = 0.5, r 0 = 0.5, α = 0 and m = 0. (b) Maximum growth rate Ci,max vs Re for two different values of shear-thinning index: (1) n = 1 and (2) n = 0.3

Fig. 4 .

 4 Fig. 4. Modified Orr-Sommerfeld spectra for Herschel-Bulkley fluid flow in a pipe at Re = 5× 10 3 , α = 1 and m = 0. (a) n = 1 and the Hb terms are canceled artificially to highlight the influence of r 0 : (•) r 0 = 0.01 and ( ) r 0 = 0.8. (b) Influence of Hb terms: (•) r 0 = 0.8, Hb = 0.96; ( ) r 0 = 0.8, Hb terms neglected

Fig. 5 .

 5 Fig. 5. Herschel-Bulkley fluid flow at Re = 5000, with r 0 = 0.5, α = 1 and m = 0. (a) Growth rate of the least stable wall mode as function of Re: (1) n = 1; (2) n = 0.7; (3) n = 5; (4) n = 0.3. (b) Phase speed of the least stable mode as function of Re: (1) n = 1; (2) n = 0.7; (3) n = 5; (4) n = 0.3.

Fig. 5 .

 5 Fig. 5. We have represented |C i,max | and Cr,max of the least stable mode, vs Re for different values of n.As n decreases, the viscosity around the plug zone becomes larger and the wall mode is squeezed against the wall. The phase velocity decreases and the growth rate comes closer to the unstable half plane. For large values of Re, the scaling with Re is the same as for a Newtonian fluid. Similarly, for interface and mean modes, the scaling with Re is the same as for a Newtonian fluid.

Fig. 6 .

 6 Fig. 6. (a) Eigenvalue spectrum at Re = 5 × 10 3 , α = 0, m = 1, n = 0.5 and r 0 = 0.5. (b) Evolution of the maximum growth rate as a function of Re for n = 1, curve (1) and n = 0.3, curve (2).

Fig. 7 .Fig. 8 .

 78 Fig. 7. Eigenvalues distribution at α = 1, m = 1, Re = 5×10 3 , r 0 = 0.1 and n = 0.5.

Fig. 9 .

 9 Fig. 9. Spectral portrait of the modified Orr-Sommerfeld operator at Re = 5000, α = 1, for a Herschel-Bulkley fluid: influence of the radius plug zone. Continuous lines are the boundaries of the ǫ-pseudospectra and the dashed line is the numerical range. (a) n = 0.5, r 0 = 0.01 ( Hb = 0.017), (b) n = 0.5, r 0 = 0.8 ( Hb = 6.93) .

Fig. 10 .

 10 Fig. 10. Spectral portrait of the linear operator for streamwise homogeneous perturbation at Re = 5000, n = 0.5, α = 0 and m = 1. The dashed curve is the numerical range of the linear operator considered. Continuous line curves delimit the isovalues of the ǫ-pseudospectra. (a) r 0 = 0.01: (1) ǫ = 2.5 × 10 -3 , (2) ǫ = 1.5 × 10 -3 , (3) ǫ = 7 × 10 -4 , (4) ǫ = 3 × 10 -4 ; (b) r 0 = 0.8: (1) ǫ = 0.25, (2) ǫ = 0.13, (3) ǫ = 9 × 10 -2 , (4) ǫ = 5 × 10 -2 , (5) ǫ = 3 × 10 -2

Fig. 11 .

 11 Fig. 11. Norm of the minimal destabilizing perturbation as a function of the modified Reynolds number for two values of the shear-thinning index, n = 0.1 and 0.3 and two values of r 0 : (a) r 0 = 0.01; (b) r 0 = 0.8. The axial and azimuthal wavenumbers are α = 1 and m = 1 respectively.)

Fig. 12 .Fig. 13 .

 1213 Fig. 12. kinetic energy amplification at Re = 3000, n = 1 and different values of the plug radius r 0 : r 0 = 0.01 ( Hb = 0.02), r 0 = 0.1 ( Hb = 0.22) and r 0 = 0.2 ( Hb = 0.5).(a) α = 0, m = 1. The Newtonian curve is given by Schmid and Henningson[START_REF] Schmid | Optimal energy density growth in Hagen-Poiseuille flow[END_REF]. (b) α = 1, m = 1. The Newtonian curve is given by Meseguer and Trefethen[START_REF] Meseguer | A spectral Petrov-Galerkin formulation for pipe flow I: Linear stability and transient growth[END_REF].

Fig. 14 .

 14 Fig. 14. Optimal perturbation and optimal streaks at Re = 5000, n = 0.5, r 0 = 0.01 ( Hb = 0.0175): α = 0 and m = 2. (a) Velocity vectors ue r + ve θ of the optimal perturbation at t = 0. (b) Axial velocity w contours at t = t opt = 165

  Fig. 15. Temporal evolution of the optimal perturbation. n = 0.5, r 0 = 0.5 ( Hb = 1.732), Re = 5000, with αopt = 0.28, m opt = 7, topt = 38.6. Axial velocity contours at (a) t = 0; (b) t = topt /2; (c) t = topt ; (d) t = (3/2) topt .

Fig. 16 .

 16 Fig. 16. Contribution of the different Reynolds-stress components to the energy amplification for an optimal oblique perturbation at Re = 5000, with (a) n = 0.5, r 0 = 0.5, α = 0.28, m = 7; (b) n = 0.5, r 0 = 0.8, α = 1.54, m = 12; (c) n = 0.5, r 0 = 0.9, α = 1.715, m = 17; (d) Case where µ t -µ b terms and their derivatives are canceled, n = 0.5, r 0 = 0.9, α = 0.6, m = 44. (1) -γrz u d w * d ;(2) -γrz u q w * q ; (3) -γrz u d w * q ; (4) -γrz u q w *

  Fig. 17. Temporal evolution of the optimal perturbation. n = 0.5, r 0 = 0.96 ( Hb = 41.47), Re = 5000, with αopt = 0.94, m opt = 0, topt = 16.2. Axial velocity contours at (a) t = 0; (b) t = topt /2; (c) t = topt ; (d) t = (3/2) topt .

Fig. 18 .

 18 Fig. 18. Variation of G opt / Re 2 and T opt / Re as a function of the radius of the plug zone for n = 0.5 and different values of the Re: (1) Re = 3000, (2) Re = 4000, (3) Re = 5000, (4) Re = 7000.

Fig. 19 .

 19 Fig. 19. Newtonian Couette-Poiseuille flow. Boundaries in the (α, Re) delimiting zones of global stability. (a) r 0 = 0.001, Re cE = 91.08 ; (b) r 0 = 0.2, Re cE = 85.06; (c) r 0 = 0.5, Re cE = 81.84; (d) r 0 = 0.9, Re cE = 80.06. The curves are labeled with the azimuthal wavenumberViscosity stratification without anisotropy termsThe influence of the viscosity stratification, on the global stability is illustrated in Fig.20. The terms arising from the viscosity perturbations are canceled. It is worthy to note that: (i) the most dangerous perturbation is a streamwise perturbation with an azimuthal wavenumber m c larger that obtained for a Newtonian Couette-Poiseuille flow; (ii) the critical Reynolds number Re CE increases strongly with increasing r 0 (or Hb), because of the increase of the viscous dissipation. From Hb ≈ 30, the numerical results show that Re cE ∝ Hb.

Fig. 20 .Fig. 21 .

 2021 Fig. 20. Influence of viscosity stratification without anisotropy terms. Boundaries in the (α, Re) delimiting zones of global stability. (a) r 0 = 0.5, n = 0.5, Hb = 1.73, m c = 7, Re cE = 237.9 ; (b) r 0 = 0.9, Hb = 15.59, m c = 42, Re cE = 1342. The curves are labeled by the azimuthal wavenumber.

2 Fig. 22 .

 222 Fig. 22. Evolution of the critical axial and azimuthal wavenumbers (α c and m c as function of Hb. The shear-thinning index is fixed: n = 0.5.

Fig. 23 .

 23 Fig. 23. Critical Reynolds numbers, Re CE versus Hb. The shear-thinning index is fixed n = 0.5. CPF denotes the limit of a Newtonian Couette-Poiseuille flow between coaxial cylinder when r 0 → 0.

  iλ (W b -ω) φ ′′φ = ru, λ = α Re, ω = C/ α and F = |DW b |. Equation A.[START_REF] Landhall | A note on an algebraic instability of inviscid parallel shear flows[END_REF] is solved by the same method as that described in §5.

Figure A. 1

 1 shows the evolution of the eigenvalue ω as a function of λ. We have not found any instability.

Fig. A. 1 .

 1 Fig. A.1. Maximum temporal amplification rate of the perturbation as function of λ = α Re at n = 0.5 and two values of the plug zone radius: ( ) r 0 = 0.8 and (•) r 0 = 0.1 .

Table 1

 1 Fitting coefficients in the scaling laws G opt / Re 2 and t opt / Re for different values of the shear-thinning index n.

		3	0.5	0.7	1
	A n × 10 6	66.1	68.1	68.1	62.9
	B n	-2.22 -3.28 -3.55 -3.62
	A ′ n × 10 3	24.6	32.1	35.4	32.9
	B ′ n	-2.72 -2.81 -2.96 -2.8

where u ′ r is the real part of u(r, t) exp [i(α z + mθ)] and similarly for v ′ r and w ′ r . Hence, E can be written as

The temporal variation of the mean kinetic energy dE dt is obtained by multiplying the linearized perturbation equations ( 29), ( 30) and ( 31) by u * , v * and w * respectively and then integrated over the domain Ω, with the condition u = 0 at the wall and at the interface. The result is a modified Reynolds-Orr equation, which after some algebra, is given by

where I (u), V(u) and A(u) denote the inertial, viscous and anisotropy terms, defined by

In the above equations

There is no energy growth of the perturbation, i.e, the basic flow remains "energy stable", if dE dt < 0. Define Re 1 as the largest value of Re for which this condition is satisfied

where u is an admissible perturbation satisfying the continuity equation and the boundary conditions. This constrained optimization problem is solved using variational calculus. The corresponding Euler-Lagrange equations are