ABSTRACT

A Fourier-Chebyshev Petrov-Galerkin spectral method is described for high accuracy computation of linearized dynamics for flow in a circular pipe. The code used here is based on solenoidal velocity variables and is written in FORTRAN. Systematic studies are presented of the dependence of eigenvalues and other quantities on the axial and azimuthal wave numbers; the Reynolds’ number of up to 10⁷ and the Weissenberg’s number that is considered lower here. The flow will be considered stable if all the real parts of the eigenvalues obtained are negative and unstable if only one of these values is positive.

Keywords:

Viscoelastic Fluids, Oldroyd-B Model, Linear Stability, Petrov-Galerkin, Generalized Eigenvalue Problem

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 \overrightarrow{u}}{\partial t}+\left(\overrightarrow{u}\cdot \nabla \right)\overrightarrow{u}\right)=-\nabla P+\overrightarrow{f}+div\left\{2\eta \left(1-\omega \right)D\left(\overrightarrow{u}\right)\right\}+div\sigma \\ \frac{\partial \rho }{\partial t}+div\left(\rho \overrightarrow{u}\right)=0\\ \frac{D\sigma }{Dt}+\frac{\sigma }{\lambda }=\frac{2\eta \omega }{\lambda }D\left(\overrightarrow{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:

$\tilde{p}=\frac{p}{\rho W_0^2};\tilde{\sigma }=\frac{\sigma }{\rho W_0^2};\overrightarrow{\tilde{u}}=\frac{\overrightarrow{u}}{W_0};\tilde{r}=\frac{r}{R};\tilde{z}=\frac{z}{R};\tilde{t}=\frac{t{W}_0}{R},\tilde{\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 \overrightarrow{\tilde{u}}}{\partial t}+\left(\overrightarrow{\tilde{u}}\cdot \overrightarrow{\tilde{\nabla }}\right)\overrightarrow{\tilde{u}}\right\}=-\overrightarrow{\tilde{\nabla }}\tilde{p}+\frac{\left(1-\omega \right)}{R_e}\tilde{\Delta }\overrightarrow{\tilde{u}}+\overrightarrow{\tilde{\nabla }}\cdot \tilde{\sigma }\\ \overrightarrow{\tilde{\nabla }}\cdot \overrightarrow{\tilde{u}}=0\\ R_e\frac{D\tilde{\sigma }}{D\tilde{t}}+R_e\frac{\tilde{\sigma }}{W_e}=\frac{2\omega }{W_e}\tilde{D}\left(\overrightarrow{\tilde{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({\overrightarrow{W}}_b\cdot \nabla \right){\overrightarrow{u}}^{\prime }-\left({\overrightarrow{u}}^{\prime }\cdot \nabla \right){\overrightarrow{W}}_b-\nabla p^{\prime }+\frac{\left(1-\omega \right)}{R_e}\Delta {\overrightarrow{u}}^{\prime }+div{\sigma }^{\prime }$ (2)

$div{\overrightarrow{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({\overrightarrow{u}}^{\prime },{\sigma }_b\right)+g_a\left({\overrightarrow{W}}_b,{\sigma }^{\prime }\right)+\frac{{\sigma }^{\prime }}{W_e}=\frac{2\omega }{R_e{W}_e}D\left({\overrightarrow{u}}^{\prime }\right)$ (4)

where

$g_a\left(\sigma ,\overrightarrow{u}\right)=\sigma W\left(\overrightarrow{u}\right)-W\left(\overrightarrow{u}\right)\sigma +a\left(D\left(\overrightarrow{u}\right)\sigma +\sigma D\left(\overrightarrow{u}\right)\right)$ (5)

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

Is the tensor of deformation rates and

$W\left(\overrightarrow{u}\right)=\frac{\left(\nabla \overrightarrow{u}-{\nabla }^t\overrightarrow{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({\overrightarrow{u}}^{\prime }\right)$ , it comes:

$\frac{\partial u^{\prime }}{\partial t}=-\left({\overrightarrow{W}}_b\cdot \nabla \right){\overrightarrow{u}}^{\prime }-\left({\overrightarrow{u}}^{\prime }\cdot \nabla \right){\overrightarrow{W}}_b-\nabla p^{\prime }+\frac{\omega }{R_e}\Delta {\overrightarrow{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(\overrightarrow{W_b^s}\right)=-\frac{2\omega }{R_e}\left\{W_b^e{\partial }_{z}D\left({\overrightarrow{u}}^{\prime }\right)+D\left({\overrightarrow{u}}^{\prime }\right)\nabla \left(\overrightarrow{W_b^e}\right)\right\}$ (7)

By $\overrightarrow{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{if}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 $\overrightarrow{u}$ in the following form.

$u_s\left(r,\theta ,z;t\right)={\displaystyle \sum }_{l=-L}^L{\displaystyle \sum }_{n=-N}^N{\displaystyle \sum }_{m=0}^M{\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)={\displaystyle \sum }_{l=-L}^L{\displaystyle \sum }_{n=-N}^N{\displaystyle \sum }_{m=0}^M\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)={\displaystyle \sum }_{l=-L}^L{\displaystyle \sum }_{n=-N}^N{\displaystyle \sum }_{m=0}^M\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 {\displaystyle \sum }_{j=0}^{M}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{and}g={\psi }_{m,n,l}^{\left(1,2\right)}$ (15)

We define,

$h_m\left(r\right)=\left(1-r^2\right)T_{2m};g_m\left(r\right)={\left(1-r^2\right)}^2{T}_{2m}$ (16)

$T_{2m}$ is the Chebyshev polynomial defined by $T_{2m}=\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);{\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]sil\ne 0\\ h_m\left(r\right)sil=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);{\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{if}n\text{is}\text{even}\\ \text{1}\text{if}n\text{is}\text{odd}\end{array}$ (19)

Test functions

1^stcase: $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]sil\ne 0\\ r{h}_m\left(r\right)sil=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]sil\ne 0\\ r^{1-\beta }{h}_m\left(r\right)sil=0\end{array}\end{array}\end{array}\right)\frac{1}{\sqrt{1-r^2}}$ (23)

$\beta =\{\begin{array}{l}0\text{if}n\text{is}\text{even}\\ \text{1}\text{if}n\text{is}\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 \overrightarrow{u}}{\partial t}=-\left(\overrightarrow{u}\cdot \overrightarrow{\nabla }\right)\overrightarrow{W_b}-\left(\overrightarrow{W_b}\cdot \overrightarrow{\nabla }\right)\overrightarrow{u}-\overrightarrow{\nabla }p+\frac{1}{R_e}\Delta \overrightarrow{u}+\overrightarrow{\nabla }\cdot {{\sigma }^{\prime }}_e$ (27)

Replacing $\overrightarrow{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 \overrightarrow{\nabla }\right)\overrightarrow{W_b}-\left(\overrightarrow{W_b}\cdot \overrightarrow{\nabla }\right)u_s-\overrightarrow{\nabla }p+\frac{1}{R_e}\Delta u_s+\overrightarrow{\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}{\dot{\alpha }}_{m,n,l}^1\\ {\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 \overrightarrow{\nabla }\right)\overrightarrow{W_b}-\left(\overrightarrow{W_b}\cdot \overrightarrow{\nabla }\right)u_s,{\psi }_{m,n,l}^{1,2}\right);d^{1,2}=\left(\overrightarrow{\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{Δ}{\varphi }_{m,n,l}^1,{\psi }_{m,n,l}^1\right)& \left(\text{Δ}{\varphi }_{m,n,l}^2,{\psi }_{m,n,l}^1\right)\\ \left(\text{Δ}{\varphi }_{m,n,l}^1,{\psi }_{m,n,l}^2\right)& \left(\text{Δ}{\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)

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⁷ (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)

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

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

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).

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.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

Cite this paper

Kama, I., Sow, M.L., Kpode, K. and Mbow, C. (2018) Effects of the Reynolds and Weissenberg’s Numbers on the Stability of Linear Pipe Flow of Viscous Fluid. Modern Mechanical Engineering, 8, 264-277. https://doi.org/10.4236/mme.2018.84018

References

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=\tilde{A}\\ QBZ=\tilde{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 $\tilde{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=\tilde{B}QA$

$\tilde{B}$ and QA are upper triangular. We obtain the eigenvalues by calculating the ratios of diagonal elements of QA and $\tilde{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₀: Axial wave number

n: Azimuthal mode

r: radial coordinate

R: radius of the cylinder

Re: Reynolds number

W₀: 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