The flows in the pipes are presented in several industrial applications such as heat exchangers, chemical reactors, gas turbine cooling and heating systems and combustion chambers and mixing systems.
The properties of a flowing fluid vary with the type of fluid. By adding, for example polymers or micro-metric particles that can interact with each other to a liquid under flow; its properties are greatly modified. It can react to hydrodynamic friction forces by changing equilibrium stator conformation through reorientation and deformation. A colloidal suspension can react to hydrodynamic forces by modifying the spatial distribution of the particles. It is these deformities and this reorganization which are at the origin of the non-Newtonian properties of the polymer solutions. When the polymer chains are quite long and flexible, they are added to a fluid, which acquires viscoelastic properties. Viscoelastic fluids have elastic properties characterized by a relaxation time. This elasticity plays a fundamental role and would be of the same order of magnitude as the characteristic time of the flow.
When a viscoelastic fluid is in flow, the polymer chains are stretched, which generates elastic stresses. These are the function of the history of the movement and the deformations undergone by the fluid particles along their trajectory. It follows a non-linear relationship between the elastic stresses thus generated and the rate of deformation. The conservation equations of the momentum of such a fluid differ from those of Navier-Stokes newtonian fluids by an additional term which involves elastic stresses due to the stretching of polymer chains by hydrodynamic friction forces.
The motivation for this study is largely related to the fact that viscoelastic fluids are the basis of many industrial applications, particularly in the polymer industry, the paper industry, the food industry, etc. and the stability of the flow can have consequences on the final product (quality of the product, quantity of production, etc.).
Three classical incompressible shear flows have received attention in connection with phenomenon of instability. Plane poiseuille and plane-Couette channel flows are the most accessible to theoretical analysis and numerical computation. Pipe poiseuille or Hagen-poiseuille flows is so accessible to laboratory experiments. It was the context of flow in a pipe that Reynolds in 1883 identified the basic problem of this field: when and how do high-speed flows undergo transition from the laminar state to more complicated states as puffs, slugs, etc. The laminar flow of a fluid through a finite dimensional cylindrical conduit is mathematically stable, but in practice such a flow may become unstable if the Reynolds is high enough. It is generally accepted that one of the explanations for this phenomenon is that, although the laminar flow is stable for infinitesimal perturbations of the velocity, certain amplitudes of the perturbations are sufficient to generate a large number of Reynolds.
The present paper concerns the numerical simulation of pipe flow. We aim to focus on some of the principal mechanisms by which small perturbations in high speed flow may grow as they flow downstream. Our focus is on the elucidation of key mechanisms involved in the phenomenon of instabilities of the flow of a viscoelastic fluid. On the other hand experience shows that the flow of a viscoelastic fluid can become even with a low Reynolds number.
With these aims in mind, and with an acute awareness of the continuity advance of desktop and laptop computers, we have written a Petrov-Galerkin spectral code for pipe flow in FORTRAN. In principle we are thus engaged in direct numerical simulation (DNS) of the Navier-Stokes equations.
Meseguer & Trefethen  conducted this study with a newtonian fluid. For this kind of fluid the flow is stable up to . Esmael  and then Carranza 2012  have in turn studied this problem with another type of fluid that described by the rheofluidifying fluid model and have comed to the same conclusion. For a viscoelastic fluid described by the Oldroyd-B model  , we try to explain the origin of the instabilities by studying the terms relating respectively to viscous effects Re and elastic effects .
Previous numerical simulation of (nonlinear) pipe flow has included works of Boberg and Brosa  , Komminaho and Johasson  Zhang et al.  and Zikanov  . In addition, various previous authors have simulated linearized pipe flow     . Our methods could be summarized as a solenoidal scheme for the pipe of the sort proposed by Leonard and Wray  .
At our level the problem we are studying has an additive term ( ). We try to show the influence of this term on the structure of this flow. To our knowledge, a theoretical study has not been done before for this type of fluid.
2. Equations Governing the Problem
We study the flow of a viscoelastic fluid along a cylinder of circular section of horizontal axis (Oz). Such a flow can be described by a cylindrical coordinates system .
The flow of a viscoelastic fluid can be described using three main equations: the momentum conservation equation, the mass conservation Equation (or continuity equation), and the behavioral equation of the extra-constraint (Chupin 2003  and Oldroyd  ). We consider the flow an incompressible viscoelastic fluid with dynamic viscosity and density driven by an extremal constant axial pressure gradient.
The equations of conservation of momentum, continuity, and constitutive form in dimensionless form are used by using the maximum speed in the established threshold regime, the radius R of the pipe and the quantity as a speed scale, length and Pessure/stress respectively:
By introducing these dimensionless parameters into the equations governing the problem, we obtain the following dimensionless equations:
We will omit the sign ~ voluntarily to lighten the scriptures.
We focus here on studying the flow of a viscoelastic fluid described by the Oldroyd-B model and subjected to disturbances.
This disturbed flow is considered as the superposition of a disturbance and of basic flow . Disturbance equations are obtained by subtracting the written conservation equations for the disturbed flow , those satisfied by the basic flow. We obtain after calculation:
Is the tensor of deformation rates and
is the vorticity tensor. is defined as a parameter which value is between −1 and +1.
Considering in addition that the elastic contribution of the extra stress is a perturbation of the Newtonian one which is equal to , it comes:
By , we define the none-dimensional velocity field of the fluid whose radial, azimuth and axial components are given in the order by the triple , is an additional term called extra stress and is the parameter of delay defined by
The system is closed with the following boundary and initial conditions:
where is defined as the Reynolds number with the kinematic viscosity of the fluid, is the Weissenberg’s number.
3. Procedure and Numerical Bases
We consider that the variation of the perturbation of the fields of velocity, pressure and extra stress is periodic along the azimuthal and axial directions. These considerations make it possible to approximate in the following form.
This approximation was used by Meseguer & Trefethen (2001)  .
The continuity equation leads to a linear dependence between the three components of leading to a system with two degrees of freedom. Therefore we note:
For long periods, this approximation can be written in the form:
For the substitution of Equation (12) in the Equation (6) the problem result in a system of ordinary differential equations of coefficients . This substitution is followed by a projection based on the scalar product:
where designates the conjugated function of F.
The basic functions will be chosen so that integrand is even, which will result in:
The basic functions are also chosen so that the projection of the pressure term is zero. For this purpose functions based on Chebyshev polynomial were considered. This choice imposes two essential criteria. The first consists of choosing weights associated with Chebyshev polynomials, which makes it easier to calculate the scalar product. The second criterion relates to the approximation of the integral by a Gauss-Chebyshev-Lobatto quadrature of the form:
is the Chebyshev polynomial defined by
The basic functions and tests were proposed by Leonard & Wray 1982  then used by Meseguer & Trefethen 2001  and are defined as follows:
Using the Fourier representation, the continuity equation becomes
4. Numerical Implementation
Replacing by the basic funct:ions and projecting on the basis of the test functions, it comes:
The numerical method chosen is adapted from work in the early 1980 by Leonard & Wray  .
With this Petrov-Galerkin procedure that was then used by Meseguer & Trefethen 2000  and Meseguer & Mellibovesky 2007  , we obtain:
The matrices and derive from the projection of the nonlinear terms and can be calculated with a pseudo-spectral method by Fast Fourier Transform (FFT).
The problem obtained is a problem with the generalized eigenvalues that we will solve numerically with QZ algorithm used by J. P. Berlioz  . To implement this method, we will use Housholder’s unitary reflection matrices and Givens rotation matrices (A. Quarteroni, R. Sacco and F. Saleri  ) and (L. Amodei and J.P. Dedieu  ).
5. Results and Discussions
We will treat this problem according to four cases that are:
5.1. Case of One-Dimensional Disturbance (n = 0; l = 0)
Let us first note that for such a flow, the radial component of the velocity is zero . On the other hand, the Figure 1 clearly shows that the fluid is Newtonian or viscoelastic , the eigenvalues are purely real and negative that it’s worth the Reynolds’ number Re. For both types of fluid, the flow remains stable whatever the Reynolds number.
5.2. Case of Homogeneous Disturbance (n ≠ 0; k0 = 0)
The eigenvalues related to this flow are real and negative for the case of the Newtonian fluid. The result does not change when increasing the number of Reynolds up to values close to 107 (Figure 2(b)). On the other hand, this is not always the case for a viscoelastic fluid. For this type of fluid the real part of some eigenvalues is positive for even for a Reynolds number equal to . However, instability occurs with a higher Reynolds number when the Weissenberg’s number We decreases (Figure 2(c)). The flow of a Newtonian fluid is stable whereas that of a viscoelastic fluid is unstable for a Reynolds number ( ).
5.3. Case of an Axisymmetric Disturbance (n = 0; k0 ≠ 0)
The real parts of the eigenvalues remain negative for the Newtonian fluids, whereas for the same flow conditions, the real parts of the eigenvalues are not all negative for a viscoelastic fluid, which testifies to the appearance of instabilities in flow. It is also noted that beyond a certain value of the Reynolds number, instabilities appear in the flow itself for the Newtonian fluid, which is not the case
Figure 1. Spectrum of the eigenvalues of a one-dimensional disturbance , case of Newtonian fluid (in red), case of viscoelastic fluid (black), with Reynolds’ number .
(a) (b) (c)
Figure 2. (a) Spectrum of the eigenvalues of a homogeneous disturbance along the axis with , case of a Newtonian fluid (black), case of a viscoelastic fluid (red) with Reynolds number and Weissenberg’s number ; (b) Spectrum of the eigenvalues of a homogeneous disturbance along the axis with for different values of the Reynolds number: case of a Newtonian fluid ; (c) Spectrum of the eigenvalues of a homogeneous disturbance along the axis with for different values of the Reynolds number: case of a viscoelastic fluid with .
for one-dimensional and homogeneous flows. Indeed, the term convection , responsible for the exchanges between the base flow and the disturbed flow, is a source of instabilities, a term which is also nil for the Newtonian one-dimensional and homogeneous flows. In addition, it is noted (Figure 3(b)) that the greater the elasticity of the fluid, the greater the instability is important. From this it can be deduced that elasticity is also a source of instability.
5.4. Case of Three-Dimensional Perturbation (n ≠ 0; l ≠ 0)
The graph clearly shows a spectrum of eigenvalues whose real parts are all negative for the Newtonian fluid. For the viscoelastic fluid, however, the real parts of
Figure 3. (a) Spectrum of the eigenvalues of an axisymmetric disturbance with , case of the Newtonian fluid (black), case of the viscoelastic fluid (red) with and ; (b) Evolution of the most unstable eigenvalue according to the delay parameter with Reynolds’ number and Weissenberg’s number ; (b) shows that the one-dimensional flow of a viscoelastic fluid with a Reynolds number is stable for a delay parameter ranging from 0 (Newtonian fluid) to 0.35. The homogeneous flow of a viscoelastic fluid remains stable as long as the retardation parameter is less than 0.2 and then becomes unstable beyond this value, but this flow remains stable until for . However at beyond , this threshold of instability is reached from .
the eigenvalues are not negative and would indicate that the fluid flows unstable for the same Reynolds number (Figure 4).
Figure 4. Spectrum of the eigenvalues of a three-dimensional perturbation case of the Newtonian fluid (black), case of viscoelastic fluid (red) with Reynolds number and Weissenberg’s number .
The analysis of the different results shows that the linear flow of a Newtonian fluid is stable, whereas that of a viscoelastic fluid with a number of Weissenberg and for a delay parameter ( ) is unstable at from except for the one-dimensional case. Indeed, for the one-dimensional and homogenous cases, the term relating to the convection which is responsible for the exchange of energies between the base flow and the disturbed flow is almost nil. It can therefore be said that viscoelasticity and the term convection are at the origin of the instabilities observed.
The algorithm QZ makes possible to determine the eigenvalues and optionally the eigenvectors of the problem with generalized eigenvalues . The idea is to transform the matrices A and B into similar upper triangular matrices. If A and B are respectively upper hessenberg and triangular matrices respectively. We call the
QZ decomposition of pair , the double factorization with Q unitary matrix, Z upper Hessenberg, R and T upper triangular. We will note later on we can construct a unitary decomposing matrix A, as a product of
matrices of elementary rotation: where each is a complex rotation matrix. modifies the first and second lines of B possibly creating a non-zero element in position which will have to be cancelled by a complex rotation matrix modifying only the first and the second column of . Note and .
is triangular superior.
Let be the upper triangular, let be the rotation matrix modifying only the kth and (k + 1)th row of possibly creating a non-zero element in position which will have to be canceled by a matrix complex rotation only modifying the kth and (k + 1)th columns of .
Let’s say is upper triangular, the process stops for and we can write , .
If B is invertible, it is obvious that Q realizes a unitary and triangular decomposition of and a unitary decomposition of , in fact:
and QA are upper triangular. We obtain the eigenvalues by calculating the ratios of diagonal elements of QA and .
θ: Azimuthal coordinate
λ: Relaxation time of the fluid
: delay time
ν: Kinematic viscosity of the fluid
ρ: Density of the fluid
σ: Tensor of extra stress
: Disturbance of extra-stress
: Newtonian contribution of the disturbance of the extra-stress
: Elastic contribution of the disturbance of the extra-stress
η: Dynamic viscosity of the fluid
ϕ: Basis function
ψ: Test function
ω: delay parameter
Ci: imaginary part of the eigenvalue
Cr: real part of the eigenvalue
f: density force
l: axial mode
k0: Axial wave number
n: Azimuthal mode
r: radial coordinate
R: radius of the cylinder
Re: Reynolds number
W0: Maximum speed of the base flow, it has for direction that of the axis of the cylinder
We: Weissenberg’s number
Wb: Axial velocity of the base flow
: Newtonian contribution of basic flow velocity
: Elastic contribution of basic flow velocity
z: axial coordinate
 Komminaho, J. and Johansson, A.V. (2000) Development of a Spectrally Accurate DNS Code for Cylindrical Geometries, Part V of J. Komminaho, Direct Numerical Simulation of Turbulent Flow in Plane and Cylindrical Geometries. Docoral Thesis, Royal Institute of Technology, Stockholm.
 Zhang, Y., Gandhi, A., Tomboulides, A.G. and Orzag, S.A. (1994) Simulation of Pipe Flow. In: Application of Direct and Large Eddy Simulation to Transition and Turbulence, AGARD Conference, Crete, 17.1-17.9.
 Leonard, A. and Wray, A. (1982) A New Numerical Method for the Simulation of Three-Dimensional Flow in a Pipe. In: Krause, E., Ed., Proceedings of the 8th International Conference on Numerical Methods in Fluid Dynamics on Numerical Methods in Fluid Dynamics, Springer, Berlin, 335-342.
 Trefethen, A.E., Trefethen, L.N. and Schmid, P.J. (1999) Spectra and Pseudospectra for pipe Poiseuille Flow. Computer Methods in Applied Mechanics and Engineering, 175, 413-420.
 Meseguer, A. and Trefethen, L.N. (2000) A Spectral Petrov-Galerkin Formulation for Pipe Flow (I): Linear Stability and Transient Growth. Numerical Analysis Group, Computing Lab, Oxford University, Rep. 00/18.
 Meseguer, A. and Mellibovsky, F. (2007) On a solenoidal Fourier-Chebyshev Spectral Method for Stability Analysis of the Hagen-Poiseuille Flow. Applied Numerical Mathematics, 57, 920-938.