Linear Stability Analysis of Viscoelastic Flows

(1)

May 9, 2012

Linear Stability Analysis of

Viscoelastic Flows

Mengqi Zhang

Thesis submitted to Royal Institute of Technology

for Master’s degree

Supervisor:

Prof. Luca Brandt

(2)

Abstract

The goal of present work is to investigate the instability of jet flow, mixing layer and Poiseuille flow of viscoelastic fluids. According to Boffeta et al. [4], small elastic effect in Kolmogorov flow will result in increasing critical Re for hydrodynamic instability, however, at high elasticity a new instability of elastic nature occurs even at vanishing Re. This thesis aims to test this result in jet and mixing layer. In addition, linear stability analysis (modal and non-modal) of viscoelastic Poiseuille flow is carried out to understand the elastic effects on the flows of both Oldroyd-B and FENE-P model fluids. Energy analysis is used to reveal the instability mechanism.

(3)

Acknowledgements

First of all, I would like to thank my supervisor Prof. Luca Brandt for offer-ing me the opportunity to do this project and supervisoffer-ing me with his great enthusiasm and knowledge on non-Newtonian flow. Without him, I would not appreciate the beauty of research.

I would also like to thank Prof. Tamer Zaki at Imperial College London for reviewing the thesis drafts in detail during the writing and offering enlightening suggestions.

Finally, I want to thank the Ph.D. students in the department for creating such a good atmosphere, especially Iman for his patience in discussion and jokes in leisure time, and Lailai for help.

(4)

1 Introduction

1.1 Viscoelastic Instability

Flow viscoelastic instability is due to the elasticity of additives in a viscous solvent. It is more complex than that of its Newtonian counterpart, and can lead to a chaotic flow named as elastic turbulence [18]. Unlike the conventional turbulence state which is due to the non-linear inertial force at large Reynolds number (Re), the non-linear effect giving rise to elastic turbulence stems from stretching and recoiling of the polymer molecules, manifested at large Weis-senberg number (W ). The elastic turbulence resistance is due to the elastic stress, rather than the viscous shear stress in conventional turbulence. (Re is a non-dimensional number describing the relative importance between the in-ertial forces and the viscous forces. W is the ratio between the characteristic relaxation time of polymer molecules and the flow characteristic time.)

The research on viscoelastic instability started from the middle of last cen-tury. Giesekus [22] first reported instability of polymer solution in circular Couette flow at a Taylor number much lower than the critical value for New-tonian flow but he considered it as an unusual fluid property. Taylor number is a dimensionless ratio between the inertial forces due to fluid rotation, and the viscous forces. Petrie & Denn [34] then pointed out that this is a new in-stability mechanism, now known as elastic inin-stability. After that, researchers tried to reveal the mechanism behind this instability using polymer solutions or melts. But since these solutions are inherently shear thinning (the fluid viscos-ity will decrease with increasing rate of stress), it is very difficult to distinguish whether the observed non-Newtonian effect is brought about by shear thinning or elasticity of the polymer, or both. The advent of Boger fluid in 1977 (James 2009 [22]) solved this dilemma in that this fluid has constant viscosity and at the same time high elasticity, the two ideal properties for rheologists. Boger thought that since shear thinning results from high concentration of the poly-mer, truly dilute fluid can yield a constant viscosity fluid. The first Boger fluid was introduced as a 0.08% solution of polyacrylamide in a concentrated aque-ous sugar solution. With the Boger fluid, in principle, the experiment could be conducted with two fluids, a Boger fluid and a Newtonian fluid, therefore it is possible to differentiate flow elastic effect from viscous effect.

The viscoelastic instability in various flow regimes, such as parallel shear flow, plate-and-plate flow and Taylor-Couette flow, is well explored. For the plane Poiseuille flow, its viscoelastic instability was first documented by Por-teous & Denn (1972 [36]). They found that in the Poiseuille flow of second order fluid, weak viscoelasticity destabilizes the flow, namely, the critical Re in the linear stability analysis decreases as the Weissenberg number increases. But Porteous, Denn and others then showed that when the W increases further (strong viscoelastic effect), the critical Re increases, that is, the viscoelasticity stabilizes the flow. These findings were discovered in the framework of classical linear stability theory, but since transition to turbulence is subcritical in channel flows, non-modal analysis could provide us with more insight (Jovanovi´c et. al.

(6)

[26][27]). Thus, probing transient growth of plane Poiseuille flow will be one goal of this thesis.

In the recent decades, the landmark elastic instability research is due to Phan-Thien in 1983 [35]. In his experiment, Phan-Thien analyzed the instabil-ity of an Oldroyd-B fluid in steady shear flow between two concentric disks. He showed that when the value of λωr is less than π

q

2

5β, the flow is stable,

other-wise unstable. λ is the characteristic relaxation time of the polymer molecules, ωris the angular speed of the rotating plate, β is the viscosity ratio µs/µ, where

µsand µ are solvent viscosity and total viscosity, respectively. The subsequent

experiments by others agreed Phan-Thien’s theory. In the early 1990s, Larson, Muller and Shaqfeh exhaustively investigated the elastic instability with iner-tialess Taylor-Couette flow of dilute polymer solutions. In their work, Larson et al. (1990 [29]) used linear stability analysis to predict the onset of purely elastic instability both experimentally and theoretically. They performed two sets of experiments, one at constant rotation rate of the outer cylinder and the other at constant stress imposed on the inner cylinder. They also formualted the eigenvalue problem and found its solution with various methods, all giving consistent results. In both cases, the onset of an instability occurs at very small Taylor number and the critical Deborah number should vary inversely with the gap. Deborah number is the ratio between the stress relaxation time and the time scale of the experimental observation.

As for viscoelastic jet and mixing layer, Rallison & Hinch (1995 [37]) inves-tigated the former and Azaiez & Homsy (1994 [1] [2]) performed linear stability analysis and nonlinear simulations on the latter. Rallison & Hinch consider submerged elastic jet charactered by parabolic flow profile at high Re. They found that the sinuous mode is fully stabilized by large elasticity, while varicose mode is partially stabilized. On the other hand, Azaiez & Homsy studied the inviscid free shear layer (high Re) with high elastic effect (high W ), so that the ratio between these two dimensionless numbers is of order unity. In their linear stability analysis, the flow modeled by Oldroyd-B model was found to be less unstable with increasing elastic effect, but not totally stable. The unstable wavelength becomes longer compared to Newtonian flow (see figure (2) in their work [?], which is also reproduced by our code in figure (6)). In their nonlinear simulation, they identified that the deficiency of Oldroyd-B model allowing in-finite polymer chain extension stems from that in the braid regions the product of W and local extensional rate are larger than 1. Their FENE-P model simu-lation results reveal that the global structure of the flow is unchanged from the Newtonian flow. But the local vorticity is intensified because of the increasing elastic polymer normal stress.

Readers are referred to the work by Shaqfeh (1996 [42]) and Larson (1972 [28]) for more discussion and history on this subject.

(7)

1.2 Drag Reduction

It has been noted that a small amount of polymer additives added in a turbulent flow would lead to a remarkable drag reduction up to 70% ∼ 80% (measured by the difference of the pressure gradient in channel flows). According to several previous works, the drag reduction is believed to be connected with the flow instability, therefore some new insight into the drag reduction mechanism could be gained when we investigate transient energy growth of the flows.

The drag reduction was first observed by Toms in 1949. Since then, nu-merous investigations are conducted to unravel the mechanism behind this phe-nomena. However, the mechanism remains poorly understood because of its complex nature. Before the direct numerical simulations, we could find two main attempts to explain the phenomenon. The first one was the time criterion (Lumley 1969 [30]), arguing that in the case of drag reduction the time scale of turbulence should be less than the characteristic relaxation time scale of the polymer molecules, i.e. λ > ν/u2

τ, where uτ is the friction velocity in the

near-wall region and ν is the kinematic viscosity. The presence of polymer additives also contributes to the extra elongational viscosity, which was proposed by Lum-ley (1973 [31]) and Hinch (1977 [20]) to account for the drag reduction by the so called ’coil -stretch’ transition. The other theory was elastic theory, appealing to the elastic memory nature of the polymer (de Gennes 1990 [11]). The main argument is that the turbulent energy cascade is suppressed by the polymer at some small scales. In the recent decades, new experiments and powerful direct numerical simulations saw further advances in this subject, discussed below.

1.2.1 Experiments

One of the important experiments was conducted by Cadot, Bonn and Douady in 1998 [6]. They demonstrated that the wall effect is essential for drag re-duction. The experiment consisted of two closed cells. In one of them, the turbulence is forced by the smooth disk, which provides boundary layer effects, while in the other one, the turbulence is created by the baffles, which solely of-fers inertial effects. Their results showed that the drag reduction is only found in the case forced by wall effect (boundary layer), indicating that the near-wall region is important in the generation of drag reduction.

