Analysis and numerical simulation of the flow of viscoelastic fluids

OF VISCOELASTIC FLUIDS

THE FLOW

OF VISCOELASTIC FLUIDS

Proefschrift

ter verkrijging van de graad van doctor

aan de Technische Universiteit Delft,

op gezag van de Rector Magnificus,

prof.drs. P.A. Schenck,

in het openbaar te verdedigen

ten overstaan van een commissie

aangewezen dóór het College van Dekanen

op 8 november 1988 om 14.00 uur

door

MARTINUS ANTONIUS HULSEIM

geboren te Sint-Oedenrode

werktuigbouwkundig ingenieur

Delft University Press/1988

TR diss

behorende bij het proefschrift

ANALYSIS AND NUMERICAL SIMULATION OF THE FLOW OF VISCOELASTIC FLUIDS

1. Bij gebruik van een niet-Newtons vloeistofmodel moet men zich goed op de hoogte stellen van de veronderstellingen die aan de afleiding ten grondslag liggen. Indien men de daaruit volgende geldigheidsgrenzen overschrijdt en toch een goede overeenkomst vindt tussen experiment en simulatie, berust dit op toeval.

2. De bilineaire termen in visco-elastische modellen van het difFerentiaaltype hebben alleen tot doel het materiaalgedrag onafhankelijk van de waarnemer te maken. Vanuit dit oogpunt gezien, bestaat er geen enkele voorkeur voor de precieze vorm van deze termen en is de grote aandacht die hieraan in de literatuur wordt besteed sterk overdreven en verwarrend voor de nieuweling in het vak reologie.

3. Het resultaat van Hoofdstuk 3 van dit proefschrift dat aan de instroomrand de defor­ matiegeschiedenis van een visco-elastische vloeistof van het Maxwell-type niet volledig kan worden voorgeschreven, is vanuit een fysisch oogpunt gezien merkwaardig en zal door reologen niet snel worden aanvaard.

4. De veronderstelling dat een vloeistof onsamendrukbaar is maakt de druk p tot een uit­ sluitend mathematische grootheid, die bovendien nog van het gekozen vloeistofmodel afhangt en daarom ook niet direkt meetbaar is.

5. De door Evans and Walters experimenteel gevonden grote invloed van kleine veran­ deringen in de geometrie van een scherpe hoek op het globale stromingsbeeld in een kontraktiestroming, zou de niet erg bemoedigende theorie kunnen ondersteunen, dat het zgn. 'high-Weissenberg-number problem' nooit volledig bevredigend kan worden opgelost.

R.E. Evans and K. Walters, "Flow characteristics associated with abrupt changes in geometry in the case of highly elastic liquids", J. of Non-Newt. Fluid Mech., 20(1986)11-29.

6. Het algoritme gebaseerd op stroomlijnintegratie, zoals beschreven in Hoofdstuk 4 van dit proefschrift, is uitermate geschikt voor implementatie op een computersysteem bestaande uit zeer vele onafhankelijk werkende parallele processoren.

7. Mechanica en stromingsleer hebben meer raakvlakken dan men op grond van de huidige samenwerking zou kunnen konkluderen.

8. Diegenen die de kunsthoofd-opnametechniek als de meest ideale beschouwen, vergeten dat een belangrijk mechanisme voor het richtingshoren van de mens bestaat uit de interpretatie van drukveranderingen in de oren t.g.v. kleine willekeurige bewegingen van het hoofd.

en onderzoek voor ogen heeft.

10. Aangezien er klachten zijn over het slechte taalgebruik van afgestudeerde ingenieurs en zowel voor wetenschapsbeoefening als de voor ingenieurspraktijk geldt dat de presentatie net zo belangrijk is als de prestatie, is het noodzakelijk om meer aandacht té besteden aan het correct gebruik van de Nederlandse taal in afstudeerverslagen. 11. Het is aannemelijk dat een verband bestaat tussen de ontkerkelijking van de samen­

leving en de toegenomen raadpleging van psychiaters en psychologen.

12. Het steeds grotere gebruik van PC's door wetenschappers heeft als onbedoeld neven-effekt dat typewerk te duur betaald wordt.

Hulsen, Martinus Antonius

Analysis and numerical simulation of the flow of viscoelastic fluids / Martinus Antonius Hulsen. - Delft:

Delftse Universitaire Pers. - III.

Proefschrift Delft. - Met lit. opg. - Met samenvatting in het Nederlands. ISBN 90-6275-494-5

SISO 533.3 UDC 532(043.3) NUGI 841 Trefw.: hydrodynamica.

Copyright © 1988 by M.A. Hulsen. All rights reserved.

No part of the material protected by this copyright notice may be reproduced or utilized in any form or by any means, electronic or mechanical, including photocopying, recording or by any

Information storage and retrieval system, without written permission from the publisher: Delft University Press.

page

Summary.

1. INTRODUCTION. 1

2. GOVERNING EQUATIONS OF THE FLOW OF VISCOELASTIC FLUIDS. 5

2.1 Basic equations from continuum mechanics. 5 2.2 Constitutive equations for viscoelastic fluids. 7 2.3 Plane steady linear flows and viscometric functions. 14

2.4 Remarks on the system of equations. 22

2.5 Concluding remarks. 25

3. ANALYSIS OF THE SYSTEM OF EQUATIONS. 26

3.1 Introduction. 26 3.2 Convection equation. 27 3.3 Type of the system of equations. 30

3.3.1 Maxwell-type models. 31 3.3.2 Jeffreys-type models. 39 3.4 Boundary conditions. 43 3.4.1 Maxwell-type models. 44 3.4.2 Jeffreys-type models. 52 3.5 Discontinuities. 54 3.5.1 Maxwell-type models. 54 3.5.2 Jeffreys-type models. 60 3.6 Comparison of Maxwell- and Jeffreys-type models. 62

4. A NUMERICAL METHOD FOR COMPUTING STEADY PLANE VISCOELASTIC FLOW. 64

4.1 Introduction. 64 4.2 Discrete form of the system of equations. 65

4.2.1 Decoupling of the system of equations. 65 4.2.2 Momentum and continuity equation. 66

4.2.3 Constitutive equation. 70 4.3 Iterative scheme for solving the equations. 79

4.4 Some results on condition, stability and accuracy. 86 4.4.1 Sensitivity to small perturbations of the

4.4.3 Stability of integration. 4. 4. 4 Accuracy. 4.5 Concluding remarks. 95 98 100 106

5 . COMPUTATION OF THE FLOW OVER A SQUARE CAVITY. 5.1 Introduction.

5.2 Problem definition.

5.3 Computation of the fully developed flow. 5.4 Results of computations.

5.4.1 Discretization. 5.4.2 The Leonov model. 5.4.3 The Giesekus model. 5.4.4 Convergence behaviour. 5.5 Concluding remarks. 107 107 107 109 115 115 116 122 122 131 6 . CONCLUDING REMARKS. L i s t of symbols. R e f e r e n c e s . 132 134 138

APPENDIX A: Some properties of the Leonov and Glesekus models. 143

APPENDIX B: Analytical solutions of the constitutive equations in plane flow.

148

APPENDIX C: Systems of quasi-linear first-order partial differential equations.

155

APPENDIX D: S t r e a m l i n e o r i e n t e d base v e c t o r s .

APPENDIX E: Generalized s o l u t i o n s and d i s c o n t i n u i t i e s .

Samenvatting.

159

162

167

This t h e s i s gives an analysis of the equations for v i s c o e l a s t i c flow.

Much a t t e n t i o n i s paid to the mathematical p r o p e r t i e s of v a r i o u s

material models and how these properties affect the choice of numerical

methods. A numerical method for solving t h e s e e q u a t i o n s i s described

and a p p l i e d t o the computation of v i s c o e l a s t i c flow over a square

cavity.

The governing equations from continuum mechanics are reviewed. Examples

of viscoelastic material models in d i f f e r e n t i a l form are given : the

classical quasi-linear models and the internal deformation models. Some

properties of the internal deformation models are d e r i v e d . The system

of equations i s put i n t o a dimensionless form and the dimensionless

numbers characterizing the flow are discussed.

With r e s p e c t t o the d e r i v a t i o n of the mathematical properties of the

system of e q u a t i o n s , the v i s c o e l a s t i c models are divided i n t o two

groups: the Maxwell-type and the J e f f r e y s - t y p e models. The l a t t e r

models differ from the former by the addition of an extra viscous term.

The mathematical p r o p e r t i e s t h a t are investigated are : type of the

system of p a r t i a l d i f f e r e n t i a l e q u a t i o n s ( c l a s s i f i c a t i o n ) , c o r r e c t

boundary conditions and possible d i s c o n t i n u i t i e s . I t appears that with

