Effects of the Reynolds and Weissenberg’s Numbers on the Stability of Linear Pipe Flow of Viscous Fluid

Show more

1. Introduction

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 [1] conducted this study with a newtonian fluid. For this kind of fluid the flow is stable up to ${R}_{e}={10}^{7}$ . Esmael [2] and then Carranza 2012 [3] 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 [4] , we try to explain the origin of the instabilities by studying the terms relating respectively to viscous effects Re and elastic effects $\omega $ .

Previous numerical simulation of (nonlinear) pipe flow has included works of Boberg and Brosa [5] , Komminaho and Johasson [6] Zhang et al. [7] and Zikanov [8] . In addition, various previous authors have simulated linearized pipe flow [9] [10] [11] [12] . Our methods could be summarized as a solenoidal scheme for the pipe of the sort proposed by Leonard and Wray [10] .

At our level the problem we are studying has an additive term ( $div\sigma $ ). 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 $\left(r,\theta ,z\right)$ .

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 [13] and Oldroyd [4] ). We consider the flow an incompressible viscoelastic fluid with dynamic viscosity $\eta $ and density $\rho $ driven by an extremal constant axial pressure gradient.

$\{\begin{array}{l}\rho \left(\frac{\partial \stackrel{\to}{u}}{\partial t}+\left(\stackrel{\to}{u}\cdot \nabla \right)\stackrel{\to}{u}\right)=-\nabla P+\stackrel{\to}{f}+div\left\{2\eta \left(1-\omega \right)D\left(\stackrel{\to}{u}\right)\right\}+div\sigma \\ \frac{\partial \rho}{\partial t}+div\left(\rho \stackrel{\to}{u}\right)=0\\ \frac{D\sigma}{Dt}+\frac{\sigma}{\lambda}=\frac{2\eta \omega}{\lambda}D\left(\stackrel{\to}{u}\right)\end{array}$ (1)

The equations of conservation of momentum, continuity, and constitutive form in dimensionless form are used by using the maximum speed ${W}_{0}$ in the established threshold regime, the radius R of the pipe and the quantity $\rho {W}_{0}^{2}$ as a speed scale, length and Pessure/stress respectively:

$\stackrel{\u02dc}{p}=\frac{p}{\rho {W}_{0}^{2}};\text{\hspace{0.17em}}\text{\hspace{0.17em}}\stackrel{\u02dc}{\sigma}=\frac{\sigma}{\rho {W}_{0}^{2}};\text{\hspace{0.17em}}\text{\hspace{0.17em}}\stackrel{\to}{\stackrel{\u02dc}{u}}=\frac{\stackrel{\to}{u}}{{W}_{0}};\text{\hspace{0.17em}}\text{\hspace{0.17em}}\stackrel{\u02dc}{r}=\frac{r}{R};\text{\hspace{0.17em}}\text{\hspace{0.17em}}\stackrel{\u02dc}{z}=\frac{z}{R};\text{\hspace{0.17em}}\text{\hspace{0.17em}}\stackrel{\u02dc}{t}=\frac{t{W}_{0}}{R},\text{\hspace{0.17em}}\text{\hspace{0.17em}}\stackrel{\u02dc}{\nabla}=R\nabla $

By introducing these dimensionless parameters into the equations governing the problem, we obtain the following dimensionless equations:

$\{\begin{array}{l}\left\{\frac{\partial \stackrel{\to}{\stackrel{\u02dc}{u}}}{\partial t}+\left(\stackrel{\to}{\stackrel{\u02dc}{u}}\cdot \stackrel{\to}{\stackrel{\u02dc}{\nabla}}\right)\stackrel{\to}{\stackrel{\u02dc}{u}}\right\}=-\stackrel{\to}{\stackrel{\u02dc}{\nabla}}\stackrel{\u02dc}{p}+\frac{\left(1-\omega \right)}{{R}_{e}}\stackrel{\u02dc}{\Delta}\stackrel{\to}{\stackrel{\u02dc}{u}}+\stackrel{\to}{\stackrel{\u02dc}{\nabla}}\cdot \stackrel{\u02dc}{\sigma}\\ \stackrel{\to}{\stackrel{\u02dc}{\nabla}}\cdot \stackrel{\to}{\stackrel{\u02dc}{u}}=0\\ {R}_{e}\frac{D\stackrel{\u02dc}{\sigma}}{D\stackrel{\u02dc}{t}}+{R}_{e}\frac{\stackrel{\u02dc}{\sigma}}{{W}_{e}}=\frac{2\omega}{{W}_{e}}\stackrel{\u02dc}{D}\left(\stackrel{\to}{\stackrel{\u02dc}{u}}\right)\end{array}$ (1')

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 $\left\{\left({u}^{\prime},{v}^{\prime},{w}^{\prime}\right);{p}^{\prime};{\sigma}^{\prime}\right\}$ and of basic flow $\left\{\left(0,0,{W}_{b}\right),{P}_{b},{\sigma}_{b}\right\}$ . Disturbance equations are obtained by subtracting the written conservation equations for the disturbed flow $\left\{{u}^{\prime},{v}^{\prime},{W}_{b}+{w}^{\prime},{P}_{b}+{p}^{\prime},{\sigma}_{b}+{\sigma}^{\prime}\right\}$ , those satisfied by the basic flow. We obtain after calculation:

$\frac{\partial {u}^{\prime}}{\partial t}=-\left({\stackrel{\to}{W}}_{b}\cdot \nabla \right){\stackrel{\to}{u}}^{\prime}-\left({\stackrel{\to}{u}}^{\prime}\cdot \nabla \right){\stackrel{\to}{W}}_{b}-\nabla {p}^{\prime}+\frac{\left(1-\omega \right)}{{R}_{e}}\Delta {\stackrel{\to}{u}}^{\prime}+div{\sigma}^{\prime}$ (2)

$div{\stackrel{\to}{u}}^{\prime}=0$ (3)

$\frac{\partial {\sigma}^{\prime}}{\partial t}+{u}^{\prime}\frac{\partial {\sigma}_{b}}{\partial r}+{W}_{b}\frac{\partial {\sigma}^{\prime}}{\partial z}+{g}_{a}\left({\stackrel{\to}{u}}^{\prime},{\sigma}_{b}\right)+{g}_{a}\left({\stackrel{\to}{W}}_{b},{\sigma}^{\prime}\right)+\frac{{\sigma}^{\prime}}{{W}_{e}}=\frac{2\omega}{{R}_{e}{W}_{e}}D\left({\stackrel{\to}{u}}^{\prime}\right)$ (4)

where

${g}_{a}\left(\sigma ,\stackrel{\to}{u}\right)=\sigma W\left(\stackrel{\to}{u}\right)-W\left(\stackrel{\to}{u}\right)\sigma +a\left(D\left(\stackrel{\to}{u}\right)\sigma +\sigma D\left(\stackrel{\to}{u}\right)\right)$ (5)

$D\left(\stackrel{\to}{u}\right)=\frac{\left(\nabla \stackrel{\to}{u}+{\nabla}^{t}\stackrel{\to}{u}\right)}{2}$

Is the tensor of deformation rates and

$W\left(\stackrel{\to}{u}\right)=\frac{\left(\nabla \stackrel{\to}{u}-{\nabla}^{t}\stackrel{\to}{u}\right)}{2}$

is the vorticity tensor. $a$ is defined as a parameter which value is between −1 and +1.

Considering in addition that the elastic contribution of the extra stress ${\sigma}^{\prime}$ is a perturbation of the Newtonian one ${{\sigma}^{\prime}}_{s}$ which is equal to $\frac{2\omega}{{R}_{e}}D\left({\stackrel{\to}{u}}^{\prime}\right)$ , it comes:

$\frac{\partial {u}^{\prime}}{\partial t}=-\left({\stackrel{\to}{W}}_{b}\cdot \nabla \right){\stackrel{\to}{u}}^{\prime}-\left({\stackrel{\to}{u}}^{\prime}\cdot \nabla \right){\stackrel{\to}{W}}_{b}-\nabla {p}^{\prime}+\frac{\omega}{{R}_{e}}\Delta {\stackrel{\to}{u}}^{\prime}+div{{\sigma}^{\prime}}_{e}$ (6)

$\frac{\partial {{\sigma}^{\prime}}_{e}}{\partial t}+{W}_{b}^{s}\frac{\partial {{\sigma}^{\prime}}_{e}}{\partial z}+\frac{{{\sigma}^{\prime}}_{e}}{{W}_{e}}+{{\sigma}^{\prime}}_{e}\nabla \left(\stackrel{\to}{{W}_{b}^{s}}\right)=-\frac{2\omega}{{R}_{e}}\left\{{W}_{b}^{e}{\partial}_{z}D\left({\stackrel{\to}{u}}^{\prime}\right)+D\left({\stackrel{\to}{u}}^{\prime}\right)\nabla \left(\stackrel{\to}{{W}_{b}^{e}}\right)\right\}$ (7)

By $\stackrel{\to}{u}$ , we define the none-dimensional velocity field of the fluid whose radial, azimuth and axial components are given in the order by the triple $\left(u,v,w\right)$ , $\sigma $ is an additional term called extra stress and $\omega $ is the parameter of delay defined by

$\omega =1-\frac{{\lambda}^{\prime}}{\lambda}$ (8)

The system is closed with the following boundary and initial conditions:

$u=0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{if}\text{\hspace{0.17em}}r=1$ (9)

where $\frac{R{W}_{0}}{\nu}=Re$ is defined as the Reynolds number with $\nu =\frac{\eta}{\rho}$ the kinematic viscosity of the fluid, $We=\frac{\lambda {W}_{0}}{R}$ 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 $\stackrel{\to}{u}$ in the following form.

${u}_{s}\left(r,\theta ,z;t\right)=\underset{l=-L}{\overset{L}{{\displaystyle \sum}}}\underset{n=-N}{\overset{N}{{\displaystyle \sum}}}\underset{m=0}{\overset{M}{{\displaystyle \sum}}}{\alpha}_{mnl}\left(t\right){\varphi}_{mnl}\left(r\right){\text{e}}^{i\left(n\theta +l{k}_{0}z\right)}$ (10)

This approximation was used by Meseguer & Trefethen (2001) [14] .

The continuity equation leads to a linear dependence between the three components of ${u}_{s}\left(r,\theta ,z\right)$ leading to a system with two degrees of freedom. Therefore we note:

${u}_{s}\left(r,\theta ,z;t\right)=\underset{l=-L}{\overset{L}{{\displaystyle \sum}}}\underset{n=-N}{\overset{N}{{\displaystyle \sum}}}\underset{m=0}{\overset{M}{{\displaystyle \sum}}}\left[{\alpha}_{mnl}^{\left(j\right)}\left(t\right){\varphi}_{mnl}^{\left(j\right)}\left(r\right)\right]{\text{e}}^{i\left(n\theta +l{k}_{0}z\right)}$ (11)

For long periods, this approximation can be written in the form:

${u}_{s}\left(r,\theta ,z,t\right)=\underset{l=-L}{\overset{L}{{\displaystyle \sum}}}\underset{n=-N}{\overset{N}{{\displaystyle \sum}}}\underset{m=0}{\overset{M}{{\displaystyle \sum}}}\left[{\text{e}}^{\chi t}{\varphi}_{mnl}^{\left(j\right)}\left(r\right)\right]{\text{e}}^{i\left(n\theta +l{k}_{0}z\right)}$ (12)

For $j=\left(1,2\right)$ the substitution of Equation (12) in the Equation (6) the problem result in a system of ordinary differential equations of coefficients ${\alpha}_{m}^{\left(j\right)}\left(t\right)$ . This substitution is followed by a projection based on the scalar product:

$\left(F,g\right)=\underset{-1}{\overset{1}{{\displaystyle \int}}}w\left(z\right){F}^{*}\cdot g\text{d}z$

where ${F}^{*}$ designates the conjugated function of F.

The basic functions will be chosen so that integrand is even, which will result in:

$\underset{-1}{\overset{1}{{\displaystyle \int}}}G\left(r\right)\text{d}r=2\underset{0}{\overset{1}{{\displaystyle \int}}}G\left(r\right)\text{d}r$ (13)

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:

$\left(F,g\right)=\underset{0}{\overset{1}{{\displaystyle \int}}}w\left(r\right){F}^{*}\cdot g\text{d}r\approx \underset{j=0}{\overset{M}{{\displaystyle \sum}}}\text{\hspace{0.05em}}w\left({r}_{j}\right){F}^{*}\left({r}_{j}\right)g\left({r}_{i}\right)$ (14)

$F={\varphi}_{m,n,l}^{\left(1,2\right)}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{and}\text{\hspace{0.17em}}\text{\hspace{0.17em}}g={\psi}_{m,n,l}^{\left(1,2\right)}$ (15)

We define,

${h}_{m}\left(r\right)=\left(1-{r}^{2}\right){T}_{2m};\text{\hspace{0.17em}}\text{\hspace{0.17em}}{g}_{m}\left(r\right)={\left(1-{r}^{2}\right)}^{2}{T}_{2m}$ (16)

${T}_{2m}$ is the Chebyshev polynomial defined by ${T}_{2m}=\mathrm{cos}\left(2m\theta \right)$

The basic functions and tests were proposed by Leonard & Wray 1982 [10] then used by Meseguer & Trefethen 2001 [14] and are defined as follows:

Basic functions

1^{st} case:
$n=0$

${\varphi}_{m,l,n}^{\left(1\right)}=\left(\begin{array}{c}0\\ \begin{array}{c}r{h}_{m}\left(r\right)\\ 0\end{array}\end{array}\right);\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\varphi}_{m,l,n}^{\left(2\right)}=\left(\begin{array}{c}i{k}_{0}lr{g}_{m}\left(r\right)\\ \begin{array}{c}0\\ \{\begin{array}{l}{D}_{+}\left[r{g}_{m}\left(r\right)\right]\text{\hspace{0.17em}}\text{\hspace{0.17em}}si\text{\hspace{0.17em}}l\ne 0\\ {h}_{m}\left(r\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}si\text{\hspace{0.17em}}l=0\end{array}\end{array}\end{array}\right)$ (17)

2^{nd} case:
$n\ne 0$

${\varphi}_{m,l,n}^{\left(1\right)}=\left(\begin{array}{c}-in{r}^{\delta -1}{g}_{m}\left(r\right)\\ \begin{array}{c}D\left[{r}^{\delta}{g}_{m}\left(r\right)\right]\\ 0\end{array}\end{array}\right);\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\varphi}_{m,l,n}^{\left(2\right)}=\left(\begin{array}{c}0\\ \begin{array}{c}-i{k}_{0}l{r}^{\delta +1}{h}_{m}\left(r\right)\\ in{r}^{\delta}{h}_{m}\left(r\right)\end{array}\end{array}\right)$ (18)

$\delta =\{\begin{array}{l}\text{2}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{if}\text{\hspace{0.17em}}n\text{\hspace{0.17em}}\text{is}\text{\hspace{0.17em}}\text{even}\\ \text{1}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{if}\text{\hspace{0.17em}}n\text{\hspace{0.17em}}\text{is}\text{\hspace{0.17em}}\text{odd}\end{array}$ (19)

Test functions

1^{st}case:
$n=0$

${\psi}_{m,l,n}^{\left(1\right)}=\left(\begin{array}{c}0\\ \begin{array}{c}{h}_{m}\left(r\right)\\ 0\end{array}\end{array}\right)\frac{1}{\sqrt{1-{r}^{2}}}$ (20)

${\psi}_{m,l,n}^{\left(2\right)}=\left(\begin{array}{c}-i{k}_{0}l{r}^{2}{g}_{m}\left(r\right)\\ \begin{array}{c}0\\ \{\begin{array}{l}{D}_{+}\left[{r}^{2}{g}_{m}\left(r\right)+{r}^{3}{h}_{m}(r)\right]\text{\hspace{0.17em}}\text{\hspace{0.17em}}si\text{\hspace{0.17em}}l\ne 0\\ r{h}_{m}\left(r\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}si\text{\hspace{0.17em}}l=0\end{array}\end{array}\end{array}\right)\frac{1}{\sqrt{1-{r}^{2}}}$ (21)

2^{nd} case:
$n\ne 0$

${\psi}_{m,l,n}^{\left(1\right)}=\left(\begin{array}{c}-in{r}^{\beta}{g}_{m}\left(r\right)\\ \begin{array}{c}D\left[{r}^{\beta +1}{g}_{m}\left(r\right)+{r}^{\beta +2}{h}_{m}\left(r\right)\right]\\ 0\end{array}\end{array}\right)\frac{1}{\sqrt{1-{r}^{2}}}$ (22)

${\psi}_{m,l,n}^{\left(2\right)}=\left(\begin{array}{c}0\\ \begin{array}{c}i{k}_{0}l{r}^{\beta +2}{h}_{m}\left(r\right)\\ \{\begin{array}{l}\left[-in{r}^{\beta +1}{h}_{m}\left(r\right)\right]\text{\hspace{0.17em}}\text{\hspace{0.17em}}si\text{\hspace{0.17em}}l\ne 0\\ {r}^{1-\beta}{h}_{m}\left(r\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}si\text{\hspace{0.17em}}l=0\end{array}\end{array}\end{array}\right)\frac{1}{\sqrt{1-{r}^{2}}}$ (23)

$\beta =\{\begin{array}{l}0\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{if}\text{\hspace{0.17em}}n\text{\hspace{0.17em}}\text{is}\text{\hspace{0.17em}}\text{even}\\ \text{1}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{if}\text{\hspace{0.17em}}n\text{\hspace{0.17em}}\text{is}\text{\hspace{0.17em}}\text{odd}\end{array}$ (24)

Using the Fourier representation, the continuity equation becomes

$Du\left(r\right)+\frac{in}{r}v\left(r\right)+i{\lambda}_{0}w\left(r\right)=0$ (25)

${D}_{+}=D+\frac{1}{r}$ (26)

4. Numerical Implementation

Let

$\frac{\partial \stackrel{\to}{u}}{\partial t}=-\left(\stackrel{\to}{u}\cdot \stackrel{\to}{\nabla}\right)\stackrel{\to}{{W}_{b}}-\left(\stackrel{\to}{{W}_{b}}\cdot \stackrel{\to}{\nabla}\right)\stackrel{\to}{u}-\stackrel{\to}{\nabla}p+\frac{1}{{R}_{e}}\Delta \stackrel{\to}{u}+\stackrel{\to}{\nabla}\cdot {{\sigma}^{\prime}}_{e}$ (27)

Replacing $\stackrel{\to}{u}$ by the basic funct:ions and projecting on the basis of the test functions, it comes:

$\left(\frac{\partial {u}_{s}}{\partial t},\psi \right)=\left(-\left({u}_{s}\cdot \stackrel{\to}{\nabla}\right)\stackrel{\to}{{W}_{b}}-\left(\stackrel{\to}{{W}_{b}}\cdot \stackrel{\to}{\nabla}\right){u}_{s}-\stackrel{\to}{\nabla}p+\frac{1}{{R}_{e}}\Delta {u}_{s}+\stackrel{\to}{\nabla}\cdot {{\sigma}^{\prime}}_{e},\psi \right)$ (28)

The numerical method chosen is adapted from work in the early 1980 by Leonard & Wray [10] .

With this Petrov-Galerkin procedure that was then used by Meseguer & Trefethen 2000 [15] and Meseguer & Mellibovesky 2007 [16] , we obtain:

$\begin{array}{l}\left(\begin{array}{cc}\left({\varphi}_{m,n,l}^{1},{\psi}_{m,n,l}^{1}\right)& \left({\varphi}_{m,n,l}^{2},{\psi}_{m,n,l}^{1}\right)\\ \left({\varphi}_{m,n,l}^{1},{\psi}_{m,n,l}^{2}\right)& \left({\varphi}_{m,n,l}^{2},{\psi}_{m,n,l}^{2}\right)\end{array}\right)\left(\begin{array}{c}{\stackrel{\dot{}}{\alpha}}_{m,n,l}^{1}\\ {\stackrel{\dot{}}{\alpha}}_{m,n,l}^{2}\end{array}\right)\\ =\frac{1}{{R}_{e}}\left(\begin{array}{cc}\left(\Delta {\varphi}_{m,n,l}^{1},{\psi}_{m,n,l}^{1}\right)& \left(\Delta {\varphi}_{m,n,l}^{2},{\psi}_{m,n,l}^{1}\right)\\ \left(\Delta {\varphi}_{m,n,l}^{1},{\psi}_{m,n,l}^{2}\right)& \left(\Delta {\varphi}_{m,n,l}^{2},{\psi}_{m,n,l}^{2}\right)\end{array}\right)\left(\begin{array}{c}{\alpha}_{m,n,l}^{1}\\ {\alpha}_{m,n,l}^{2}\end{array}\right)+\left(\begin{array}{c}{c}^{1}\\ {c}^{2}\end{array}\right)+\left(\begin{array}{c}{d}^{1}\\ {d}^{2}\end{array}\right)\end{array}$ (29)

where

${c}^{1,2}=\left(-\left({u}_{s}\cdot \stackrel{\to}{\nabla}\right)\stackrel{\to}{{W}_{b}}-\left(\stackrel{\to}{{W}_{b}}\cdot \stackrel{\to}{\nabla}\right){u}_{s},{\psi}_{m,n,l}^{1,2}\right);\text{\hspace{0.17em}}\text{\hspace{0.17em}}{d}^{1,2}=\left(\stackrel{\to}{\nabla}\cdot {{\sigma}^{\prime}}_{e},{\psi}_{m,n,l}^{1,2}\right)$ (30)

The matrices ${c}^{1,2}$ and ${d}^{1,2}$ derive from the projection of the nonlinear terms and can be calculated with a pseudo-spectral method by Fast Fourier Transform (FFT).

By posing:

$A=\frac{1}{{R}_{e}}\left(\begin{array}{cc}\left(\text{\Delta}{\varphi}_{m,n,l}^{1},{\psi}_{m,n,l}^{1}\right)& \left(\text{\Delta}{\varphi}_{m,n,l}^{2},{\psi}_{m,n,l}^{1}\right)\\ \left(\text{\Delta}{\varphi}_{m,n,l}^{1},{\psi}_{m,n,l}^{2}\right)& \left(\text{\Delta}{\varphi}_{m,n,l}^{2},{\psi}_{m,n,l}^{2}\right)\end{array}\right)\left(\begin{array}{c}{\alpha}_{m,n,l}^{1}\\ {\alpha}_{m,n,l}^{2}\end{array}\right)+\left(\begin{array}{c}{c}^{1}\\ {c}^{2}\end{array}\right)+\left(\begin{array}{c}{d}^{1}\\ {d}^{2}\end{array}\right)$ (31)

And

$B=\left(\begin{array}{cc}\left({\varphi}_{m,n,l}^{1},{\psi}_{m,n,l}^{1}\right)& \left({\varphi}_{m,n,l}^{2},{\psi}_{m,n,l}^{1}\right)\\ \left({\varphi}_{m,n,l}^{1},{\psi}_{m,n,l}^{2}\right)& \left({\varphi}_{m,n,l}^{2},{\psi}_{m,n,l}^{2}\right)\end{array}\right)$ (32)

The problem obtained is a problem with the generalized eigenvalues $\left(Av=B\chi v\right)$ that we will solve numerically with QZ algorithm used by J. P. Berlioz [17] . To implement this method, we will use Housholder’s unitary reflection matrices and Givens rotation matrices (A. Quarteroni, R. Sacco and F. Saleri [18] ) and (L. Amodei and J.P. Dedieu [19] ).

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
$\left(u=0\right)$ . On the other hand, the Figure 1 clearly shows that the fluid is Newtonian
$\left(\omega =0\right)$ or viscoelastic
$\left(\omega >0\right)$ , the eigenvalues are purely real and negative that it’s worth the Reynolds’ number R_{e}. For both types of fluid, the flow remains stable whatever the Reynolds number.

5.2. Case of Homogeneous Disturbance (n ≠ 0; k_{0} = 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 10^{7} (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
${R}_{e}=5\times {10}^{3}$ . 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 (
${R}_{e}=5\times {10}^{3}$ ).

5.3. Case of an Axisymmetric Disturbance (n = 0; k_{0} ≠ 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 $\left(n=0,{k}_{0}=0\right)$ , case of Newtonian fluid (in red), case of viscoelastic fluid (black), with Reynolds’ number ${R}_{e}=5\times {10}^{3}$ .

(a) (b) (c)

Figure 2. (a) Spectrum of the eigenvalues of a homogeneous disturbance along the axis with $n=1;{k}_{0}=0$ , case of a Newtonian fluid (black), case of a viscoelastic fluid (red) with Reynolds number ${R}_{e}=5\times {10}^{3}$ and Weissenberg’s number $We={10}^{-2}$ ; (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 $\left(\omega =0\right)$ ; (c) Spectrum of the eigenvalues of a homogeneous disturbance along the axis with $n=1;{k}_{0}=0$ for different values of the Reynolds number: case of a viscoelastic fluid with $We={10}^{-3}$ .

for one-dimensional and homogeneous flows. Indeed, the term convection $\left(\stackrel{\to}{u}\cdot \nabla \right)\stackrel{\to}{u}$ , 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

(a)(b)

Figure 3. (a) Spectrum of the eigenvalues of an axisymmetric disturbance with $n=0;{k}_{0}=0.1$ , case of the Newtonian fluid (black), case of the viscoelastic fluid (red) with ${R}_{e}=5\times {10}^{3}$ and $We={10}^{-2}$ ; (b) Evolution of the most unstable eigenvalue according to the delay parameter $\omega $ with Reynolds’ number $Re={10}^{4}$ and Weissenberg’s number $We={10}^{-2}$ ; (b) shows that the one-dimensional flow of a viscoelastic fluid with a Reynolds number $Re={10}^{4}$ 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 $\omega =0.25$ for ${k}_{0}=0.1$ . However at beyond ${k}_{0}=0.1$ , this threshold of instability is reached from $\omega =0.2$ .

the eigenvalues are not negative and would indicate that the fluid flows unstable for the same Reynolds number ${R}_{e}=5\times {10}^{3}$ (Figure 4).

Figure 4. Spectrum of the eigenvalues of a three-dimensional perturbation $n=1;{k}_{0}=0.1$ case of the Newtonian fluid (black), case of viscoelastic fluid (red) with Reynolds number $Re=5\times {10}^{3}$ and Weissenberg’s number $We={10}^{-2}$ .

6. Conclusion

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 $We={10}^{-2}$ and for a delay parameter ( $\omega =0.1$ ) is unstable at from ${R}_{e}=5\times {10}^{3}$ 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.

Appendix Code

The algorithm QZ makes possible to determine the eigenvalues and optionally the eigenvectors of the problem with generalized eigenvalues $\left(Av=\chi Bv\right)$ . 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 $\left(A,B\right)$ , the double factorization $\{\begin{array}{l}QA=R\\ \left(QB\right)Z=T\end{array}$ with Q unitary matrix, Z upper Hessenberg, R and T upper triangular. We will note later on $\{\begin{array}{l}QAZ=\stackrel{\u02dc}{A}\\ QBZ=\stackrel{\u02dc}{B}\end{array}$ we can construct a unitary decomposing matrix A, as a product of

$n-1$ matrices of elementary rotation: $Q={Q}^{\left(n-1\right)}\cdots {Q}^{\left(2\right)}{Q}^{\left(1\right)}$ where each ${Q}^{\left(i\right)}$ is a complex rotation matrix. ${Q}^{\left(i\right)}$ modifies the first and second lines of B possibly creating a non-zero element in position $\left(2,1\right)$ which will have to be cancelled by a complex rotation matrix ${Z}_{1}$ modifying only the first and the second column of ${Q}^{\left(1\right)}B$ . Note ${B}^{\left(1\right)}=B$ and ${B}^{\left(2\right)}={Q}^{\left(1\right)}{B}^{\left(1\right)}{Z}^{\left(1\right)}$ .

${B}^{\left(2\right)}$ is triangular superior.

Let ${B}^{\left(k\right)}$ be the upper triangular, let ${Q}^{\left(k\right)}$ be the rotation matrix modifying only the kth and (k + 1)th row of ${B}^{\left(k\right)}$ possibly creating a non-zero element in position $\left(k+1,k\right)$ which will have to be canceled by a matrix complex rotation ${Z}^{\left(k\right)}$ only modifying the kth and (k + 1)th columns of ${Q}^{\left(k\right)}{B}^{\left(k\right)}$ .

Let’s say ${B}^{\left(k+1\right)}$ is upper triangular, the process stops for $k=n-1$ and we can write $\stackrel{\u02dc}{B}=QBZ={B}^{\left(n\right)}$ , $Z={Z}^{\left(1\right)}\times {Z}^{\left(2\right)}\times \cdots \times {Z}^{\left(n-1\right)}$ .

If B is invertible, it is obvious that Q realizes a unitary and triangular decomposition of $A{B}^{-1}$ and ${Z}^{*}$ a unitary decomposition of ${B}^{-1}A$ , in fact:

${Z}^{*}{B}^{-1}A={Z}^{*}{B}^{-1}{Q}^{*}QA=\stackrel{\u02dc}{B}QA$

$\stackrel{\u02dc}{B}$ and QA are upper triangular. We obtain the eigenvalues by calculating the ratios of diagonal elements of QA and $\stackrel{\u02dc}{B}$ .

Nomenclature

Greek letters

θ: Azimuthal coordinate

λ: Relaxation time of the fluid

${\lambda}^{\prime}$ : delay time

ν: Kinematic viscosity of the fluid

ρ: Density of the fluid

σ: Tensor of extra stress

${\sigma}^{\prime}$ : Disturbance of extra-stress

${\sigma}^{s}$ : Newtonian contribution of the disturbance of the extra-stress

${\sigma}^{e}$ : Elastic contribution of the disturbance of the extra-stress

η: Dynamic viscosity of the fluid

ϕ: Basis function

ψ: Test function

ω: delay parameter

Latin letters

C_{i}: imaginary part of the eigenvalue

C_{r}: real part of the eigenvalue

f: density force

l: axial mode

k_{0}: Axial wave number

n: Azimuthal mode

r: radial coordinate

R: radius of the cylinder

Re: Reynolds number

W_{0}: Maximum speed of the base flow, it has for direction that of the axis of the cylinder

We: Weissenberg’s number

W_{b}: Axial velocity of the base flow

${W}_{b}^{s}$ : Newtonian contribution of basic flow velocity

${W}_{b}^{e}$ : Elastic contribution of basic flow velocity

z: axial coordinate

References

[1] Meseguer, A. and Trefethen, L.N. (2003) Linearized Pipe Flow to Reynolds Number 107. Journal of Computational Physics, 186, 178-197.

https://doi.org/10.1016/S0021-9991(03)00029-9

[2] Esmael, A. (2008) Transition to Turbulence for a Flow through Fluid in a Cylindrical Pipe. PhD thesis, INPL, Nancy.

[3] Carranza, L. (2012) Laminar-Turbulent Transition in Cylindrical Pipe for a Non-Newtonian Fluid. PhD Thesis, INPL, Nancy.

[4] Oldroyd, J.G. (1950) On the Formulation of Rheological Equations of State. Proceedings of the Royal Society A, 200, 523-541.

https://doi.org/10.1098/rspa.1950.0035

[5] Bergstrom, L. and Brosa, U. (1988) Onset of Turbulence in a Pipe. Zeitschrift für Naturforschung, 43a, 696.

[6] 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.

[7] 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.

[8] Zikanov, O.Yu. (1996) On the Instability of Pipe Poiseuille Flow. Physics of Fluids, 8, 2923.

https://doi.org/10.1063/1.869071

[9] Bergstrom, L. (1997) Optimal Growth of Small Disturbances in Pipe Poiseuille Flow. Physics of Fluids, 9, 1043.

[10] 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.

https://doi.org/10.1007/3-540-11948-5_40

[11] Schmid, P.J. and Henningson, D.S. (1994) Optimal Energy Grothin Hagen-Poiseuille Flow. Journal of Fluid Mechanics, 277, 197.

https://doi.org/10.1017/S0022112094002739

[12] 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.

https://doi.org/10.1016/S0045-7825(98)00364-8

[13] Chupin, L. (2003) Contribution to the Study of Mixtures of Viscoelastic Fluids. Doctoral Thesis, The Mathematical Institute of Bordeaux, 177 p.

[14] Meseguer, A. and Trefethen, L.N. (2001) A Spectral Petrov-Galerkin Formulation for Pipe Flow II: Nonlinear Transitional Stages. Numerical Analysis Group, Oxford University, Rep. 01/19.

[15] 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.

[16] 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.

https://doi.org/10.1016/j.apnum.2006.09.002

[17] Berlioz, J.P. (1979) A propos de l’algorithme QZ. RAIRO-Analyse Numérique, 13, 21-30.

https://doi.org/10.1051/m2an/1979130100211

[18] Quarteroni, A., Sacco, R. and Saleri, F. (2004) Méthodes Numériques: Algorithmes, analyse et applications. Springer-Verlag, Italia, Milano.

[19] Amodei, L. and Dedieu, J.P. (2008) Analyse numérique matricielle: Cours et exercices corrigés détaillés. Dunod, Paris 2008, SMAI (Société de Mathématiques Appliquées et Industrielle).