Another important experiment was due to Warholic, Massah and Hanratty (1999 [45]), in which they probed extensively the turbulence statistical proper-ties in the near-wall region. In their work, they differentiated low drag reduction (LDR, < 35% ∼ 40%) and high drag reduction (HDR, 40% ∼ 70%). In both cases, it was found that the non-dimensional velocities profile changes due to the presence of the polymer molecules. In LDR, the non-dimensional mean ve-locity complies U+_{= 2.41 ln y}+_{+ B in the buffer layer with different intercept}

coefficients B for different drag reduction percentages. That larger B is found in the larger drag reduction clearly indicated that the buffer layer shifts away from the wall when drag reduction happens. The coefficient 2.41 remains unchanged, meaning that the lines in the semi-log plot has the same slope for different drag

(8)

reductions. Also, the streamwise fluctuating velocity increases in the near-wall region and its peaking value shifts outward from the wall, which indicates that the viscous sublayer thickness increases. Besides, the experiment showed a de-crease of normal and spanwise fluctuating velocities as well as Reynolds stress with increasing drag reduction percentages. In the HDR, the same trend of the aforesaid changes was observed but with some exceptions. The non-dimensional streamwise velocity in the buffer layer now does not comply a constant slope in the semi-log plot and the slope is bound to 11.7 in the maximum drag reduction (MDR) case. The normal and spanwise fluctuating velocities are suppressed much more, so is the Reynolds stress, even being zero in the near-wall region, therefore the turbulence should be sustained by other force. The authors then proposed that because the Reynolds stress is so small the conventional turbu-lence regeneration mechanism may be modified in the presence of the polymer in HDR and the interaction between the polymer and the flow may create and sustain the turbulence.

1.2.2 Numerical Simulation

Around the beginning of the new century, direct numerical simulation was widely employed to investigate the drag reduction mechanism. The numerical simula-tion gives people with powerful tools of probing and manipulating the targeted quantities in the research. Important simulations were conducted by Min et al . [32], Dubief et al. [13] Dimitropoulos et al. [12] and Dubief et al. [14], among

others. In Min’s paper, the turbulence statistics in dilute polymer solution

was obtained in agreement with the experiments, namely, the increase of vis-cous sublayer thickness, the upward shift of the buffer layer and the decrease of Reynolds shear stress. They also computed the energy budget for the mean flow, the turbulence and the polymer. Their results of energy analysis showed that the turbulence production decreases with increasing W at a certain Re. More importantly, they reported that the production of the elasticity energy decrease locally in the near-wall region while increases in the buffer layer, which leads the authors to conjecture that the polymer molecules take less energy away from the near-wall region and release more energy in the buffer and log-law layer when drag reduction happens, therefore, there is a energy transfer of elastic energy from the near-wall region to the buffer and log-law layer.

Another important direct numerical simulation was reported by Dubief et al.. They investigated the phenomenon by numerical experiment, which manipu-lates the flow conditions by deliberately enhancing or suppressing one possible generating process or quantity to see the global effect of it. Their work was meaningful in a sense that they interpreted the drag reduction mechanism in terms of the near-wall turbulence regeneration cycle proposed by Jimenez and Pinelli in [24]. According to the latter, the autonomous cycle for the turbu-lence involves that streamwise vortices extract energy from the mean flow to create alternating steaks of streamwise velocity in the near-wall region and that these streaks in turn break down to give rise to the formation of streamwise vortices, completing a cycle for regenerating the turbulence incessantly (useful

(9)

discussion can also be found in Hamilton, Kim & Waleffe [19]). Dubief et al. focused on the regime of HDR in a minimal computational domain to remove the effect of larger vortices in the turbulent core (Jim´enez & Moin [23]). The numerical experiments of manipulating the polymer stress (fx, fy and fz) lead

them to argue that fx sustains turbulence in the near-wall region by enhancing

streamwise fluctuation velocity and fy and fz reduce the drag by damping the

streamwise vortices. This result is consistent with the experimentally observed fluctuation velocity changes in the dilute polymer solution since the major con-tribution to wall-normal and spanwise fluctuating velocities comes from quasi-streamwise vortices and the quasi-streamwise fluctuating velocity is associated with the streaks, therefore the authors conclude that the polymers reduce the turbu-lence by damping the streamwise vortices, while at the same time enhance the streamwise velocity streaks in the near-wall region.

This thesis in itself is not directly related to drag reduction by polymers, but by exploring the viscoelastic instability in various flows it may serve as a small step towards our understanding of such phenomenon.

(10)

2 Problem Formulation

2.1 Polymeric Constitutive Equations

For dilute viscoelastic solution, the simplest linear constitutive equations of the polymer are the upper convected maxwell model (UCM) and the Oldroyd-B model proposed by Oldroyd in 1950s. The UCM model does not take into ac-count the hydrodynamic interaction between polymer molecules and solution. This model is the limiting case of the Oldroyd-B model, which overcomes that shortage by considering the total viscosity as the sum of Newtonian flow viscos-ity and additional viscosviscos-ity due to the presence of polymer (and it predicts no shearing thinning). However, both models are limited in some aspects. Firstly, these two models only consider one relaxation time while the polymer chains actually contain a spectrum of relaxation times. Secondly, the models fail to impose a restriction on the extensibility of the polymer molecules. It is unphys-ical to allow for an infinite extensibility though it may not always happen in these models (Min et al. 2003 [32]). The Finitely Extensible Nonlinear Elastic (FENE) model predicts finite extension of the polymer, but entails statistical closure for the restoring force, for example the Peterline closure (FENE-P).

2.1.1 Polymer stress

In our analysis, both Oldroyd-B model and FENE-P model are employed as we linearize the governing equation. In these models, the polymer chain is depicted as a linear elastic dumbbell with the two beads at each end connected with an entropic spring. The conformation tensor ¯C_ij∗ =< ¯R∗_iR¯∗_j >, where ¯R is the end-to-end vector of molecule, is the main parameter to describe the dynamics of the polymer molecule. Here, the superscript star∗represents dimensional quantity,

the overbar−indicates an undisturbed quantity and the angular brackets mean

the average over white noise. The conformation tensor is a symmetric tensor with six independent components in three dimensions and its trace tr( ¯C) is a measure of the squared polymer elongation, proportional to the elastic energy stored in the molecule. From a perspective of thermal equilibrium, the polymeric stress tensor is related to the conformation tensor by defining the stress in proportion to the deviation of the conformation tensor from its equilibrium state, e.g., in the Oldroyd-B model (Chokshi & Kumaran 2009 [9]),

¯

τ_p∗= µpH λkBT

( ¯C∗− ¯C∗eq). (1)

In this definition, H is the spring constant of the elastic dumbbell, kB is the

Boltzman constant and T is the temperature. λ, as mentioned before, is the relaxation time of the polymer molecules and µp is the additional flow viscosity