respect to these properties the Jeffreys-type models are simple exten­

sions of the Navier-Stokes equations. The Maxwell-models behave in a

more complicated way, for example in plane steady flow the shear s t r e s s

t a n g e n t i a l to the streamlines s a t i s f i e s an e l l i p t i c equation, whereas

the normal stresses are governed by hyperbolic equations.

A numerical method for plane flow based on the finite-element method

and streamline integration i s presented. The numerical method i s par­

t i c u l a r l y s u i t e d for Jeffreys-type models. Some i t e r a t i v e schemes for

solving the equations are d i s c u s s e d . The numerical s t a b i l i t y of the

i n t e g r a t i o n of the constitutive equations i s studied. I t appears that

the internal deformation models have a much b e t t e r s t a b i l i t y than the

discussed.

The numerical method is used to compute the flow over a square cavity.

The convergence behaviour of the i t e r a t i v e scheme i s s t u d i e d . The

r e s u l t s show t h a t i t i s p o s s i b l e to compute viscoelastic flows with

moderately high Deborah numbers, but t h a t the convergence of t h e

i t e r a t i v e scheme i s s t i l l very slow. Hence, better i t e r a t i v e schemes

for the solution of the equations will have to be developed in order to

avoid excessive computing times.

1. INTRODUCTION.

For many f l u i d s t h e c o n s t i t u t i v e e q u a t i o n , i . e . t h e r e l a t i o n between s t r e s s and d e f o r m a t i o n , i s more c o m p l i c a t e d t h a n t h a t g i v e n by t h e Newtonian m o d e l . Examples of such non-Newtonian f l u i d s a r e : molten polymers, polymer s o l u t i o n s , b l o o d , many f l u i d s p r o c e s s e d by t h e food i n d u s t r y , p a i n t s . Some t y p i c a l non-Newtonian phenomena a r e f o r example:

- S h e a r - r a t e - d e p e n d e n t v i s c o s i t y . Most f l u i d s show a s h e a r - t h i n n i n g b e h a v i o u r : t h e s h e a r v i s c o s i t y d e c r e a s e s w i t h i n c r e a s i n g shear -3 -4 r a t e , sometimes d e c r e a s i n g t o 10 or 10 times t h e z e r o - s h e a r - r a t e v i s c o s i t y (see F i g u r e 1 . 1 ) . • v e l o c i t y p r o f i l e v i s c o s i t y o = ' x y n o r m a l —stress d i f f e r e n c e N = <rx)( - <r„y

Figure 1.1 Typical s t r e s s behaviour of non-Newtonian f l u i d s in shear flow.

Non-zero normal-stress differences in steady shear flows (see Figure 1 . 1 ) . Well-known i s the Weissenberg e f f e c t ( s e e Figure 1.2) : a f l u i d climbing up a r o t a t i n g rod. This e f f e c t i s g e n e r a l l y a t ­ t r i b u t e d t o the d i f f e r e n c e between the r a d i a l and t h e angular normal-stress components.

N e w t o n i a n

non —Newtonian

Figure 1.2 The Weissenberg effect.

Increase of the diameter of an extrudate leaving a d i e , s o - c a l l e d

"die swell" (see Figure 1.3).

<=>-Figure 1.3 Die swell.

S o l i d - l i k e behaviour of polymeric liquids subjected to high-speed

impact t e s t s . For example, a ball of " s i l l y putty" will bounce on a

r i g i d t a b l e l i k e an e l a s t i c rubber b a l l . However, if the ball is

placed on a surface i t will s t a r t to flow and slowly spread i t s e l f

on the s u r f a c e due to g r a v i t y forces, similar to a highly viscous

f l u i d .

One of the main goals of the d i s c i p l i n e of rheology i s t o define

c o n s t i t u t i v e e q u a t i o n s (models) of s u f f i c i e n t a c c u r a c y f o r

non-Newtonian f l u i d s . Many non-non-Newtonian fluids are v i s c o e l a s t i c . They ex­

h i b i t p a r t i a l e l a s t i c recovery a f t e r deformation ( f l u i d memory). A

g r e a t number of models are described in the l i t e r a t u r e , ranging from

the classical and rather simple q u a s i - l i n e a r Maxwell models to very

elaborate integral models. Most models describe viscometric steady flow

well, but have often l i m i t e d p r e d i c t i o n c a p a b i l i t y in more complex

transient flows.

Considering the increasing use of synthetic i n d u s t r i a l products, i t i s

not s u r p r i s i n g that after successful numerical simulation of Newtonian

flows, numerical methods for s o l v i n g v i s c o e l a s t i c flow problems are

being developed as w e l l . S t a r t i n g from 1975, many papers appeared on

the subject of numerical simulation of v i s c o e l a s t i c flows in complex

g e o m e t r i e s . A good survey i s given by Crochet, Davies and Walters

(1981). There has been a strong bias towards finite-element methods, in

p a r t i c u l a r the s o - c a l l e d mixed methods, pioneered by Kawahara and

Takeuchi (1977) and further refined by Crochet and h i s coworkers. The

model t h a t has been used in the majority of papers i s the r e l a t i v e l y

simple upper-convected Maxwell model.

The a l g o r i t h m s developed for v i s c o e l a s t i c flow problems often suffer

from severe convergence problems for higher degrees of e l a s t i c i t y in

the flow ( i n d i c a t e d by a dimensionless number: Deborah or Weissenberg

number). However, in the recent papers of Shen (198H), which i s based

on the work by Upadhyay, Upadhyay and Isayev (1986) and Luo and Tanner

(1986b), high Deborah numbers a r e r e p o r t e d by u s i n g r e a l i s t i c

c o n s t i t u t i v e models and numerical methods employing streamline-in­

tegration methods, which take into account the f a c t t h a t s t r e a m l i n e s

are r e a l c h a r a c t e r i s t i c s of the e q u a t i o n s . Therefore, we think that

most of the encountered instability problems are not related to instabilities in the underlying physical problem, but are caused by non-realistic fluid models.and/or deficiencies of the numerical method that is used. The aim of the work presented in this thesis is to gain more insight in these problems by:

a. analysing the basic properties of the system of equations depending on the various constitutive equations that are used, and

b. developing a numerical scheme for solving viscoelastic flows, aided by the results of a.

The outline of this thesis is as follows. In Chapter 2 we decribe the basic equations for the flow of viscoelastic fluids. In Chapter 3 we give an analysis of the basic properties of the equations : type, boun­ dary conditions and discontinuities. In Chapter 1 we present a numerical scheme based on the finite-element method and streamline in­ tegration. In Chapter 5 we compute a complex flow problem and analyse the convergence behaviour of the iterative methods to solve the equations. Finally, in Chapter 6 some concluding remarks are made.

2. GOVERNING EQUATIONS OF THE FLOW OF VISCOELASTIC FLUIDS.

The governing equations are the basic laws from continuum mechanics and the constitutive equation describing a particular fluid. In this chap­ ter we review the equations from continuum mechanics for non-polar media and consider several constitutive models in differential form. We discuss the behaviour of these models in steady linear flows. In sec­ tion 2. 4 we derive the dimensionless numbers governing the steady flow problem.

2.1 Basic equations from continuum mechanics.

We assume that the continuum hypothesis is valid. Therefore it is al­ lowed to introduce macroscopic field quantities defined pointwise, such as velocity and density. For a discussion of the continuum hypothesis see Batchelor (1967). Kuiken (1976).

The following principles from continuum mechanics are used for describing fluid flow (Truesdell, 1977) :

i. Conservation of mass, ii. Balance of linear momentum, iii. Balance of angular momentum,

iv. Balance of energy.

In order to simplify the basic equations, we make the following as­ sumptions :

- incompressible flow, - non-polar media,

- all variables and parameters in i., ii. and iii. are independent of temperature.

The last assumption leads to a decoupling of the balance of energy from the other equations. Hence, the balance of energy can be discarded from the system of equations. However, this assumption is made for reasons of simplicity only; in many practical problems thermal effects have to be taken into account.

With the above a s s u m p t i o n s , the equations in an inertial frame of reference take the following form :

conservation of mass:

divv = 0 in p., (2.1 .1)

balance of linear momentum:

pv = p£ + divo in n, (2.1.2)

balance of angular momentum:

o = gT in fi, (2.1.3)

where Q is a fixed bounded subspace of the Euclidian space, p is the density, v is the velocity, f is a body force per unit mass and o is the Cauchy stress tensor. A dot above a symbol denotes the material derivative, defined by :

C > - 75t< > - f t

(

>

+

I-8

rad

( >■ (2.1.1)

Constitutive equations specifying a should be such that (2.1.3) is satisfied; henceforth this will be assumed.

The incompressibility condition (2.1.1) is an internal constraint on the deformation of the medium. In order to satisfy this constraint, a hydrostatic stress tensor p1_ is necessary (Truesdell and Noll, 1965). Therefore, the Cauchy stress tensor o is decomposed into a hydrostatic part and a deviatoric part :

o = - p ^ + t, (2.1.5)

where p i s the h y d r o s t a t i c pressure and t the stess tensor caused by

the flow, the so-called e x t r a - s t r e s s tensor. In equilibrium the tensor

t vanishes.

Substitution of (2.1.5) into (2.1.2) yields for the balance of linear momentum :

pv + gradp = pf + divt. (2.1.6)

The extra-stress tensor t is determined by the deformation history of a material particle and is specified by the constitutive equation of the particular fluid. Constitutive equations for viscoelastic fluids are discussed in the next section.

The above equations have to be supplemented by the boundary conditions which will be discussed in Chapter 3.

2.2 Constitutive equations for viscoelastic fluids.

In the literature many constitutive equations have been proposed. For an extensive discussion and derivation of models the reader is referred to books on this subject, for example Astarita and Marucci (1971), Bird et. al. (1977), Janeschitz-Kriegl (1983), Tanner (1985).

To get some feeling for linear viscoelastic behaviour, it is useful and instructive to first consider two classical models from the linear one-dimensional viscoelastic theory, namely the Maxwell and Jeffreys models:

Maxwell Ao + o = n<, (2.2.1)

Jeffreys a = n K + o., Ad. + o. = n<, (2.2.2)

where o and o, are stresses, K is a strain, n and n are viscosities

1 s and A is a relaxation time. The mechanical representation of these

models, consisting of springs and dashpots, are given in Figure 2.1. Displacement is the analogue of strain and force is the analogue of stress.

M a x w e l

J e f f r e y s

n/x

t^ww

- » . CT

n/x

CMAAArn

—tt

-*■ CT

Figure 2.1 Mechanical representation of the Maxwell and Jeffreys models. The elastic moduli of the springs and the vis­ cosities of the dashpots are indicated.

For large time scales (low frequencies, low deformation rates) both the Maxwell and the Jeffreys model have a viscous behaviour. For short time scales (high frequencies, high deformation rates) the Maxwell model shows an elastic behaviour and the Jeffreys model shows a viscous behaviour. This is the fundamental difference between the Maxwell and the Jeffreys models. Combining several elements in a parallel fashion, the so-called generalized Jeffreys model can be constructed:

V ' ^ V

Vk

V -

( 2 . 2 . 3 )

where o. , A. , n. are the s t r e s s , the relaxation time and the v i s c o s i t y

of mode k r e s p e c t i v e l y . The m e c h a n i c a l r e p r e s e n t a t i o n of the

generalized Jeffreys model i s shown in Figure 2 . 2 . The g e n e r a l i z e d

Maxwell model i s obtained by taking n

0 .

The Maxwell and J e f f r e y s models are c h a r a c t e r i z e d by a d i s c r e t e

r e l a x a t i o n spectrum ( a f i n i t e number of d i s t i n c t relaxation times A

k

,

k = 1 , 2 , . . . , n ) . In many models for viscoelastic f l u i d s a small number

of d i s t i n c t r e l a x a t i o n mechanisms i s assumed. These models have the

form of a raulti-dimensional ( g e n e r a l i z e d ) Maxwell or J e f f r e y s model

extended with non-linear terms.

-tt^WVH

'n n n

t^WVW

Figure 2.2 Mechanical representation of the g e n e r a l i z e d J e f f r e y s

model.

Truesdell and Noll (1965) have given a g e n e r a l theory of c o n s t i t i v e

e q u a t i o n s in which the e x t r a - s t r e s s tensor t may be expressed as a

functional of the e n t i r e deformation history of a p a r t i c l e ,

t ( t ) = H ( C

t

(v) ),

( 2 . 2 . 4 )

where t i s the current time and C.(v) i s a tensor describing the defor­

mation of a m a t e r i a l p a r t i c l e at the time v r e l a t i v e to the time t .

Although i t i s possible to develop c o n s t i t u t i v e e q u a t i o n s from t h i s

general approach, one often o b t a i n s impractical equations with many

coefficients. Therefore, most c o n s t i t u t i v e models are derived from

other p r i n c i p l e s , which w i l l not be discussed h e r e . The r e s u l t i n g

constitutive equations, which s t i l l have to satisfy (2.2.H), may be in

the form of integral equations or d i f f e r e n t i a l equations. These models

are called integral and differential models respectively.

I n t h i s t h e s i s we s h a l l consider only d i f f e r e n t i a l models t h a t can be w r i t t e n i n the f o l l o w i n g form : k=n t = 2 n d + T , 1 = 1 I. , ( 2 . 2 . 5 ) s= = - k = 1- K where n i s a v i s c o s i t y (n £ 0 ) , d i s t h e s y m m e t r i c p a r t of t h e s T s

v e l o c i t y g r a d i e n t h = gradv_ and x. i s the e x t r a - s t r e s s t e n s o r of mode k . The s t r e s s t e n s o r s T. , k = 1 , 2 , . . , n , e a c h s a t i f y a d i f f e r e n t i a l e q u a t i o n with a t l e a s t two parameters : a r e l a x a t i o n time A (A 2 0 ) and a v i s c o s i t y n. (n. 2 0 ) . We s h a l l give some e x a m p l e s i n t h e f o l ­ l o w i n g . The modes a r e ordered such t h a t A > A for i < j . For b r e v i t y , we s h a l l often omit k = 1,2 n. For example, we s h a l l often w r i t e T. f o r the complete s e t of t e n s o r s x. , k = 1,2 n .

There i s a fundamental d i f f e r e n c e between n * 0 and n = 0 . If n * 0, s s s t h e ( d i f f e r e n t i a l ) model w i l l be s a i d t o be of J e f f r e y s or J type and i f n = 0 of Maxwell or M t y p e . In Chapter 3 we w i l l show t h a t t h e a d ­ d i t i o n of t h e term 2n d i n e q u a t i o n ( 2 . 2 . 5 ) changes the type of the system of e q u a t i o n s . The n e c e s s i t y of i n c l u d i n g the Newtonian term 2ri d i s most o b v i o u s f o r m o d e l l i n g d i l u t e polymer s o l u t i o n s . In t h a t c a s e , n accounts for t h e v i s c o s i t y of t h e Newtonian s o l v e n t . For polymer m e l t s e q u a t i o n ( 2 . 2 . 5 ) w i t h n * 0 can a l s o be u s e f u l for a l i m i t e d range of s h e a r i n g r a t e s , in o r d e r t o take i n t o account small r e l a x a t i o n t i m e s . T h i s may be c l a r i f i e d a s f o l l o w s . Polymer l i q u i d s c o n s i s t of networks of molecular c h a i n s with a l a r g e r a n g e of l e n g t h s c a l e s and t h e r e f o r e r e l a x a t i o n t i m e s , spanning s e v e r a l o r d e r s of magnitude. Let T be a c h a r a c t e r i s t i c t i m e of t h e f l o w . R e l a x a t i o n mechanisms h a v i n g r e l a x a t i o n t i m e s much s m a l l e r t h a n T b e h a v e Newtonian and t h e c o n t r i b u t i o n t o t h e t o t a l s t r e s s of t h e s e m e c h a n i s m s can b e a p ­ proximated by a s i n g l e term 2n d , n * 0. Depending on t h e value of T, n can be v e r y s m a l l w i t h r e s p e c t t o t h e s h e a r v i s c o s i t y f o r z e r o s t r a i n r a t e n . For T approximately equal or s m a l l e r than the s m a l l e s t r e l a x a t i o n t i m e s , n should be taken zero (M type model).

We s h a l l now g i v e some examples of c o n s t i t u t i v e models for the s t r e s s t e n s o r s x . .

i. Quasi-linear models.

These models are of the following form:

x k!k + ïk = 2 nk = ' k = 1'2 n' (2.2.6)

where T. is a derivative of T. given by

ï

k

- ï

k

"

h

- ï

k

"

Ï K - 5

•

h = h - 5d = (1 - |)h - | hT, 0 2 5 ^ 2 .

The following special cases are well-known classical models :

5 = 0 : upper-convected or B models, (2.2.7) * =k ïk = 6t ïk " =-ïk " , T " = k ' = (2.2.8) 5 = 1 : corotational models, * 6=k T ïk = 6T = ïk " ï'-k " ïk'ï ' ( 2-2"9 ) 5 = 2 : lower-convected or A models, * 6ïk T T = — - = T + h .X + T .h, (2.2.10) ïk St ïk " =k ïk "' K '

Often used i n n u m e r i c a l s i m u l a t i o n s a r e t h e upper-convected Maxwell (UCM) model (5 = 0, n = 1, n = 0 ) and t h e s o - c a l l e d O l d r o y d - B model ( 5 = 0 , n = 1, n * 0 ) . Note t h a t t h e f i r s t i s of Maxwell type and t h e l a t t e r of J e f f r e y s t y p e . Other models t h a t w i l l be r e f e r r e d t o i n t h e next chapter a r e t h e lower-convected Maxwell (LCM) model ( 5 = 2 , n = 1, n = 0 ) and t h e c o r o t a t i o n a l Maxwell (CRM) m o d e l ( 5 = 1 , n = 1 , n = 0 ) . The model (5 * 0 , 1 , 2 , n = 1, n = 0) i s sometimes c a l l e d t h e

s s Johnson-Segalman model ( s e e Tanner 1 9 8 5 ) . T h e s e m o d e l s a r e u s u a l l y

mechanics (e.g. frame indifference). Therefore, these models have only physical relevance for small deformation rates.

ii. Leonov and Giesekus models.

The Leonov model (Leonov, 1976) is a model based on the thermodynamics of irreversible processes and some additional assumptions. For plane flow the model may be written as follows

*k 6T

+ T

=k

+

2n~ik "

2

\i'

k

"

1

'

2 n

'

(2

'

2

-

11)

where 6T, / 6 t i s t h e f r a m e - i n d i f f e r e n t t i m e - d e r i v a t i v e defined in

_=k

(2.2.8). The s t r e s s tensor T^ i s subjected to the following constraint

A

k

det(1 + —T, ) = 1. (2.2.12)

■ v

k

In Appendix A we show that (2.2.12) i s satisfied at any time t if i t i s

s a t i s f i e d i n i t i a l l y .

The Giesekus model ( G i e s e k u s , 1982a,b) in a multi-mode form can be

written as follows

A

k W-

+

Jk

+

- V i k

= 2

V1.

k

- 1 . 2 . . . n , (2.2.13)

where a. (0 £ a. < 1) i s the so-called mobility parameter. This model

i s based on microstructural considerations and has not yet been f u l l y

examined for i t s usefulness in various flow configurations. Note that

for a. = 0 the model equals the upper-convected models given by (2.2.6)

1

with ( 2 . 2 . 8 ) , and for a. = -= i t equals the Leonov model in plane flow.

We have formulated the Leonov and the Giesekus model in terras of the

e x t r a - s t r e s s t e n s o r s T, . However, both Leonov and Giesekus have

o r i g i n a l l y formulated t h e i r models in terms of s o - c a l l e d i n t e r n a l

deformation tensors, which we shall discuss b r i e f l y . The Leonov and the

Giesekus models will be called internal deformation models.

The s t r e s s e s in e l a s t i c materials are determined by the deformation of

the material p a r t i c l e s r e l a t i v e to a fixed s t r e s s - f r e e s t a t e . E l a s t i c

materials have an ideal memory, i . e . if external loads are removed from

the e l a s t i c body, the body always r e t u r n s t o the fixed s t r e s s - f r e e

s t a t e . The deformation of the m a t e r i a l p a r t i c l e s r e l a t i v e to the

s t r e s s - f r e e s t a t e may be described by the Finger deformation tensor b,

satisfying the equation

5b

-£ = 0, (2.2.11)

ot

-where b = 1_ for the s t r e s s - f r e e s t a t e . The tensor b is symmetric and

positive definite (Truesdell, 1977). A well-known model for incompres­

s i b l e e l a s t i c materials i s the neo-Hookean model, derived from network

theory (Treolar, 1975) given 'by

t = y(b - n , u > 0, (2.2.15)

where y is an e l a s t i c modulus. Note, t h a t for b = 1^ the body i s s t r e s s

free.

Viscoelastic fluids do not have an ideal memory. When e x t e r n a l loads

a r e removed form t h e body, t h e r e i s a r e l e a s e of r e s i d u a l (or

recoverable) e l a s t i c deformation ( ' r e c o i l ' ) , but the r e l a x e d ( s t r e s s

-free) s t a t e i s not the original s t a t e . This may be i l l u s t r a t e d by the

d a s h p o t - s p r i n g model in Figure 2 . 1 . If a i s removed, t h e s p r i n g

r e l e a s e s e l a s t i c energy, but the system will not return t o the s t a t e

with zero s t r a i n K.

I n t e r n a l deformation models a r e based on t h e i d e a t h a t the e x t r a

-stresses for viscoelastic fluids are s t i l l determined by e l a s t i c defor­

mation and can, for example, be written in the form (2.2.15). However,

the e l a s t i c deformation i s not the t o t a l deformation, but the r e s i d u a l

e l a s t i c deformation ( i n t e r n a l d e f o r m a t i o n ) . The deformation tensor

describing the residual deformation i s not r e l a t i v e t o the o r i g i n a l

s t a t e , but r e l a t i v e t o an imaginary intermediate s t a t e , which can be

found a f t e r r e l a x a t i o n of the c u r r e n t s t r e s s s t a t e . F u r t h e r m o r e ,

analogous to the generalized linear models shown in Figure 2.2, i t i s

assumed that we may generalize the idea to models with multiple modes,

each having a d i f f e r e n t intermediate s t a t e . The Giesekus model can be

written as follows:

\

=k

=

F"

(

=k " V '

k = 1 , 2 n

' \ * °' (2.2.16)

k

where b is the internal deformation tensor, representing the residual elastic deformation of mode k. The tensor b, can be found from

=k

\W

+

h ' i

+ a

k

( b

= k - V

2

= 2 - (2.2.17)

I t i s e a s i l y shown, that substitution of (2.2.16) into (2.2.17) yields

equation ( 2 . 2 . 1 3 ) . The Leonov model for p l a n e flow i s found by

substitution of a. = -=. The tensors b are related to the Finger tensor

and should also be positive d e f i n i t e . However, comparing (2.2.11) and

( 2 . 2 . 1 7 ) , we see t h a t contrary t o the Finger tensor b, the value of

6b /fit does not vanish but depends on b . In Appendix A we show t h a t

t h e t e n s o r s b a r e positive definite for plane flow of the Leonov and

Giesekus models if proper i n i t i a l conditions are specified.

2.3 Plane steady linear flow and viscometric functions.

In the following chapters, we mainly consider two-dimensional (plane)

f l o w s . In order t o study the behaviour of various constitutive models

in plane flows, we consider the case of a steady linear flow in an i n ­

f i n i t e domain:

This solution satisfies the continuity equation (2.1.1). Substitution of (2.3.1) into the momentum equation (2.1.6), and using (2.1.4), gives

pH.V + pH.H.x + gradp = pf. (2.3.2)

For a two-dimensional incompressible flow it is easily shown by writing into components and using trH = 0, that H.H = (-detH)K Hence, equation

(2.3.2) can be written as follows

grad(px.H.V - |pdetH||x||2 + p + pef) = 0, (2.3-3)

where we have assumed that f is a conservative force: f_ = - pgrade . Integration of (2.3.3) yields

p = -px.H.V + lpdetH||x||2 - pef + P, (2.3.4)

where P is an integration constant. Hence, a velocity field _v and a stress field t given by (2.3-1) and a pressure field given by (2.3-4) satisfy the equation of motion. However, nothing is said about the stability of the flow or whether the flow is physically acceptable. These items require separate analysis for a particular constitutive model, which gives the relation between T and H.

Examples of steady linear flows are

plane elongation H__ = èe e - èe e„, V = 0 , (2.3.5)

— ril —1 —1 "~IL—c. — —

simple shear H„_ = q£,£2' 1 = 2 • (2.3.6)

rotation HD = 6e,e„ - Se„e. , V = 0 , (2.3.7)

where ê is the stretching rate, q is the shear rate, 6 is the angular velocity and e_ , e_. are the unit base vectors of an orthogonal Cartesian coordinate system x , x . The streamlines of the three flows are shown in Figure 2.3. The deformation of a material particle along its trajectory, which is initially square , has also been indicated.

plane elongational flow simple shear flow

rotational flow

Figure 2.3 S t r e a m l i n e s for plane e l o n g a t i o n a l flow (è > 0 ) ,

simple shear flow (q > 0) and r o t a t i o n a l flow (6 > 0 ) .

I t is useful to i n t r o d u c e a c l a s s i f i c a t i o n of plane steady l i n e a r

flows. Let us decompose the velocity gradient into a symmetric part D

and an anti-symmetric part ft :

H = D + ft, (2.3.8)

I T 1 T

where D = ^(H + H ) and ft = -x(h - H ) . Introduce the deformation rate D

= 2 = = = 2 = = and vorticity ft according to :

D = /-detD, (2.3.9)

If m , rn? are the principal directions of D and e. _{e? f o r m an}

( a r b i t r a r y ) o r t h o g o n a l p a i r of u n i t base v e c t o r s , then D and ft can be w r i t t e n as f o l l o w s :

p = ±(Dm m.- Dnyrip),

ft = ±(fte e_2- fte e1 ) ,

( 2 . 3 . 1 1 )

( 2 . 3 . 1 2 )

where t h e s i g n s a r e such t h a t D i 0 and f! £ 0. S u b s t i t u t i o n of t h e s e e q u a t i o n s i n t o ( 2 . 3 - 7 ) and comparison with t h e s t a n d a r d flows ( 2 . 3 . 5 ) , ( 2 . 3 . 6 ) and ( 2 . 3 7 ) shows t h a t D * 0, JJ = 0 r e p r e s e n t s a p l a n e e l o n -g a t i o n a l f l o w , D = ft * 0 a s i m p l e s h e a r f l o w and D = 0 , ft * 0 a r o t a t i o n a l flow. Therefore, i t i s u s e f u l to d i s t i n g u i s h t h e f i v e cases r e p r e s e n t e d i n Figure 2 . 4 . I: elongation II: dominance of elongation III: shear IV: dominance of rotation V: rotation

Q

Figure 2.4 C l a s s i f i c a t i o n of two-dimensional steady l i n e a r f l o w s .

The c l a s s i f i c a t i o n g i v e n above i s f o r s t e a d y l i n e a r f l o w s o n l y . However, i t may be useful f o r o t h e r t y p e s of f l o w a s w e l l . In t h e s e c a s e s , t h e c l a s s i f i c a t i o n i s only l o c a l and H should be r e p l a c e d by t h e

T

l o c a l value h = (gradv) and D2 and ft2 by d2 = - d e t d and ID2 = d e t u , u =

1 T = = =

If the constant e x t r a - s t r e s s t e n s o r in plane elongation and simple

shear are c a l l e d T

p

„ and T r e s p e c t i v e l y , then the following v i s ­

cometric functions are defined

TPE11 - TPE22 TSS12 _ _ . . ,

nPE 1 ' nSS " I - ' N1 = TSS11 TSS22 ( 2-3-1 3 )

where n_„ i s the plane elongational v i s c o s i t y , n

0 0

i s the shear v i s

-c o s i t y and N. i s the f i r s t normal-stress differen-ce. The se-cond normal

s t r e s s difference N_ is equal to 1„„

2?

, because T „ „ „ is zero for plane

flow. The v i s c o m e t r i c f u n c t i o n s n , n , N , N. of a constitutive

model can be used to f i t experimental v i s c o m e t r i c data of the f l u i d .

The function r\ i s not used (and measured) very often. Instead, the

elongational viscosity n

r

i s used, which i s defined s i m i l a r l y t o n ,

E. r e.

b u t now f o r t h r e e - d i m e n s i o n a l f l o w s . In t h a t c a s e , t h e steady l i n e a r flow with H„ = ê e . e . - -^e(e_e- + e _ e . ) i s c o n s i d e r e d .

=E —1—1 2 —2—2 — 3—3

The s h e a r v i s c o s i t y for z e r o s t r a i n r a t e ru = lira riqpi "iay be expressed in the v i s c o s i t y p a r a m e t e r s n and n by

n

n0 = ng + I nk, ( 2 . 3 . 1 4 )

u s k = l K

which i s valid for a l l the models to be considered in t h i s t h e s i s .

We s h a l l now c o n s i d e r plane steady l i n e a r flows and v i s c o m e t r i c

functions for the constitutive models given in section 2.2.

i . Quasi-linear models.

With the r e s u l t of appendix B, we can show that a plane constant s t r e s s

field T can be derived from (2.3-D and (2.2.6) if :

Fk = 1 - 1 X * ( 1 - C)2D2 + JU^n2 * 0, k = 1 , 2 , . . . , n ( 2 . 3 . 1 5 )

where D and f! are defined in (2.3.9) and (2.3.10), respectively. In Figure 2.5, the curve F. = 0 is shown. This curve divides the (!5,D)

plane into two regions. On the dividing curve the constant stress field T does not exist. Hence, the stress solution beyond this (singular) curve is unacceptable from a physical point of view. Note that the dividing curve moves to infinity for the corotational models (5 = 1) and asymptotically approaches the line D = (1 (simple shear) for the upper-convected models (£ = 0).

D = O

u n a c c e p t a b l e s o l u t i o n

1 - i I D = n

Q

Figure 2.5 (fl,D) diagram with r e s p e c t t o mode k of the steady

linear flow of a quasi-linear model.

The agreement of quasilinear models with experiments in various v i s

-cometric flows of polymeric fluids i s only moderately good. Fitting one

of the viscometric functions leads to a bad approximation of the other

functions.

In Figure 2.6 the viscometric function Ti

nc

. i s shown for the Oldroyd-B

model. Note the s i n g u l a r i t y in the function n

p F

. which denotes the

entrance of the unacceptable region in Figure 2 . 5 . Tanner (1985) gives

a short review of the very few biaxial stretching experiments that have

been carried out. These experiments indicate that for the polymers that

are used, r\ i s a decreasing function of é. I t i s not r e p o r t e d t h a t

PE

singularity in n

Ep

and generally for the unacceptable region of Figure

2 . 5 .

Figure 2.6 Viscometric function n

D

„ of the 01droyd-B model.

Because of the p o s s i b i l i t y of singular behaviour in steady linear flows

and the unstable behaviour in other flows (see Chapter 1 ) , the q u a s i

-l i n e a r mode-ls a r e not v e r y s u i t a b -l e for numerica-l s i m u -l a t i o n s .

Unfortunately, the majority of the numerical schemes in the l i t e r a t u r e

use q u a s i - l i n e a r models. That is why we have included these models and

show that other (more r e a l i s t i c , stable) models should be used.

i i . Leonov and Giesekus models.

In Appendix B we have shown that a plane constant s t r e s s field T can be

derived from ( 2 . 3 . 1 ) and ( 2 . 2 . 1 3 ) for any velocity gradient H. Good

f i t s for the viscometric functions n„q

ar

>d N. can be made with t h e s e

models. For other v i s c o m e t r i c functions, p a r t i c u l a r l y n , the models

h

are less s a t i s f a c t o r y . The viscometic functions n

D

ir

a n d n

oo of a

one-mode Giesekus one-model are shown in Figure 2.7 for 0 < a. < -=, a. = ■=

(Leonov) and ^ < o < 1. Note that contrary to the quasi-linear models,

t h e Giesekus and Leonov models have r e g u l a r viscometric functions.

Hence, from a numerical point of view t h e s e models are b e t t e r s u i t e d for numerical simulations.

Figure 2.7 Viscometric f u n c t i o n s r\ and nc o for the one-mode

re, OO

Giesekus model.

The results of this section concerning the regularity of T in a linear steady flow, equally apply to small regions in a steady flow around points with y_ = 0^ (fixed walls, center of a region of recirculation, interior stagnation point). In that case, the tensor H should be replaced by the local value of the tensor h.

2 . 4 Remarks on t h e system of e q u a t i o n s . Combining e q u a t i o n s ( 2 . 1 . 1 ) , ( 2 . 1 . 6 ) , ( 2 . 2 . 5 ) , ( 2 . 2 . 6 ) , ( 2 . 2 . 7 ) and ( 2 . 2 . 1 3 ) ! we have t h e f o l l o w i n g system of e q u a t i o n s i n t h e r e g i o n il : n pv^ + gradp = pf_ + n Av + £ d i v i . . ( 2 . 1 . 1 a ) 3 k=1 = K divv = 0, (2.U.1b)

V i

k - s.sk - ik-ST) + ik + !

^ i k = 2

\ ^ (2-4-1c>

where A = d i v g r a d and h = h f o r a * 0. Note, t h a t equation ( 2 . 4 . 1 ) i n c l u d e s a l l c o n s t i t u t i v e e q u a t i o n s t h a t have been d i s c u s s e d i n s e c t i o n 2 . 2 . I n t r o d u c t i o n of V = V V ' , p = - j - p ' , Tk = — Tk. h - Eh ' , V 1 1 1 d = 7 d ' , A = j -2 A', d i v = - d i v ' , grad = - g r a d ' , ( 2 . 4 . 2 )

where L is a characteristic length and V is a characteristic velocity, gives the following dimensionless form of (2.4.1):

ns n \ Rev'.grad'v' + grad'p' = —A'v' + £ —div'x' ,

~ n0 k-1 n0 = K

div'v'= 0, (2.4.3)

Dek(v'.grad'x|J. - h'.jk - Ik.h'T) + T^ + o ^ D e ^2 = 2d',

k = 1,2,...,n,

where we have assumed f_ - £ and steady flow. The Reynolds number Re and the Deborah numbers De. are defined by

R e - H t , (2.4.4)

V

Dek = ^ - , k = 1,2 n. (2.4.5)

The length L is a characteristic measure for the size of the flow region (e.g. the width of a channel) or of an obstacle in the flow. The velocity V is a characteristic measure for the velocity (e.g. the averaged velocity across the width of a channel).

We see from (2.4.3) that the steady flow problem is governed by the following dimensionless numbers

Re, Dek, nk/n0> ak> k = 1,2 n. (2.4.6)

For a particular fluid model (£, n , n , a , p have been specified) only two independent dimensionless numbers remain, for which it is com­ mon to choose the Reynolds number Re and the Deborah number De based on a mean relaxation time A :

N n

A0 = l i m 2ÏT— = ( l Xknk) / n0 ' (2-"-7)

q-0 ^qTSS12 k=1 K K U

De = -£-, (2.4.8)

where the expression for X. is valid for all models we have discussed. For the definition of T„„._ and N. in a simple shear flow, see section 2.3.

Equation (2.4.3) shows that for De << 1, the flow actually behaves Newtonian. In this case, the Reynolds number Re is the magnitude of the inertia terms relative to the viscous terms.

Although (2.4.3) is valid for any positive value of Re and De , it does not show explicitly the relative importance of the various terms in flows for high Deborah numbers. In that case, the stresses T. do not show viscous behaviour and t' is not 0 ( 1 ) , but may be of much lower

o r d e r , depending on the particular model and the s t r e s s component that

i s considered. Hence, the i n e r t i a terms and the extra-viscous terms may

be more important than i n d i c a t e d by Re and n Au respectively. If we

take the Leonov model as an example, then in steady shear flow the

shear s t r e s s of mode k i s approximately equal to n . A

k

for De. >> 1 , so

that the dimensionless shear s t r e s s of mode k i s of 0(De. ) for De >>

1. In t h a t c a s e , t h e magnitude of the i n e r t i a terms r e l a t i v e to the

t o t a l shear s t r e s s terms may be estimated by

n V n n n

n n

u i

i n e r t i a / s h e a r = pV

2

/[-j- + I -r^) - Re/(— + I — De, ) , (2.4.9)

L k=1Ak n0 k-I^O k

which can be much larger than Re. In the next chapter we will show that

the r a t i o given in equation ( 2 . 4 . 9 ) plays an important role in the

determination of the type of the p a r t i a l differential equations (2.4.1)

for Maxwell-type models (n = 0 ) .

The magnitude of the extra-viscous terms r e l a t i v e t o the v i s c o e l a s t i c

shear s t r e s s terms becomes

n _

1

extra-viscous/viscoelastic shear = n /( I n, De,. ] , (2.4.10)

s \ . k k '

which shows that the extra-viscous terra may be more important for high

Deborah numbers than suggested by the r a t i o n / n - in ( 2 . 4 . 3 ) . The

r e l a t i v e importance of various terms i s very much dependent on the par­

t i c u l a r v i s c o e l a s t i c model and one should be very c a r e f u l in e . g .

neglecting the i n e r t i a terms on the basis of Re << 1.

Near geometrical d i s c o n t i n u i t i e s ( e . g . sharp c o r n e r s ) , the v e l o c i t y

g r a d i e n t s are much l a r g e r than i n d i c a t e d by V/L. The velocity V and

length L have t o be r e p l a c e d by l o c a l values V and L in ( 2 . 4 . 3 ) ,

(2.4.4) and (2.4.5) with usually

Note, that although V < V the local velocity gradients are much l a r g e r .

Also, the local Reynolds number i s smaller and the l o c a l Deborah num­

bers are l a r g e r than the global numbers. Hence, near such discon­

t i n u i t i e s the flow i s much more ' e l a s t i c ' than predicted by the global

Deborah n u m b e r s . If t h e v e l o c i t y g r a d i e n t s a r e s i n g u l a r ( e . g .

| |gradv| | •+ » if x ■* cornerpoint) the flow i s always highly e l a s t i c

near the s i n g u l a r p o i n t , even if the global Deborah numbers De are

very small.

2.5 Concluding remarks.

In t h i s chapter, we have considered the basic equations for v i s c o e l a s

-t i c f l u i d flow, including various cons-ti-tu-tive models. We have s-tudied

the behaviour in plane l i n e a r steady flows and have shown t h a t the

s t r e s s tensor can become s i n g u l a r for q u a s i - l i n e a r models. Hence,

quasi-linear models are not very a p p e a l i n g , because d i f f i c u l t i e s a r e

expected i n numerical simulations. Models having a regular s t r e s s ten­

sor in linear steady flows are preferred. Examples of such models a r e

t h e Leonov and G i e s e k u s m o d e l s . A l t h o u g h t h e s e models are not

completely satisfactory from a physical point of view, we w i l l not d i s ­

cuss more s u i t a b l e models and use only the Leonov and Giesekus models

in our numerical simulations. The quasi-linear models will only be used

to show the d i s t i n c t advantages of the Leonov and Giesekus models.

In the next chapter, we will discuss various mathematical properties of

the system of equations.

3 . ANALYSIS OF THE SYSTEM OF THE EQUATIONS.

3.1 I n t r o d u c t i o n .

Most of the numerical methods for v i s c o e l a s t i c flow have c o n v e r g e n c e problems a t high e l a s t i c i t y . I t i s p o s s i b l e t h a t some of t h e s e problems a r e r e l a t e d t o t h e b a s i c p r o p e r t i e s of t h e e q u a t i o n s , which a r e not w e l l u n d e r s t o o d . For example, Joseph, Renardy and Saut (1985) suggest t h a t some problems might be r e l a t e d t o change of type of the e q u a t i o n s from e l l i p t i c t o h y p e r b o l i c , s i m i l a r t o t h e s i t u a t i o n i n gas dynamics.

T h e r e a r e o n l y few p a p e r s t h a t pay a t t e n t i o n t o t h e m a t h e m a t i c a l p r o p e r t i e s of the e q u a t i o n s for v i s c o e l a s t i c flow. However, s t i m u l a t e d by the d i f f i c u l t i e s i n f i n d i n g s u i t a b l e numerical m e t h o d s , t h e number of such p a p e r s i s g r o w i n g . R u t k e v i c h and R e g i r e r ( 1 9 6 8 ) , Rutkevich ( 1 9 6 9 ) and R u t k e v i c h (1970) have been t h e f i r s t t o a n a l y s e b a s i c p r o p e r t i e s such a s t y p e of t h e e q u a t i o n s and t h e n a t u r e of p o s s i b l e d i s c o n t i n u i t i e s . Ultman and Denn (1970) e x p l a i n t h e o n s e t of some a n o m a l o u s h e a t and momentum t r a n s f e r by a c h a n g e of t y p e of t h e e q u a t i o n s .

R e c e n t l y , t h e r e h a s been more i n t e r e s t i n t h i s f i e l d and some new papers have been p u b l i s h e d . Luskin (1981) a n a l y s e s a s p e c i a l l i n e a r i z e d form of the upper-convected Maxwell model by t r a n s f o r m i n g t h e system to a c a n o n i c a l form. V.d. Zanden e t . a l . (1984) determine t h e t y p e of the s y s t e m of e q u a t i o n s f o r an u p p e r - c o n v e c t e d Maxwell model for steady f l o w . J o s e p h , Renardy and S a u t (1985) e x t e n d t h i s f o r a c l a s s of O l d r o y d m o d e l s . Yoo, Ahrens and Joseph (1985) c o n s i d e r t h e problem of s t e a d y f a s t flow of a family of Maxwell f l u i d s i n t o a h o l e . They show t h a t t h e flow i s p a r t i t i o n e d i n t o e l l i p t i c and h y p e r b o l i c r e g i o n s , s i m i l a r t o the s i t u a t i o n in g a s d y n a m i c s . Yoo and J o s e p h (1985) and A h r e n s , Yoo and J o s e p h (1987) a n a l y s e t h e flow of an upper-convected Maxwell model through a channel with wavy w a l l s . For high s p e e d s t e a d y f l o w , t h e e q u a t i o n governing the v o r t i c i t y has a h y p e r b o l i c c h a r a c t e r i n the c e n t r a l r e g i o n and an e l l i p t i c c h a r a c t e r i n t h e w a l l r e g i o n .

Renardy (1985) studies the existence of slow steady flow of viscoelas­

t i c fluids for homogeneous boundary conditions. Joseph and Saut (1986),

Dupret and Marchal (1986a,b) consider the change of type for various

viscoelastic fluid models.

In o r d e r t o o b t a i n a b e t t e r u n d e r s t a n d i n g of the mathematical

properties of the equations, t h e i r t y p e , the r e q u i r e d boundary con­

d i t i o n s and the nature of possible discontinuities are investigated in

t h i s chapter. Apart from the theoretical analyses of Rutkevieh (1969)

and Dupret and Marchal (1986a), we do not know of any a r t i c l e paying

attention to discontinuities in viscoelastic flow. Presumably, no ex­

perimental evidence that discontinuities occur e x i s t s so far.

In section 3.2 we analyse the convection equation, which resembles the

c o n s t i t u t i v e e q u a t i o n s , and show on the basis of that equation that

under certain conditions, i t i s likely that d i s c o n t i n u i t i e s may occur

in viscoelastic flows. We have included the analysis of discontinuities

in order to r a i s e discussions on whether they e x i s t or n o t , both in

r e a l flows and in the solution of the p a r t i a l d i f f e r e n t i a l equations.

Moreover, the jump conditions can provide information on c o r r e c t boun­

dary conditions.

3.2 Convection equation.

Before we proceed with the analysis of v i s c o e l a s t i c flow, we analyse a

special convection equation, in order to gain some i n s i g h t . We consider

t h e f o l l o w i n g e q u a t i o n , which shows much r e s e m b l a n c e to t h e

constitutive equations of a Maxwell-type or Jeffreys-type model:

At + T = f ( x ) , A > 0, (3.2.1)

where T i s the transported quantity (such as a component of the e x t r a

-s t r e -s -s ten-sor) and where a velocity field v(x) i -s given. We -shall only

consider plane flow. Writing the m a t e r i a l d e r i v a t i v e according to

Euler's formula, (3.2,1) becomes

*[fr + v.gradx] + T = f(x). (3.2.2)

o t — —

This is a time-dependent convection equation where T is a quantity both transported by the flow as well as influenced by a relaxation mechanism and a source f.

If the characteristic surface is written as $(x,y,t) = 0 we find for the characteristic equation of (3.2.2)

Q U ) = A(<(>t + v.grad*) = 0, (3-2.3)

which i s e a s i l y shown with the r e s u l t s of Appendix C. This means, t h a t p a r t i c l e t r a j e c t o r i e s of t h e b a s i c flow a r e c h a r a c t e r i s t i c c u r v e s . For s t e a d y x t h e s t r e a m l i n e s a r e c h a r a c t e r i s t i c s . The s o l u t i o n of ( 3 . 2 . 1 ) i s given by x ( x , t ) = x e x p [ - ( t - t )/X] + — o o (3.2.1)) t + Y j f e x p [ - ( t - t ' ) / X ] d t ' , fc0

where the integration should be carried out along particle trajectories and x is the value of x at time t £ t of a particle being at x at

o o ^ ° -time t . The -time t has t o be chosen such t h a t x i s known, or i f t h a t

o o i s not p o s s i b l e , t - = - » .

The c o n d i t i o n f o r a jump i n x a c r o s s a s u r f a c e <(>(x,y,t) = 0 i s given by ( s e e Appendix E ) :

X(<t>t + v . g r a d <(>)['x] = 0, ( 3 . 2 . 5 )

where the jump in x is denoted by [x]. From (3.2.5) it follows that discontinuities can only occur along characteristics.

Assume t h a t the given v e l o c i t y field v i s of the type t h a t occurs in the flow over a step (Figure 3.1)•

/ / / / / / / / / / / / / / / / / / / / /

; ; ; / ; ; / / > / ;

/ / / / / / / / / / /

V'

Figure 3-1 Flow over a step.

Let the initial condition at time t=0 be T=0, while the value of T at the inlet r is given as a function of time. The solution is given by (3-2.10 where tQ=0, TQ= 0 for particles that have been in the domain

from time = 0 to time = t. For particles that have entered the flow domain in this time period, time t (0 £ t & t) is the entry time at

+ ° + u

r and T is the specified value at V .

A recirculation zone exists. It is realistic to assume that the separation line between this zone and the main flow is attached to the sharp corner A, as shown in Figure 3-1. We have further assumed that the time spent at the corner A by particles that are on streamlines arbitrary close to the corner is finite. For example, this is the case if the radial velocity v at the corner behaves like

m v - r ,

r 0 < m < 1, (3.2.6)

where r i s the r a d i a l coordinate. Newtonian flow appears to be close to such behaviour (Moffatt, 1961).

Now, consider a streamline I in the mainflow and a streamline I I in the r e c i r c u l a t i o n zone (see Figure 3 . 1 ) . Since the time spent in the corner i s f i n i t e , the h i s t o r y of p a r t i c l e s a t both sides of the separation l i n e i s d i f f e r e n t . This means, t h a t the (history) i n t e g r a l in ( 3 . 2 . 4 ) i s d i f f e r e n t ; h e n c e , t h e value of T i s d i s c o n t i n u o u s a c r o s s the

separation l i n e . We may conclude, that the occurrence of discontinuous

solutions i s not unlikely for equations of type (3-2.1) with a velocity

f i e l d containing r e c i r c u l a t i o n zones, even if the boundary and i n i t i a l

c o n d i t i o n s are smooth. I t i s expected t h a t discontinuities may also

occur in real viscoelastic flows, where ( p a r t s of) the c o n s t i t u t i v e

e q u a t i o n s behave l i k e equation ( 3 . 2 . 1 ) . We n o t e t h a t recirculation

zones are quite common in domains with sharp c o r n e r s . Furthermore, we

n o t e t h a t van der Zanden e t . a l . (1984) have shown that the

upper-con-vected Maxwell model may have an inherent discontinuity of the boundary

s t r e s s e s at sharp corners.

3.3 Type of the system of equations.

Evaluation of the type of the system of p a r t i a l d i f f e r e n t i a l equations

i s essential for studying various properties of the equations (Courant

and H i l b e r t , 1962). Furthermore, the type can give us information on

which numerical methods should be used in order to make successful com­

putations.

The general form of the functional of the r e l a t i v e deformation tensor C

of a f l u i d p a r t i c l e , given in equation ( 2 . 2 . 4 ) , suggests a hyperbolic

character of the constitutive equations. This i s confirmed by the dif­

f e r e n t i a l equations when the velocity is supposed to be given: the par­

t i c l e t r a j e c t o r i e s are real c h a r a c t e r i s t i c s (or streamlines in the case

of steady flow) and the Riemann i n v a r i a n t s are composed of the com­

ponents of the e x t r a - s t r e s s tensor. However, for the coupled system of

e q u a t i o n s t h e r e i s a fundamental difference between models of Maxwell

type and of Jeffreys type.

I t a p p e a r s t h a t f o r M-type models t h e p a r t i c l e t r a j e c t o r i e s

(streamlines in steady flow) are real c h a r a c t e r i s t i c s of the system of

equations with a multiplicity that i s equal to the number of components

of T. minus 1. For J-type models t h i s m u l t i p l i c i t y i s always equal to

the number of s t r e s s components. The equations of M-type models can

change type due to flow conditions and J-type models cannot.

3.3.1 Maxwell-type_models.

The s y s t e m of e q u a t i o n s f o r M a x w e l l - t y p e m o d e l s h a s been g i v e n i n ( 2 . 4 . 1 ) w i t h n = 0 . We only c o n s i d e r p l a n e flow. The unknowns a r e u,

k k S k v, p , T , T , T and t h e s y s t e m of 3n + 3 e q u a t i o n s becomes f o r xx yy xy X k *0 : * k ^ k ,, - ^ - n 3T n 3T

<%

+

- t x

+

^

+

ft " l ^ -J

=

^ - PV (3.3.la)

*

+

P-g

+

-fy

*%

" l ^ - 1 ^ -PV

i

+

g = °' ".MO

3T 3T 3T , « , » , „ T a, xx ^ xx xx ^ k 3u _,_ k 3u . k 3v ^ xx k. , . . " 3 T + U~3x- + V- 3 F + 3x 3ÏÏ + ay 3 ^ + bx 3x + ~ + ^ k W °' 8 K - K ,. K K T 3T 3T , . , „ , . T a ( 3 - 3 . 1 d ) ¥■ + u ^ ♦ v ^ S + c,k | H ♦ dk | X + d,k | ï + ^ * ^ ( I k > y y = °> ( 3 - 3 . 1 e ) 3t 3x 3y y 3y x 3x y 3y Ak nk =k'yy 3T 3T 3T , . , . T a, ( 3 . 3 . I f ) 3t 3x 3y y 3y x 3x A ii »k xy k k k k k k k k

where the c o e f f i c i e n t s a , a , b , c , d , d , e and f are defined by

x ' y ' x ' y x ' y ' y x J

ak = - 2 ^ - 2 ( 1 - 0 4 , ak = dk = - (2 - O xk y.

bx = cy = < > dy = - 2 ^ - 2 ( 1 - ^ y y ' ( 3'3'2 )

ek

= - ^Ü - (1 - £ )

T k +

I

T k f k

= - \ - o - i

)xk +

I

T k

E q u a t i o n s ( 3 - 3 . 1 ) form a q u a s i - l i n e a r system of 3n + 3 p a r t i a l d i f ­ f e r e n t i a l e q u a t i o n s . The c h a r a c t e r i s t i c e q u a t i o n may be found by u s i n g t h e r e s u l t s of Appendix C. I f t h e c h a r a c t e r i s t i c s u r f a c e i s w r i t t e n as <f>(x,y,t) = 0, we f i n d

QUJ

/

B B . . I 0 , ( 3 . 3 - 3 )

where t h e 3X3 m a t r i c e s A, B, G. and I a r e defined by

G, = P<t> 0 * 0 p<j) <(i 0 * 0 x Ty a * _{x x} , B = -<j> 0 - * x yy 0 - * - * y x ( 3 . 3 . 4 ) + a % y y c * y y ki e d> y y b % X X dk* + dk* x x y y X X 0 0 0 , I = <f> 0 0 0 0 <(> 0 0 è w i t h <f> « <f> + u<J> + v<t> a n d w h e r e t h e C a r t e s i a n c o m p o n e n t s o f t h e n o r -' Srh Ï^A %A

mal vector are denoted by (*,*,$.) = ( ^ > ^ > ^ ) - After evaluating the x y t dx dy dt d e t e r m i n a n t , e q u a t i o n ( 3 . 3 - 3 ) r e d u c e s t o : Q C ^ ) =

13n-3

A n k=1 K ( 3 - 3 - 5 ) (<j>2 + <j.2)(a *2 vx Ty x xTx a d>2 + 2a * * - pd>2) = 0 , yy y xy x y r T

where a , a and a are the components of a tensor a defined by xx yy xy n

I

k=1

a = I a

k

. (3.3.6)

with k nk . £. k C k k _ k axx = Ak U 2; TX X " 2Tyy' xy " Tx y ' k k ,, g- k E, k a = T— + (1 - $ ) T - £ T . yy *„ 2 yy 2 xx (3.3-7)

An e q u a t i o n s i m i l a r t o ( 3 - 3 . 5 ) f o r a one-mode model h a s a l s o been d e r i v e d by v . d Zanden e t . a l . ( 1 9 8 4 ) , Joseph e t . a l . (1985) a n d Hulsen ( 1 9 8 6 ) .

The c h a r a c t e r i s t i c e q u a t i o n ( 3 . 3 - 5 ) i s composed of t h r e e f a c t o r s , which w i l l be t r e a t e d s e p a r a t e l y b o t h for u n s t e a d y and s t e a d y f l o w . In t h e l a t t e r case t h e r e a r e o n l y two independent v a r i a b l e s x and y . The c o r ­ r e s p o n d i n g c h a r a c t e r i s t i c e q u a t i o n can be found from ( 3 . 3 . 5 ) by p u t t i n g

♦t

=

°

-The f i r s t f a c t o r in ( 3 - 3 . 5 ) g i v e s :

<(. = <(.,. + u<j>x + v<j>y = 0 , ( 3 - 3 . 8 )

which shows that surfaces consisting of particle trajectories in the (x,y,t)-space are real characteristic surfaces of multiplicity 3n-1. In the steady case, streamlines (=particle trajectories) are real charac­ teristic curves of multiplicity 3n-1. It is remarkable that particle trajectories are real characteristics of multiplicity 3n-1 only, since 3n stress components are determined by the deformation history. The constitutive equation considered separately has 3n real characteristics which are all particle trajectories. The reduction in the number of real characteristics indicates that there is a strong coupling between the momentum equation and the constitutive equation.

The second f a c t o r in ( 3 . 3 . 5 ) y i e l d s :

<f>* + <t>y = 0 . ( 3 . 3 . 9 )

which l e a d s to t h e c h a r a c t e r i s t i c s u r f a c e <j> = <f> = 0, <j> ^ 0, which i s t h e x , y - p l a n e . This i s t y p i c a l for p a r a b o l i c e q u a t i o n s ( e . g . t h e h e a t e q u a t i o n ) . For t h e s t e a d y c a s e , e q u a t i o n ( 3 . 3 . 9 ) does not have r e a l s o l u t i o n s .

The most complicated f a c t o r of ( 3 . 3 . 5 ) i s the l a s t o n e , which g i v e s :

a <(>2 + a <(>2 + 2 a <t> <f> - p<J>2 = 0 . ( 3 - 3 . 1 0 )

x xTx y yry x yTxYy

T h i s e q u a t i o n i s the c h a r a c t e r i s t i c e q u a t i o n of the following p r o t o t y p e e q u a t i o n :

Ü Ü + a 3 2 w + ? 3'W = 2!«

*xx 3x2 yy ïy"2" axy 3x3y p Dt2

a „ „ ^ 7 ? + aw i, -r-i + 2 av i r T - T - = p TTTT . ( 3 - 3 . 1 1 )

I f the t e n s o r a i s p o s i t i v e d e f i n i t e , e q u a t i o n (3-3-11) i s a h y p e r b o l i c wave e q u a t i o n w i t h c o n v e c t i o n and n o n - i s o t r o p i c wave speed p r o p e r t i e s .

We can s i m p l i f y ( 3 . 3 . 1 0 ) by t r a n s f o r m a t i o n t o t h e p r i n c i p a l a x e s of a and by t h e i n t r o d u c t i o n of t h e p r i n c i p a l wave speeds c. and c_:

c2<(.2 + c2*2 - (4>t + vl<(.1 + v2*2)2 = 0 , ( 3 - 3 . 1 2 )

where the principal directions of the tensor a/p are denoted by the indices 1 and 2, and the corresponding principal values by c2 and c2,

hence

c2 = aj/p i = 1, 2, (3-3-13)

a b e i n g t h e p r i n c i p a l v a l u e s of a . R e a l - v a l u e d wave speeds a r e r e ­ q u i r e d f o r a w e l l - p o s e d i n i t i a l v a l u e problem f o r e q u a t i o n ( 3 - 3 - 1 1 )

(Courant and H i l b e r t , 1962). Equation (3-3.13) shows t h a t for c . t o be r e a l , i t i s r e q u i r e d t h a t a Z 0 , i = 1 , 2 .

If a. < 0 for at l e a s t one value of' i=1,2, an i n i t i a l value problem for

equation (3.3.11) i s not well-posed (ill-posed) and so-called Hadamard

i n s t a b i l i t i e s exist (Courant and Hilbert, 1962). These are very strong

i n s t a b i l i t i e s with an i n f i n i t e number of i n s t a b i l i t y modes, having a r ­

b i t r a r i l y short wavelenghts. The growth rates tend to i n f i n i t y as the

wavelengths tend to z e r o . Problems t h a t are i l l - p o s e d a r e normally

r e j e c t e d as p h y s i c a l models, as has been pointed out by Courant &

Hilbert (1962). However, the same authors give examples of meaningful

but ill-posed problems.

Recently, there has been some discussion on the p h y s i c a l relevance of

of the system of equations for viscoelastic flows becoming ill-posed.

For example, Ahrens, Joseph, Renardy and Renardy (1984) try t o explain

c e r t a i n abrupt changes in polymer flow ( e . g . melt f r a c t u r e ) with

i n s t a b i l i t i e s of Hadamard type. We agree with Joseph and Saut (1986)

that for some physical phenomena, Hadamard i n s t a b i l i t i e s can be related

to the onset of severe i n s t a b i l i t i e s in the physical problem. However,

we think t h a t i t i s s p e c u l a t i v e to explain i n s t a b i l i t i e s in flows of

viscoelastic fluids by the ill-posedness of a s e t of equations which

are based on such simple models as the one-mode quasi-linear ones.

These models have only a physical meaning for small deformation r a t e s .

Therefore, we regard loss of evolution as a model deficiency and agree

with Dupret and Marchal (1986b) that such models should not be used for

practical applications. Furthermore, we also believe that a. = 0, which

leads to a vanishing wave speed : c. = 0 , i s not a c c e p t a b l e . Equation

(3-3.10) represents wave phenomena related to the e l a s t i c i t y . From t h i s

viewpoint we think that a full "locking" of the e l a s t i c waves in the i

-direction i s also a model deficiency. The foregoing discussion leads to

the requirement that the tensor a should be positive d e f i n i t e .

Real f l u i d s such as polymer melts have many relaxation times, spanning

several orders of magnitude. The values of n. and A decrease for

k » 1,2 n . The value of n /A , which i s an e l a s t i c modulus,

increases for smaller A (for data see for example Tanner (1985)). From

the d e f i n i t i o n of the tensor a in equation (3-3.6) and (3-3-7), we see