due to polymer. Non-dimensionlizing ¯C∗ by kBT /H and ¯τp∗ by µpUc/Lc (Uc

(11)

obtain ¯ τp = ¯ C − I W (Oldroyd − B model), (2) ¯ τp= f ¯C − I W (F EN E − P model), (3)

where W is the Weissenberg number defined as

W = λUc

Lc

. (4)

Weissenberg number is an important parameter in rheology measuring the elas-ticity of the polymer. High Weissenberg number implies high elastic effect. Besides, in FENE-P model,

f = 1

1 −C¯kk

L2

(5) is the Peterlin function confining the polymer extensibility to be less than L, which is the maximum extensibility of polymer molecules. ¯Ckk= ¯C11+ ¯C22+

¯

C33 is the trace of the base-state conformation tensor. Einstein’s summation

notation is used in this thesis.

2.1.2 Conformation Tensor Governing Equation

The non-dimensional governing equation for Oldroyd-B and FENE-P models with polymeric conformation tensor reads

∂ ¯C

∂t + ¯u · ∇ ¯C − ¯C · (∇¯u) − (∇¯u)

T

· ¯C = −¯τp, (6)

where ¯τp is related to conformation tensor by equations (2) and (3) of the

two models, respectively. The left-hand side of the evolution equation is the upper convective derivative acting on the conformation tensor. This type of time derivative takes into account the rotating and stretching of the coordinate system in which the fluid and polymer are observed. Since the polymer molecules are constantly stretched and rotated along with the flow, the upper convective derivative is physically relevant.

For Oldroyd-B model, an explicit expression could be derived for the con-formation tensor using the governing equations above. For FENE-P model, the conformation tensor is coupled with the base flow, which will be discussed later in the Appendix A. As for a steady-state parallel base flow ¯u = ¯u(y), it is straightforward to derive that the analytical solution for conformation tensors in Oldroyd-B model ¯ C =   2W2u¯02+ 1 W ¯u0 0 W ¯u0 ₁ ₀ 0 0 1  . (7)

(12)

¯

u0 is the normal derivative of the flow velocity, or minus spanwise vorticity −w in this particular case. Using equation (2), the stress tensor for a parallel base flow could be derived

¯ τp =   2W ¯u02 u¯0 0 ¯ u0 0 0 0 0 0  . (8)

We see that the Oldroyd-B model predicts the shear stress to be proportional to the shear rate as ¯τp12 = ¯u0 and also that the first normal stresses difference

¯

τp11− ¯τp22 is a function of the squared shear rate ¯u02 while the second normal

stresses difference ¯τp22− ¯τp33 is always zero for these simple flows.

2.2 Navier-Stokes Equation

Based on the continuum hypothesis, the unsteady incompressible flow govern-ing equations involve the continuity equation and the Navier-Stokes equation, namely, ∇ · ¯u∗= 0, (9) ρ(∂ ¯u ∗ ∂t + ¯u ∗_{· ∇¯}_u∗_{) = −∇¯}_p∗_{+ ∇ · ¯}_τ∗ t. (10) ¯

p∗ is the pressure. ¯τ_t∗ is the total deviatoric stress in the fluid. In the study of Newtonian flow, ¯τ_t∗solely consists of the viscous shear stress ¯τ_s∗and is modeled by a simple assumption that shear stress is proportional to shear rate, namely, ¯

τs∗ = µs˙γ∗, where ˙γ∗ is shear rate ∇¯u∗+ (∇¯u∗)T. For the non-Newtonian

flow, in which ¯τt∗ additionally involves the viscoelastic stress, the total stress

is modeled by ¯τt∗ = µs˙γ∗+ ¯τp∗. The polymeric constitutive models discussed

above then can be used to derive an explicit expression for the polymer stress ¯

τ_p∗. In equation (10), non-dimensionalizing the velocity ¯u∗by characteristic flow velocity Uc, the time t by characteristic flow time Lc/Uc and the pressure by

ρU2

c, we obtain in the two viscoelastic constitutive models

∂ ¯u ∂t + ¯u · ∇¯u = −∇¯p + β Re∇ 2_{u +}_¯ 1 − β Re ∇ · ¯τp. (11)

In the equations above, β, as mentioned before, is the viscosity ratio between the solvent viscosity and the total viscosity. β could be viewed as the concentration of polymers in the flow; small β implies high concentration of polymer molecules,

vice versa. Re is the Reynolds number defined as ρUcLc/µ based on the base

flow state and the total viscosity. (The subscript p in the stress tensor will be omitted in the following context).

2.3 Base Flows

Different flows configurations listed in table (1) and depicted in figure (1) are considered in this thesis.

(13)

−100 −5 0 5 10 0.2 0.4 0.6 0.8 1 y u (a) −10 −5 0 5 10 −1 −0.5 0 0.5 1 y u (b) −10 −0.5 0 0.5 1 0.2 0.4 0.6 0.8 1 y u (c)

Figure 1: Different flow profiles considered in the thesis. The abscissa is the flow domain along wall-normal direction. The ordinate is the velocity. (a) is jet flow, u = sech2_{(y). (b) is mixing layer, u = tanh(y). (c) is the plane Poiseuille}

flow, u = 1 − y2_.

Table 1: Investigated flows

Flow Profile Flow domain

J et f low u = sech2_(y) _(-10,10)

M ixing layer u = tanh(y) (-10,10)

P oiseuille f low u = 1 − y2 _(-1,1)

Base flow refers to the steady solution to Navier-Stokes equation (10). For the flow geometry, x, y, z represent the streamwise, wall-normal and spanwise directions, respectively. ¯u = (¯u, ¯v, ¯w) is the velocity components. Note that

(14)

Lc is the half length of flow domain for Poiseuille flow but for jet and mixing

layer Lc is simply set to be 1. As for Poiseuille flow using Oldroyd-B model, the

base flow remains the same as Newtonian base flow. But in FENE-P model, the parabolic base flow will be modified as shown in figure (2). One can see that the base flow profile becomes a little larger than the Newtonian base flow, meaning that due to the presence of polymers the viscoelastic fluids become easier to flow than its Newtonian counterpart if pressure gradient remains the same in the two flows. The derivation of the analytical solution to the modified Poiseuille base flow in FENE-P model is shown in Appendix A.

−10 −0.5 0 0.5 1 0.5 1 1.5 y u

Newtonian base flow FENE−P base flow

Figure 2: Newtonian base flow and the FENE-P base flow for plane Poiseuille type. Parameters are L = 60, β = 0.8, W = 10, Re = 5000

2.4 Linearization

We employ the linear stability analysis to examine the instability of jet, mixing layer and Poiseuille flow. The perturbation of the variables is introduced by decomposing the flow variables into a base component and a fluctuation compo-nent, which is ¯u = U + u, ¯p = P + p, ¯C = C + c, ¯τ = T + τ . The capital letters are the base state and the small letters are the fluctuation components. Gen-erally, the classical linear analysis of temporal instability assumes a wave-like perturbation as

φ(x, y, z, t) = ˜φ(y)eiαx+iγz−iωt. (12)

The symbol tilde ˜ represents the wall-normal shape function depending only on y. φ is the general variable φ = (u, c)T_{, α is the real-valued streamwise}

wavenumber, γ is the real-valued spanwise wavenumber and ω is the complex-valued wavespeed. In this definition, the sign of imaginary part of ω plays a pivotal role in the analysis since if positive, the perturbation will increase exponentially with time, namely, instability occurs. This wave-like disturbance intrinsically involves parallelism assumption of the flow. Thus, for jet and mixing layer, the stability analysis is local.

(15)

Linearization is applied both for flow governing equation (11) and polymer conformation tensor evolution equations (6). By substitution of the decom-position into equation (11) and subtraction of the base state terms from the equation, the linearized governing equation for the flow with the polymer con-formation term in Oldroyd-B model reads

∂u ∂t + U · ∇u + u · ∇U = −∇p + β Re∇ 2_{u +} 1 − β Re W∇ · c. (13)

In the linearized equation, the quadratic terms are neglected. By eliminating the pressure term in the three scalar equations (13), we can obtain the classical Orr-Sommerfeld (O-S) and Squire (Sq) system of equations with the conformation tensor ∂ ∂t∇ 2_{v = [−U} ∂ ∂x∇ 2_{+ U}00 ∂ ∂x + β Re∇ 4_]v+ 1 − β ReW[ ∂ ∂xj ∇2_c jy− ∂3_c jx ∂xj∂x∂y − ∂ 2_c jy ∂xj∂y2 − ∂ 2_c jz ∂xj∂x∂z ], (14) ∂ ∂tη = −U ∂ ∂xη − U 0∂v ∂z + β Re∇ 2_{η +}1 − β ReW( ∂2cxj ∂xj∂z − ∂ 2_c zj ∂xj∂x ). (15)

η here is the wall-normal vorticity ∂u_∂z −∂w

∂x. The subscript j represents one of

x, y, z directions. The boundary conditions for wall-normal velocity and normal vorticity are v = ∂

∂yv = η = 0. For the streamwise and spanwise directions, we

impose the periodic boundary condition because spectral collocation method will be used.

By substituting the decomposition of variables into the polymer conforma-tion tensor evoluconforma-tion equaconforma-tion (6), the linearized governing equaconforma-tions for the polymeric conformation tensor read in the two models, respectively,

∂c

∂t + u · ∇C + U · ∇c − c · ∇U − C · ∇u − (∇U)

T

· c − (∇u)T · C = − c

W,

(16)

∂c

∂t + u · ∇C + U · ∇c − c · ∇U − C · ∇u − (∇U)

T _{· c − (∇u)}T_{· C}

= −f

0_{C + f c}

W , (17)

where f0 is defined as (Nouar, Bottaro & Brancher 2007 [33])

f0 = ∂f ∂C11 _Bc11+ ∂f ∂C22 _Bc22+ ∂f ∂C33 _Bc33. (18)

(16)

Subscript B refers to the base state values. No boundary condition is im-posed for conformation tensor or polymer stress tensor. For FENE-P model, the linearization is more complicated than that of Oldroyd-B model because the Peterlin function (it is a function of C11, C22and C33, and therefore a function

of y) should also be linearized. Since the FENE-P model returns to Oldroyd-B model at large L, it is expected that the above equation (17) will go back to

Oldroyd-B governing equation (16) when f0 = 0 and f = 1 at large L, which is

readily seen. The code of FENE-P can also be partially verified on the base of this fact. We will discuss this more in the numerical part.

Upon substituting

v(x, y, z, t) = ˜v(y)eiαx+iγz−iωt, (19)

η(x, y, z, t) = ˜η(y)eiαx+iγz−iωt, (20)

c(x, y, z, t) = ˜c(y)eiαx+iγz−iωt, (21)

into the above equations (14), (15) and (16) or (17) and writing the resultant equations in the form of

A∂

∂tφ = Bφ, (22)

we obtain an eigenvalue problem with φ = (v, η, c11, c22, c33, c12, c13, c23)T in

three dimensions. (Note that the variables v(x, y, z, t), η(x, y, z, t), c(x, y, z, t) are complex numbers with their real parts bearing physical meanings.) The way to solve this eigenvalue problem will be discussed in the numerical method chapter.

2.5 Energy Analysis

The energy analysis provides another way of analyzing the instability mechanism of the polymeric flow (Hoda, Jovanovi´c & Kumar 2008 [21]). In the spectral space, we rewrite equation (13) in the form of

∂u ∂t = −U · ∇u − ∇p + β Re∇ · ∇u + 1 − β Re ∇ · τ − u · ∇U, (23)

The hydrodynamic perturbation kinetic energy is obtained by multiplying the complex conjugate of the perturbation velocity by the above equation

u∗_i∂ui ∂t = −u ∗ iUj ∂ui ∂xj − u∗_i ∂p ∂xi + u∗_i β Re ∂2_u i ∂2_x j + u∗_i1 − β Re ∂τij ∂xj − u∗_iuj ∂Ui ∂xj . (24) Then after taking the conjugate of the above equation and recalling that the operations of conjugate and derivative could be commuted (_dxdf)∗=df_dx∗, we get

∂u∗_i ∂t ui = −Uj ∂u∗_i ∂xj ui− ∂p∗ ∂xi ui+ β Re ∂2_u∗ i ∂2_x j ui+ 1 − β Re ∂τ∗ ij ∂xj ui− u∗jui ∂Ui ∂xj . (25)

(17)

By adding the two equations and dividing the final equation by 2, we finally arrive at 1 2 ∂(u∗_iui) ∂t = − 1 2(u ∗ iUj ∂ui ∂xj + Uj ∂u∗_i ∂xj ui) − 1 2(u ∗ i ∂p ∂xi +∂p ∗ ∂xi ui)+ β 2Re( ∂2_u i ∂2_x j u∗i + ∂2_u∗ i ∂2_x j ui) + 1 − β 2Re (u ∗ i ∂τij ∂xj +∂τ ∗ ij ∂xj ui) − 1 2(u ∗ iuj ∂Ui ∂xj + u∗jui ∂Ui ∂xj ). (26) Integrating by parts and using the continuity equation, it yields

∂(e) ∂t = ∂ ∂xj −1 2uiu ∗ iUj− 1 2(uip ∗_+u∗ ip)δij+ β 2Re(u ∗ i ∂ui ∂xj +ui ∂u∗_i ∂xj )+1 − β 2Re (u ∗ iτij+uiτij∗) − β Re ∂u∗_i ∂xj ∂ui ∂xj −1 − β 2Re (τij ∂u∗_i ∂xj + τ_ij∗ ∂ui ∂xj ) −1 2(u ∗ iuj ∂Ui ∂xj + uiu∗j ∂Ui ∂xj ), (27)

where e(y, t) is defined as u∗iui

2 , the perturbation energy density. The terms

in the square brackets are transport terms, which only redistribute the energy inside the domain. Thus, if the boundary condition is Dirichlet type, as we assume here for the velocity being zero at the boundaries, the transport terms contribute nothing to the total energy rate. Therefore, integrating the above

equation in the domain Ω = ab with a = 2π_α and b = 2π_γ being the lengths in

x and z directions, we obtain the time variation of energy density in the (x, z) plane 1 Ω Z Ω ∂(e) ∂t dV = 1 ab Z 1 −1 Z a 0 Z b 0 ∂(e) ∂t dzdxdy = Z 1 −1 ∂(e) ∂t dy = − Z 1 −1 β Re ∂u∗ i ∂xj ∂ui ∂xj dy− Z 1 −1 1 − β 2Re (τij ∂u∗ i ∂xj +τij∗ ∂ui ∂xj )dy− Z 1 −1 1 2(u ∗ iuj+uiu∗j) ∂Ui ∂xj dy. (28) The last term in the equation is the production against the base flow. The first and second ones on the right hand side of the last equal sign are the dissipa-tion terms, including the viscous dissipadissipa-tion and the polymer stress dissipadissipa-tion (or polymer work rate), latter of which is the energy exchanged in the interac-tion between the polymer fluctuating stress and the fluctuating flow field. For the ease of later discussion, these different terms are given names according their natures. For example, in viscous dissipation terms (V D), V D11 refers to

−1 Ω R Ω β Re ∂u∗ 1 ∂x1 ∂u1

∂x1dV (the other terms are similarly named). The same rule

ap-plies to polymeric work rate (P D) terms and production against the base flow shear (P ) terms. For the parallel flow cases considered here, production term reduces to −1R 1(u∗u2+ u1u∗)∂UdV .

(18)

Furthermore, according to the wave-like assumption of perturbation, it can be shown that Re= 1 Ω R Ω ∂(e) ∂t dV 1 Ω R ΩedV = R Ω 1 2u˜iu˜ ∗ id(e−iwt+iw ∗_t) /dt dV R Ω 1 2u˜iu˜ ∗ ie−iwt+iw ∗_t dV =d(e −iwt+iw∗t)_/dtR Ω 1 2u˜iu˜ ∗ idV e−iwt+iw∗_tR Ω 1 2u˜iu˜ ∗ idV = −iw + iw∗= 2wi. (29)

wi is the imaginary part of least stable eigenvalue. The reason why the time

exponentials could be extracted from the integrals is that the control volume

is assumed to be time-invariant. Therefore, the value of Re should be exactly

equal to two times the imaginary part of the least stable eigenvalue. This fact bridges the results of linear stability analysis and energy analysis. In fact, energy analysis is performed based on the information of resultant eigenfunction, while the classical linear stability results concern the least stable eigenvalue. These two categories of information come from the solution to the same eigenvalue problem. Therefore, consistent results will always yield Re= 2ωi.

2.6 Non-Modal Analysis

The classical linear stability analysis begins with linearizing the Navier-Stokes equation, which yields O-S equation and Sq equation. The criterion then is to examine whether there are any eigenvalues having positive imaginary parts. However, this criterion fails to predict transition in the cases of Newtonian Poiseuille flow and Couette flow, both of which turn to turbulence at a much lower critical Reynolds number compared to the theoretical prediction of clas-sical linear stability analysis. In those cases, the instability originates from a mechanism called transient growth due to the non-orthogonality of the eigen-functions of O-S and Sq operators in these flows. Specifically, the amplification of disturbance energy in a short time (transient growth) stems from the coupling operator in the O-S and Sq system. The energy during transient growth can be thousand times larger than the initial perturbation energy leading to instability before the perturbation energy asymptotically decays to zero as predicted by the wave-like assumption above.

The linearized equation (22) discussed above could be rewritten as ∂

∂tφ = Lφ, (30)

namely, φ = etLφ0. Therefore, the largest possible energy growth in the global

time domain up to time t is the norm of the linear evolution operator T = etL, i.e., G(t) = max φ0 ||φ||E ||φ0||E = max φ0 ||etL_φ 0||E ||φ0||E = ||etL||E (31)

(19)

The initial condition which gives rise to this maximum growth at time t is called the optimal initial condition (Farrell 1988 [?]; Reddy, Schmid & Henningson 1993 [38]). For the Newtonian Poiseuille flow, the optimal initial condition is found when the streamwise wavenumber is approximately 0. The aim of this part in the thesis is to investigate the effects of polymer additives on the transient growth and the optimal initial conditions of these flows. To compute G, we should correctly formulate the energy. The kinetic energy associated with the system is E(t) = 1 2 Z Ω φ∗M φdV (32)

To ease the calculation, the energy norm defined as above is in L2_-norm.

There-fore, we have G(t) = max φ0 ||φ||E ||φ0||E = max φ0 ||F φ||2 ||F φ0||2 = max φ0 ||F etL_F−1_{F φ} 0||2 ||F φ0||2 = ||F etLF−1||2, (33)

where F is the Cholesky factorization of M = F F∗_{. The last equation holds}

because it could be understood that we search the maximum in the space of L2_{-norm, which means that ||F φ}

0||2is the initial disturbance.

In our analysis, we are interested in how the flow will response if an initial disturbance is added in the fluid velocity. In order to do this, filter matrices are

introduced such that φ = Bφin and φout = Cφ with B and C being

(Klinken-berg, de Lange & Brandt 2011 [25])

B =             1 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0             , C =1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 (34)

Matrix B allows that only the flow disturbance is brought into the system and matrix C only the transient growth of flow field is considered. Therefore, the

evolution operator from φin to φout becomes T = CetLB. In addition, the

input and output energy matrices are accordingly defined as Min= FinFin∗ and

Mout = FoutFout∗ . Thus, the optimal energy growth in L2-norm with given

(20)

G(t) = max φ0 ||φout(t)||Eout ||φin(0)||Ein = max φ0 ||T φin(0)||Eout ||φin(0)||Ein = max φ0 ||FoutT φin(0)||2 ||Finφin(0)||2 = max φ0 ||FoutT Fin−1Finφin(0)||2 ||Finφin(0)||2 = ||FoutT Fin−1||2= ||FoutCetLBFin−1||2. (35)

(21)

3 Numerical Methods

3.1 Chebyshev Points

To solve the equations (22) derived above, we use spectral collocation method based on Chebyshev points, which is suitable for bounded domains. They are defined as (Trefethen 2000 [44])

yj= cos(

(j − 1)π

n − 1 ), j = 1, 2, ..., n, (36)

where n is the number of Chebyshev modes. Approximately, any function f (y) could be interpolated by a n − 1 degree polynomial

f (y) ≈ pn−1(y) = n

X

j=1

Lj(y)f (yj), (37)

where f (yj) is the function value at yj and Lj(y) are defined as

Lj(y) = (−1)j cj 1 − y2 (n − 1)2 T_n−10 (y) y − yj , j = 1, 2, ..., n. (38)

Tn−1(y) is the Chebyshev polynomial of degree n − 1 (so Lj(y) is of degree

n − 1) and c1 = cn = 2 and c2 = ... = cn−1= 1. From another point of view,

this formulation implies that we have n interpolating equations at n Chebyshev points with n unknowns in pn−1(y). Therefore the problem is well-defined.

In terms of boundary conditions, for the velocity, we consider v = ∂

∂yv = 0 at

the boudaries, which amounts to four additional equations put into the problem. Theoretically, in the case of Chebyshev nodes with the endpoints y = ±1 deleted, this will require the interpolant polynomial to be converted to (Canuto et al. 2007 [7]; Weideman & Reddy 2000 [46])

L+_j(y) = (1 − y

2

1 − y2 j

)2Lj(y), (39)

The plus sign+_{denotes the polynomial or derivative with boundary conditions.}

Note that in this way, the order of the interpolating polynomial L+_j(y) is in-creased by four. But in our code, the Dirichlet boundary condition v = 0 is alternatively realized by directly forcing the velocity vectors v to be zero at the two ends, meaning that we do not need to consider these two equations in our interpolation formulation. Thus, the modification of the interpolant polynomial only accounts for the non-penetration boundary conditions _∂y∂ v = 0, that is,

L+_j(y) = 1 − y

2

1 − y2 j

(22)

The problem is still well-defined because the number of the equations is still equal to the number of unknowns. Then the derivatives with boundary condi-tions are obtained by applying chain rules of derivative

D_c+= −2y 1 − y2 j Lj(y) + 1 − y2 1 − y2 j DcLj(y) (41) D+2_c = − 2 1 − y2 j Lj(y) − 4y 1 − y2 j DcLj(y) + 1 − y2 1 − y2 j D2_cLj(y) (42) D_c+3= − 6 1 − y2 j DcLj(y) − 6y 1 − y2 j D_c2Lj(y) + 1 − y2 1 − y2 j D3_cLj(y) (43) D_c+4= − 12 1 − y2 j D2_cLj(y) − 8y 1 − y2 j D_c3Lj(y) + 1 − y2 1 − y2 j D4_cLj(y) (44)

Dc is the derivative matrix for the Chebyshev points. In our code, the first

three derivative matrices D+

c, Dc+2and D+3c are solely used in constructing the

forth derivative matrix D+4

c because only the last one needs the four boundary

condition to be fully imposed.

As mentioned before, we are about to solve an eigenvalue problem A_∂t∂φ = Bφ, where φ is (v, η, c11, c22, c33, c12, c13, c23)T. Matrices A and B are (8N +

12) × (8N + 12) in 3 dimensions. In Appendix B, the explicit expression of matrices A and B will be given. Left-multiplying A matrix in A_∂t∂φ = Bφ and using the assumption of wave-like perturbation, we obtain ω ˜φ(y) = iA−1B ˜φ, where ω is the eigenvalue and ˜φ(y) is the eigenfunction. The parameter space in the thesis is (α, γ, W, Re, β, L).

3.2 Stretching of Grids

In the case of jet and mixing layer flows the most significant derivative of the flow profile is found near the origin while the Chebyshev polynomial clusters the points around the two ends. Therefore we need to stretch the Chebyshev points and gather them around the origin to increase accuracy, as shown in figure (3).

−1 −0.5 0 0.5 1 −2 −1 0 1 2

Figure 3: Comparison of two sets points distribution. The blue is the Chebyshev distribution. The red is the stretched distribution that we want.

(23)

The mapping function between the original Chebyshev points and the stretched ones is a 5-ordered polynomial. For the jet profile u = sech2(y) we use a map-ping function interpolated from the one used in (Govindarajan 2004 [17])

y = 0.965588y5c− 0.162978y 3

c + 0.191655yc, (45)

where yc are the Chebyshev points. The mapping between the two distributions

can be clearly seen in figure (4), where the horizontal coordinates of the points are the Chebyshev points on (-1,1) and the vertical coordinates of the points are the mapped points, which are now gathering near the origin. Other forms of the mapping function exist if the flow profile indicates the most significant derivative at other regions. This could be done by adjusting the shape of the mapping function to place more points on that region.

−1 −0.5 0 0.5 1 −1 −0.5 0 0.5 1

Figure 4: Mapping between two points distribution. The abscissa is Chebyshev distribution of points on [-1,1]. The ordinate is the stretched distribution of points on [-1,1].

With the change of the discretized points, the derivative matrices need also to be modified. For the stress, we assume no boundary condition. Therefore, the derivatives are modified according to the chain rule as following. Dsrepresents

the stretched derivative along the normal direction. The first order derivative in the stretched domain is

Ds=

dyc

dyDc. (46)

The second order derivative in the stretched domain is D2_s= (dyc dy) 2_D2 c + d2_y c dy2Dc. (47)

The third order derivative in the stretched domain is D3s= ( dyc dy ) 3_D3 c + 3 dyc dy d2yc dy2 D 2 c+ d3yc dy3 Dc. (48)

(24)

The forth order derivative in the stretched domain is D4_s= (dyc dy ) 4_D4 c + 6( dyc dy) 2d 2_y c dy2D 3 c + 3( d2yc dy2 c )2D_c2 + 4dyc dy d3_y c dy3D 2 c + d4_y c dy4 c Dc. (49)

Replacing the above Dc with D+c we could get the stretched derivative matrices

for the flow variables with boundary condition.

3.3 Code Verification

The Oldroyd-B code is verified by comparing our results with the existing liter-ature on the same viscoelastic flow problem, while the FENE-P code is verified by setting large L to see whether it returns to Oldroyd-B model.

Table 2: Verification of our Oldroyd-B code

Our code Reference [43]

Critical eigenvalues 0.34089441 + 1.9888×10−7 i 1.9696×10−7 + 0.34089442i

In the table (2), it is shown that our Oldroyd-B code for Poiseuille flow can yield the same critical eigenvalue at Re = 3960, β = 0.5, W = 3.96, α = 1.15, γ = 0 as the other published literatures have (Sureshkumar & Beris 1995 [43]). The two critical eigenvalues have the reversed position of the real and imaginary parts because two different definitions of the wave-like perturbation. The FENE-P model code is expected to generate the same eigenspectrum as the Oldroyd-B model code when the parameter L is very large and other parameters are the same. The following figure shows this fact at Re = 3960, β = 0.5, W = 3.96, α = 1.15, γ = 0. The number of modes for discretization is 101 for velocity and normal vorticity and 103 for the conformation tensor.

0 0.5 1 1.5 −3.5 −3 −2.5 −2 −1.5 −1 −0.5 0

Figure 5: Eigenspectrum for Oldroyd-B model (red) and FENE-P model (green) when L is large at Re = 3960, β = 0.5, W = 3.96, α = 1.15, γ = 0.

(25)

Besides, the results of instability growth rates versus streamwise wavenumber α shown in Azaiez & Homsy (1994 [?]) have been regenerated in figure (6) as part of the code validation. The parameter is G defined as (1−β)W/Re, representing the relative strength of elastic effect and viscous effect. The problem considered in their paper is inviscid free shear flow, so the Re in our setting is 1000.

0 0.2 0.4 0.6 0.8 1 0 0.05 0.1 0.15 0.2 α ω i G=0 G=1 G=4

(26)

4 Results of Jet and Mixing Layer

In this section, we present the results of modal stability analysis for jet and mixing layer of dilute polymer suspensions modeled by Oldroyd-B. Since it was previously shown that the Squire’s theorem is valid also for the Oldroyd-B model (Bistagnino 2007 [3]), only 2D disturbances are considered in this part.

4.1 Neutral Curves and Confinement

To see the influence of elasticity on flow stability, the neutral curves in the (Re, α) plane are plotted according to imaginary part of the least stable eigen-value for different eigen-values of viscoelastic parameters. For the jet, the base flow is u = sech2(y) and the neutral curve in the case of Newtonian flow (viscosity ratio β = 1) is well-known and repeated in figure (7). The black line is the neutral curve, on which the imaginary part of the least stable eigenvalue is zero. The red region represents instability and blue region stability (see the color bar at right). Re α 5 10 15 0.5 1 1.5 −0.2 −0.15 −0.1 −0.05 0 0.05

Figure 7: Neutral curve for Newtonian jet flow.

In the case of viscoelastic flow, the neutral curves change at different W as shown in figure (8). When the Weissenberg number is moderate, the critical Re increases with increasing W , meaning that elasticity stabilizes the flow. How-ever, when the Weissenberg number is large enough, there is a new unstable region appearing at small Re and large α representing a pure elastic instability. When we increases the resolution, the instability at vanishing Re still exists, implying that it is not due to numerical instability. In figure (8) and other alike figures without color, the different regions are denoted by their own

(27)

char-acteristic natures, i.e., HI denotes hydrodynamic instability usually occurring at large Re, EI elastic instability and S stability. (Sometimes, the border be-tween HI and EI may blur. These two mechanisms are easier to distinguish in the energy analysis presented subsequently.) Changing β, i.e., the polymer concentration or viscosity, is another way to change the viscoelastic property of the flow. Decreasing β has a similar effect of increasing W , as shown in figure (9). Besides, the disturbance wavenumber α at which instability occurs in the above two figures is not varying significantly with different W and β.

Re α 5 10 15 0.5 1 1.5 Newtonian W=5, β=0.8 W=7, β=0.8 W=10, β=0.8 HI EI S

(28)

Re α 5 10 15 0.5 1 1.5 Newtonian W=7, β=0.8 W=7, β=0.6 HI EI S

Figure 9: Neutral curve of jet flow at different β.

In addition, different jet confinements are considered with the flow domains being 10, 14, and 20 and the jet velocity profile unaltered. There are two types of Re in the confinement research context: one Re is based on the width of jet or wake, and the other one is based on the width of channel. Our Re is the former one since the characteristic length in Re is where the velocity is 0.41Uc

(therefore the length is simply 1 in jet according to its base flow profile in ta-ble (1)). In figure (10), confinement effect is stronger when flow domain is 10 since walls are nearer to the flow centerline compared to other cases. Former investigation on jet (wake) confinement by other researches could be found in Chen, Pritchard & Tavener (1995 [8]) and Rees & Juniper (2010 [39]). Note that these investigations pertain to Newtonian flow. Our results show that in the viscoelastic jets, for the instability of hydrodynamic nature (HI) stronger confinement implies more stable flow (without developing boundary layer) and these are consistent with the results in [8], where it shows that stronger con-finement leads to larger critical Re whose characteristic length is based on the diameter of the cylinder (table (1) in their paper). Therefore, it can be re-marked that even with the presence of polymer, the confinement effect on the hydrodynamic instability (high Re) remains the same as in the Newtonian case. Interested readers are referred to the references listed here and therein for more useful discussion.

(29)

Re α 5 10 15 0.5 1 1.5 Length=20 Length=14 Length=10 HI S EI

Figure 10: Neutral curves for different confinements in jet flow at W = 10, β = 0.8, γ = 0.

Next we consider a mixing layer whose base flow is u = tanh(y); the neutral curve for Newtonian fluid is plotted in figure (11). Unlike the jet flow, the mixing layer does not show a new unstable region at vanishing Re when a stronger elastic effect is present (see figure (12)). Instead, the critical Re first slightly increases at W ≈ 1 and then decreases when W increases. In fact, as will be shown in the energy analysis, the instability at large W will turn to elastic nature (instability due to the destabilizing effect of the increasing polymer normal stress [40]) when E = W/Re & 1. Decreasing β has a destabilizing effect on the flow, as shown in figure (13). The energy analysis will be performed subsequently to reveal the reason.

(30)

Re α 5 10 15 0.2 0.4 0.6 0.8 1 1.2 −0.1 −0.05 0 0.05 0.1

Figure 11: Neutral curve for the Newtonian mixing layer flow.

Re α 5 10 15 0.2 0.4 0.6 0.8 1 1.2 Newtonian W=5 W=7 W=10 W=1 S HI

(31)

Re α 5 10 15 0.2 0.4 0.6 0.8 1 1.2 Newtonian β=0.8 β=0.6 S HI

Figure 13: Neutral curve of mixing layer flow at different β at W = 7. Before next section, we give the reason why in the mixing layer flow there is no elastic instability region emerging at vanishing Re while in jet flow this oc-curs. The main argument is that in mixing layer, the shear effect at very low Re only extends the polymer molecules at the flow interface of different velocities while in the jet, the polymer chains are stretched in the two adjacent layers, interaction between which will probably give rise to instability. To verify this argument, two new base flows are considered as

u3= tanh(y)sech2(y) y ∈ (−10, 10) (50)

u4= −2sech2(y)tanh2(y) + sech4(y) y ∈ (−10, 10) (51)

−10 −5 0 5 10 −0.4 −0.2 0 0.2 0.4 y u (a) −1 −0.5 0 0.5 1 −0.5 0 0.5 1 1.5 y du/dy (b)

(32)

−10 −5 0 5 10 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 y u u₄ Jet flow (a) −1 −0.5 0 0.5 1 −3 −2 −1 0 1 2 3 y du/dy (b)

Figure 15: ew base flow (a) u4 and jet, (b) du_dy4.

The flow u3 consists of three inflectional points and u4 four such points,

see figure(14b) and (15b). One can see that in u3 profile there is an interface

between flows of two different velocities at the centerline, similar to mixing layer, however, we expect that disturbance at very low Re will lead this flow to instability because there are adjacent shearing layers besides the interface. This is true since at β = 0.8, W = 20, R = 0.1, α = 0.4, γ = 0, the flow

becomes elastically unstable (imaginary part of least stable eigenvalue ωi =

0.003498934598863). On the other hand, u4 is compared to jet flow in figure

(15a). Because the adjacent shearing layers are closer in u4 and there are also

small adjacent shearing layers at negative velocities, we expect that the growth rate of u4is larger than jet flow at the same parameters when there is an elastic

instability region in the vanishing Re limit. This is true according to table (3).

Table 3: Disturbance growth rate ωi of u4and jet at the same parameters

F lows Growth rate

jet 0.017235583628661

u4 0.072818116187888

4.2 Eigenfunction

According to linear stability analysis, eigenfunction of the system can be readily obtained from the eigenvalue problem (22). The eigenfunction for jet flow at W = 10, β = 0.8, Re = 14.77, α = 0.82, γ = 0 is shown in figure (16). From neutral curve (8), we know that this instability is of hydrodynamic nature. The eigenfunction at W = 10, β = 0.8, Re = 0.01, α = 0.82, γ = 0 is shown in figure (17). In this case, the instability is due to the elasticity of the polymer molecules. One can note that the hydrodynamic jet mode is antisymmetric in u, whereas the elastic mode is symmetric in u. The eigenfunction for mixing layer at W = 10, β = 0.8, Re = 2, α = 0.66 is shown in figure (18). The instability is of elastic nature, thus, the u mode is also symmetric.

(33)

−100 −5 0 5 10 0.5 1 1.5x 10 −4 y u mode v mode (a) −100 −5 0 5 10 0.01 0.02 0.03 0.04 y τ 11 τ 22 τ 12 (b)

Figure 16: (a) Eigenfunction and (b) polymer stress of jet flow at W = 10, β = 0.8, Re = 14.77, α = 0.82. It is hydrodynamic instability. −100 −5 0 5 10 1 x 10−4 y u mode v mode (a) −100 −5 0 5 10 0.01 0.02 0.03 y τ 11 τ 22 τ 12 (b)

Figure 17: (a) Eigenfunction and (b) polymer stress of jet flow at W = 10, β = 0.8, Re = 0.01, α = 0.82. It is elastic instability. −100 −5 0 5 10 1 x 10−4 y u mode v mode (a) −100 −5 0 5 10 0.01 0.02 0.03 y τ 11 τ 22 τ 12 (b)

Figure 18: (a) Eigenfunction and (b) polymer stress of mixing layer flow at W = 10, β = 0.8, Re = 2, α = 0.2. It is elastic instability.

(34)

4.3 Energy Analysis

The energy analysis can be used to explain instability mechanisms in the flows. In order to make the comparison between different cases meaningful, all the terms in equation (28) are normalized with the integral of the perturbation

energy R

ΩedV , so we are dealing with the normalized energy density as in

equation (29). Energy analysis is performed with different Re and W , both of which can provide the quantitative differences between various energy terms with varying elastic effect in the flows.

Table 4: Energy analysis of Newtonian jet flow

Terms Re = 4(N ) Re=6(N) Re = 10(N ) Re = 16(N ) V D11 -0.040 -0.036 -0.030 -0.022 V D12 -0.131 -0.148 -0.148 -0.131 V D21 -0.140 -0.084 -0.042 -0.023 V D22 -0.040 -0.036 -0.030 -0.022 V D -0.350 -0.304 -0.249 -0.198 P 0.231 0.287 0.332 0.348 Total -0.119 -0.017 0.082 0.150

Table (4) shows the results of Newtonian jet flow with varying Re (The N behind Re represents Newtonian flow). Production against mean flow increases with Re whereas viscous dissipation decreases leading to instability. This is the well-known mechanism leading to instability in Newtonian flow. The result of energy analysis on viscoelastic flow at W = 10, β = 0.8, α = 0.6, γ = 0 (Chebyshev modes are 131) and different Re is shown in table (5). P D, as mentioned before, refers to elastic work rate, V D viscous dissipation (always negative) and P production against the base flow shear. The subscript numbers represent different directions.

(35)

Table 5: Energy analysis of viscoelastic jet flow with varying Re Terms Re = 4 Re = 6 Re = 10 Re = 16 P D11 1.041 0.859 0.693 0.182 P D12 -0.252 -0.227 -0.211 -0.095 P D21 -0.045 -0.034 -0.023 -0.028 P D22 -0.034 -0.024 -0.015 -0.019 V D11 -0.089 -0.062 -0.038 -0.021 V D12 -0.359 -0.299 -0.238 -0.117 V D21 -0.055 -0.034 -0.019 -0.015 V D22 -0.089 -0.062 -0.038 -0.021 P D 0.711 0.574 0.444 0.040 V D -0.591 -0.456 -0.334 -0.174 P -0.096 -0.107 -0.119 0.217 T otal 0.023 0.011 -0.009 0.082 2ωi 0.023 0.011 -0.009 0.082

The prediction is consistent with the result of linear stability analysis, i.e., time variation of total energy is positive if the flow is unstable while negative if the flow is stable. Besides, as mentioned before, the value of Rein equation

(29) should be exactly equal to two times the imaginary part of least stable eigenvalue. Figure (19) shows this fact, so that we validate our results. The relative difference between these two values is usually 10−7at Chebyshev modes being 131. 0 5 10 15 0 0.02 0.04 0.06 0.08 0.1 Re Eigenvalues*2 R e

Figure 19: The least stable eigenvalues and Reat W = 10, β = 0.8, α = 0.6, γ =

(36)

In table (5), the instability at Re = 16 is of hydrodynamic nature, while at Re = 6 we have elastic instability. Note that in hydrodynamic instability the polymer work rate P D may be positive (in this particular case W = 10, so W/Re ≈ 1, the elastic effect is still significant), but the main destabilizing effect comes from the production P ; the opposite holds for EI. With increasing Re, the production P increases, viscous dissipation V D decreases and the elastic work rate P D decreases. Interestingly, it can be seen that production against the mean shear is positive in the hydrodynamic instability region (Re = 16) while is negative in the elastic instability region (Re = 6 or 4). It is important to note that when an instability occurs, disturbance energy is generated only by the term P D11, interaction between the fluctuating polymer stress and the

fluctuating flow field, i.e., −R

Ω 1−β 2Re(τ11 ∂u∗₁ ∂x1 + τ ∗ 11 ∂u1

∂x1)dV in the energy equation

(28). Besides, comparing the results of viscoelastic and Newtonian flows, we can see that with the presence of polymers, the production against base flow shear is reduced. The data are also presented visually in figure (20).

4 6 8 10 12 14 16 −0.4 −0.2 0 0.2 0.4 0.6 Re Energy PD VD P T VD(N) P(N) T(N)

Figure 20: PD, VD, P and T in viscoelastic and Nowtonian jet flow with varying Re.

Next we present the energy analysis on jet with varying W at β = 0.8, Re = 7, α = 0.6, γ = 0 in table (6). W = 0 is the Newtonian limit. The data is as well visualized in figure (21a). One can see that at small W , the sum of P D and V D terms remains constant. When W increases to 5, the P D and P curves cross over each other, after which P D increases significantly, leading the flow to elastic instability. At the same time, production against mean shear P decreases dramatically with W around 5, but when P D reaches a saturation, P as well decreases slowly. Therefore, it is clear that at small W , the flow is stabilized (see black line T ) due to the polymer molecules reducing production against mean shear P , and destabilized at larger W because of the increasing

(37)

P D, especially P D11. The crossover point of P D and P curves at Re = 7 occurs

earlier than Re = 10 at fixed β. This indicates that at larger Re, the elastic effect of polymer should be more pronounced to destabilize the flow. The number (= W/Re), defined as characteristic elastic relaxation time over viscous diffusion time, captures this fact. One can see that at β = 0.8, the point where the black line T reaches its minimum and starts to increases is around ≈ O(1).

Table 6: Energy analysis of viscoelastic jet flow with varying W

Terms W = 0 W = 1 W = 5 W = 7 W = 10 P D11 0 -0.0092 0.2098 0.6050 0.8016 P D12 0 -0.0383 -0.0988 -0.1453 -0.2205 P D21 0 -0.0153 -0.0424 -0.0245 -0.0300 P D22 0 -0.0192 -0.0412 -0.0219 -0.0206 V D11 -0.0346 -0.0266 -0.0376 -0.0527 -0.0534 V D12 -0.1502 -0.1072 -0.1380 -0.2372 -0.2784 V D21 -0.0683 -0.0557 -0.0447 -0.0296 -0.0288 V D22 -0.0346 -0.0266 -0.0376 -0.0527 -0.0534 P D 0 -0.0819 0.0274 0.4133 0.5306 V D -0.2877 -0.2161 -0.2579 -0.3722 -0.4142 P 0.3039 0.2943 0.1691 -0.1002 -0.1111 T otal 0.0163 -0.0036 -0.0614 -0.0591 0.0053 2ωi 0.0163 -0.0036 -0.0614 -0.0591 0.0053 0 5 10 15 −1 −0.5 0 0.5 1 W Energy PD VD P T (a) 0 5 10 15 −0.4 −0.2 0 0.2 0.4 0.6 W Energy PD VD P T (b)

Figure 21: PD, VD, P and T in viscoelastic jet flow with varying W at β = 0.8, α = 0.6, γ = 0 (a) Re = 7, (b) Re = 10.

Energy analysis on the mixing layer has also been performed for varying W , see table (7) and figure (22). The result shows the same trend as we find for P D and V D in jet with small elastic effect, but the production P seems constant in the presence of polymer molecules, while at large W , though the destabilizing effect of polymer is significant, the viscous dissipation also increases by the same order of amount, therefore the instability is still brought about by the production against mean shear.

(38)

Table 7: Energy analysis of viscoelastic mixing layer at β = 0.8, Re = 2, α = 0.6, γ = 0 Terms W = 0 W = 0.5 W = 1 W = 2 W = 5 P D11 0 -0.0087 -0.0075 0.0033 0.1400 P D12 0 -0.0249 -0.0272 -0.0302 -0.0683 P D21 0 -0.0004 -0.0019 -0.0067 -0.0166 P D22 0 -0.0092 -0.0097 -0.0106 -0.0127 V D11 -0.0212 -0.0169 -0.0168 -0.0171 -0.0198 V D12 -0.1297 -0.0986 -0.0950 -0.0957 -0.1430 V D21 -0.0188 -0.0151 -0.0152 -0.0149 -0.0122 V D22 -0.0212 -0.0169 -0.0168 -0.0171 -0.0198 P D 0 -0.0431 -0.0464 -0.0442 0.0424 V D -0.1910 -0.1474 -0.1438 -0.1448 -0.1950 P 0.2084 0.2048 0.2026 0.2040 0.2078 T otal 0.0174 0.0143 0.0123 0.0149 0.0554 2ωi 0.0174 0.0143 0.0123 0.0149 0.0554 0 5 10 15 −0.4 −0.2 0 0.2 0.4 W Energy PD_VD P T

Figure 22: PD, VD, P and T in viscoelastic mixing layer with varying W at β = 0.8, Re = 2, α = 0.6, γ = 0.

To summarize, in the jet flow, but not for mixing layer, polymer stretching creates production of disturbance kinetic energy and an instability of elastic na-ture at low Re, whereas mixing layer only stretches the polymer at the interface. At larger Re and small W for both flows, the instability is of hydrodynamic na-ture, caused by production against the mean shear, and polymers have a slightly stabilizing effect at small elastic effect in the jet flow. At larger W in jet, the

(39)

main destabilizing effect is brought about by the polymer normal stress, leading to instability, but in mixing layer, the production against mean shear is still the main reason leading to instability and the polymer destabilizing effect is dissipated by the viscosity.

(40)

5 Results of Poiseuille Flow

The viscoelastic Poiseuille flow is investigated and we present here the results for neutral curves, modal energy analysis and non-modal analysis.

5.1 Neutral Curves

Neutral curves in this chapter are plotted in (α,Re) and (W, Re) planes with different viscoelastic parameters. The results generally depict FENE-P model in 2D, but results for Newtonian and Oldroyd-B fluids in flows at the same pressure gradient will be presented as well . (Note that the Squire’s theorem for FENE-P Poiseuille flow has not been proven.) There are no background color in the figures indicating stable and unstable regions, but the convention is that on the left hand side of curve (low Re), the flow is stable and on the right hand side unstable. Re α 5000 5500 6000 6500 7000 7500 8000 8500 0.8 0.85 0.9 0.95 1 1.05 1.1 1.15 W=10 W=5 W=1 W=0.5 Newtonian

Figure 23: Neutral curves at β = 0.9, L = 60, γ = 0 with different W . Figure (23) displays the effect of W at β = 0.9, L = 60, γ = 0. It can be seen that at larger W , the strong elastic effect results in stabilization for FENE-P fluids . But interestingly, at small W (around 1), the flow is less stable. This trend may be clearer in the following figure (26).

(41)

5000 5500 6000 6500 7000 7500 8000 8500 0.8 0.85 0.9 0.95 1 1.05 1.1 1.15 Re α β=0.7 β=0.8 β=0.9 Newtonian (a) 5000 5500 6000 6500 7000 7500 8000 8500 0.8 0.9 1 1.1 β=0.7 β=0.8 β=0.9 Newtonian (b)

Figure 24: Neutral curves in (α, Re) space at W = 10, L = 60, γ = 0 with different β (a) for FENE-P, (b) for Oldroyd-B.

Figure (24) shows the effect of viscosity ratio β at W = 10, L = 60, γ = 0 for FENE-P and W = 10, γ = 0 for Oldroyd-B. One can see that with decreasing β (which amounts to increasing polymer molecule concentration), the critical Re for FENE-P fluid decreases, i.e., at this W , stronger elastic effect stabilizes the flow. For Oldroyd-B fluid at W = 10, decreasing β (with β in Newtonian case being 1) first stabilizes, then destabilizes the flow. But the effect is smaller than the stabilization of FENE-P model once one compares the figures above. Note that the conclusion is drawn within the range of large β.

20000 4000 6000 8000 10000 12000 14000 5 10 15 Re W β=0.5 FENE−P β=0.7 FENE−P β=0.9 FENE−P β=0.5 Oldroyd−B β=0.7 Oldroyd−B β=0.9 Oldroyd−B

(42)

−50000 0 5000 10000 15000 5 10 15 (Re−5772)/(1−β) W β=0.5 FENE−P β=0.7 FENE−P β=0.9 FENE−P β=0.5 Oldroyd−B β=0.7 Oldroyd−B β=0.9 Oldroyd−B

Figure 26: Neutral curves in ((Re-5772)/(1-β),W) plane at L = 60, γ = 0 with different β. −50000 0 5000 10000 15000 20000 1 2 3 4 (Re−5772)/(1−β) W* ω r _{β=0.4 FENE−P} β=0.5 FENE−P β=0.6 FENE−P β=0.7 FENE−P β=0.8 FENE−P β=0.9 FENE−P

Figure 27: Neutral curves in ((Re-5772)/(1-β),W*ωr) plane at L = 60, γ = 0

with different β.

Figure (25) and (26) show the relation between W and the critical Re for different β for the two models, FENE-P and Oldroyd-B. Critical Re implies that the values of Re are obtained by searching a minimum in the (α, Re) plane for

(43)

each case, so that α might be different at each point. It is notable that at certain β and Re, with increasing W , the flow is first stable, then unstable (W around 2), and finally stable. In figure (25), we can see that at W = 10, the flow modeled by FENE-P is stabilized with decreasing β, the same conclusion deduced from

figure (24). Interestingly, at lower W , the different β curves In figure (25)

collapse into a single one once we scale the abscissa to be (Re − 5772)/(1 − β) In figure (26). This implies that the small elastic effect, leading to the moderate deviation of critical Re from the Newtonian limit (Re = 5772), is not influenced by the concentration of polymer molecules, depicted by 1 − β. The decrease of critical Re at small W may be due to the increase of effective viscosity. Figure (27) shows the result of W × ωrversus critical Re. ωris the real part of the least

stable eigenvalue (the frequency of the disturbance) at the critical condition,

so W × ωr is the polymer relaxation time over the period of the disturbance.

We can see that when the value of W × ωr is less than around 1, the flow is

destabilized with increasing W , and larger than around 1, the flow is stabilized with increasing W . 5000 5500 6000 6500 7000 7500 8000 8500 0.8 0.85 0.9 0.95 1 1.05 1.1 1.15 Re α _L=20 L=60 L=120 L=240 Newtonian Oldroyd−B

(44)

50000 5500 6000 6500 7000 7500 5 10 15 Re W L=60 L=120 L=240 Oldroyd−B

Figure 29: Neutral curves at β = 0.9, γ = 0 with different L.

The effect of L, the maximum extension of polymer molecule, is presented next. The result is shown in figure (28) where we choose W = 10, β = 0.9. Oldroyd-B model is simply the limit of FENE-P model when L goes to infinity, therefore, the result of the Oldroyd-B model is shown as a limiting case of a series of FENE-P model with increasing L. The shorter L, the stronger elastic effect the polymer molecules can display, though the stronger normal stress gradient is found in longer L case. As seen from the figure (28), smaller L shows stabilization of the flow, which in some extent is consistent with the previous results for W = 10 case (stronger elasticity stabilizes the flow). As will be shown in the energy analysis, the reason of smaller L resulting in more stable flow is that at small L, the production against the mean shear P is inhibited. The effect of L is also presented in the (W, Re) plane in figure (29), where Re is the critical Re. But unlike for β in figure (26), at all the W investigated, increasing L destabilizes the flow. This suggests that L, elasticity, always stabilizes the flow, while polymer viscosity simply enhances the effect of the viscoelastic additives. Since there is no Squire’s theorem for viscoelastic flow of non-linear consti-tute model (e.g. FENE-P in our study), we investigated (α, γ) plane as shown in figure (30). The most unstable region is around α = 1, γ = 0. At γ = 0, if γ increases a little, the least stable ωi will decrease compared to γ = 0 case (not

clear from the figure), meaning that the flow in 2D is more unstable in these cases.

(45)

γ α 0.5 1 1.5 0.5 1 1.5 −12 −10 −8 −6 −4 −2 0 x 10−3

Figure 30: Neutral curve in (α, γ) plane at β = 0.9, W = 5, L = 60, Re = 6000.

5.2 Energy Analysis

The energy analysis of FENE-P Poiseuille flow is presented in this section. As we linearize the equation, the expressions for perturbed polymer stress τij are

different in the two constitutive models due to their different definitions. In

Oldroyd-B model, τij is simply

cij

W, where the perturbed conformation tensor

could be readily extracted from the resultant eigenfunction. However, in FENE-P model, τij is

f0_C ij+f cij

W , where f and f

0 _{are defined previously in equations}

Table 8: Energy analysis of viscoelastic Poiseuille flow

Terms W = 0 W = 0.5 W = 2.5 W = 6

P D11 0 -1.378E-04 -2.817E-04 3.416E-04

P D12 0 -1.355E-03 -1.226E-03 -1.641E-03

P D21 0 1.537E-05 2.430E-06 -2.601E-05

P D22 0 -5.426E-05 -6.130E-05 -6.688E-05

V D11 -2.673E-04 -2.406E-04 -2.406E-04 -2.408E-04

V D12 -1.429E-02 -1.288E-02 -1.300E-02 -1.317E-02

V D21 -1.043E-04 -9.385E-05 -9.380E-05 -9.363E-05

V D22 -2.673E-04 -2.406E-04 -2.406E-04 -2.408E-04

P D 0 -1.531E-03 -1.567E-03 -1.392E-03

V D -1.493E-02 -1.345E-02 -1.358E-02 -1.375E-02

P 1.432E-02 1.478E-02 1.561E-02 1.400E-02

Total -6.031E-04 -2.067E-04 4.657E-04 -1.137E-03

Linear Stability Analysis of Viscoelastic Flows

May 9, 2012