IIJIT'IITIIIWmm}mmmflll‘mmWillI j

3‘ 1293 90071 7870

 

LIBRARY
Michigan State
University

 

 

 

This is to certify that the

dissertation entitled

A finite element approximation of the

flow in a toroidal tube.

presented by

Adel N. Boules

has been accepted towards fulfillment
of the requirements for

Ph.D. Mathematics

degree in

 

 

Q

dyflm;

fiajor professor /

Da Mr 7

MS U is an Affirmative Action/Equal Opportunity Institution 0-12771

 

 

MSU

LIBRARIES
.—_——.

 

 

 

RETURNING MATERIALS:
Place in book drop to
remove this checkout from
your record. FINES will
be charged if book is
returned after the date
stamped below.

 

 

35.5;‘1' l 4 l9???
./
.‘ if”, ZI? ‘

I7 ,

MAGIC 2

MAR 0 51999
r" @55-

 

 

A FINITE ELEMENT APPROXIMATION

OF THE FLOW IN A TOROIDAL TUBE

By

Adel Boules

AN ABSTRACT OF A DISSERTATION
SUBMITTED TO
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Mathematics

1987

ABSTRACT

A FINITE W APPROXIMATION

OF THE FLOW IN A TOROIDAL TUBE

BY

ADEL BOULES

The problem is to compute the velocity field of an incompressible
viscous fluid in a toroidal pipe under the effect of constant pressure
gradient in the axial direction. The complete set of momentum and
continuity equations were derived. The motion depends on two
parameters: The Reynolds number and curvature ratio.

An appropriate function space was introduced. and a variational
formulation of the problem in that space was defined.

A proof of existence and uniqueness of the solution was then
presented.

Using a finite element approximation. the continuous infinite
dimensional problem was replaced by a finite dimensional (yet
nonlinear) approximate problem.

Two methods were used to solve the finite element equations:

(a) A direct linearization iterative scheme which was shown to
converge to the solution of the approximate problem. This scheme has
the rather serious drawback that at each iteration one must compute and

solve a different set of unsynlnetric linear equations.

(b) A least squares formulation was presented. In this
formulation the solution of the problem is the minimizer of a certian
"cost" functional. A conjugate gradient scheme was found most
efficient in finding the minimizer of the cost functional and hence the
solution to the problem.

The great advantage of this scheme is that during the different
steps of the conjugate gradient scheme. one solves the same symmetric
positive definite banded matrix equation for different forcing terms.
Thus we need to compute and factor this constant matrix only once
during the entire computation.

Being symmetric and banded. the mtrix is easy to store compactly
thus allowing the use of a fine mesh. and consequently producing
greater accuracy. Because the untrix is also positive def inite. it can
be factored by the stable Choleski decomposition.

Both schemes were implemented and the velocity field was computed.
The two schemes produced almost identical results for all cases
computed. Excellent agreement is also achieved in comparison with the
results of previous investigators. I

A perturbation solution of the equations of motion was then
obtained for small curvature ratios and low Reynolds numbers. These
solutions also compare well with the numerical solutions described

above .

ACKNOWLEDGEMENTS

I would like to express my gratitude to Professor C.— Y. Vang. my
thesis adviser. for his patience and guidance during the course of this
work. I also wish to extend my thanks to the members of my thesis
committee. Professors D. Yen. J. Kurtz. B. Drachman.and C. MacCluer.

I should also acknowledge my gratitude to Mrs. B. Reiter. the
mathematics librarian. for her kind assistance. and to Miss S. Polomsky
for typing this dissertation.

Finally. I wish to express my sincerest thanks to the Case center
for Cbmputer Aided Engineering. for furnishing the computer equipment

without which this work would not have been possible.

DEDICATED TO

MYWIFEANDMYSON

ii

Chapter 1

1.1

1.2

1.3

1.4
Chapter 2

2.1

2.2

2.3

2.4
Chapter 3

3.1

3.2

3.3

3.4

3.5

3.6
Chapter 4
Chapter 5

5.1

5.2

5.3
Chapter 6

6.1

6.2

6.3

TABLE OF CONTENTS

Introduction and formulation ...........

Description of the problem ...................

Governing equations .................

Literature review ...................

Objectives and organization ..................

EXistence and uniqueness OOOOOOOOOOOOOOOOOOOOO

Preliminaries and notation ...................

Variational formulation of the problem .......

EXistence 0.0.0.000...OOOOOOOOOOOOOOOOOOOCOOOO

uniqueness 0.0.0.0....OOIOOOOOCOOOOOOOOO0.0...

The approximate problem and its solution .....

Motivation and statement of the problem ......

BreZZi-type condition OOOOOOOOOIOOOOOOOOOOOOOO

The EQUivalenCE Of (Ph) and (Ch) 0000000 0.0000

The solution of the approximate problem ......

Afixed paint SCheme OOOOOOOOOOOOOOOOOOOOOOOOO

The discrete Stokes problem ..................

Error anaIYSis O0..OICOOOOOOOOOOOOOOOOOOOO0.0.

Least squares approach .......................

IntrOdUCtion OOOOOCCOOOOOOOOOOOOOOOOOOOOOOOOOO

A least squares formulation ..................

A conjugate gradient scheme ..................

Numerical resu1ts 0....OOOOOOOOOOOOOOOOOOOOOOC

Computational procedure

some comparisons 0......OOOOOOOOOOOOOOOOOOOOC.

Results and discussion

iii

33
38

43
51
51

Chapter 7 A perturbation solution ......................
7.1 Introduction and formulation .................

7.2 Results and discussion .......................

Appendix A ..............................................

Appendixs 00.000.000.00.0..I0.00000000000000000000000000

List of references

iv

74

79

89
90
93

LIST OF FIGURES

FIGURE PAGE
1 Geometry and coordinate system ..................... 2
2 The reference elements ............................. 28
3 Mesh configurations ................................ 59
4 Contour lines of axial velocity at R=250 ........... 66

5 Secondary flow stream lines at R2250 ............... 67
6 Contour lines of axial velocity at R=500 ........... 68
7 Secondary flow stream lines at R2500 ............... 69
8 Contour lines of axial velocity at D=1000 .......... 70
9 Secondary flow stream lines at D-1000 .............. 71
10 Contour lines of axial velocity at 022000 .......... 72
11 Secondary flow stream lines at 022000 .............. 73
12 Geometry and coordinate system for elliptic cross
section ............................................ 74
13 Stream lines for a,= 0 (perturbation solution) .. 83
14 Stream lines for a = 45. (perturbation solution) .. 84
15 Stream lines for a = 90 (perturbation solution) .. 85
16 Effect of curvature on flow ratio (Rs6) ............ 86
17 Effect of curvature on flow ratio (R-20) ........... 87
18 Effect of curvature on flow ratio (R-SO) ........... 88

Table

LIST OF TABLES

Variation of u3max with mesh size .............. 60
Variation of u3(0) with mesh size .............. .60
Comparison of u3max ................... ... 62
Comparison of u3(0) ............................ 62
The flow ratio 0/ Q0 .......................... 65
The flow ratio 0/ QD ........................... 65
Coefficients in equation (7.17) ................ 81
Coefficients in equation (7.17) ................ 81

vi

CHAPTER 1: INTRODUCTION AND FORMULATION
1.1 Description of the problem

This work is concerned with the flow of an incompressible viscous
fluid in a toroidal pipe. Duct flows are important in many engineering
and biological applications. such as piping systems. the design of heat
exchangers. and the study of blood flow in arterial systems. The
geometry and coordinate system are shown in Figure 1. where O is the
center of the torus. and C is the center of the circular cross section.
Let one of the independent variables be the arc length X3 of the center
line of the pipe measured from a fixed point F. and Let T. N. B be unit
vectors along the tangent. normal and binomial of the center line. thus
N and B lie in the plane of the cross section. and for a typical point

Q in that plane we have
Q = xln + x213

The flow is assumed to be laminar. steady. fully developed. and
driven by a constant pressure gradient in the axial direction. We
therefore assume that the pressure P is a linear function of X3. and

that the velocity components are functions of X1 and X2 only.

 

 

~16

 

 

 

 

omoannn< mam ooonmwamnm m<mnma

mwosnm H

1.2 Governing equations

The major portion of this work is concerned with a circular cross
section. Let 2 be the radius of the torus. a be the radius of the
cross section. and let x = é-be the curvature of the tube. As was done
in Vang (1981). we computed the metric tensors (gij) and (gij). and the

Christoffel symbols F1 For a torus the system is orthogonal. The

jk'
metric tensors are
- 11 _ -1
O 1 O O 1 0
O 0 J O O 1/J
where J = (1 - KX1)2

There are only two non-zero Christoffel symbols

For j = 1. 2. 3. let uJ be the velocity component in the xJ direction
and let VJ be the tensorial components in the same directions. thus
V1 = U1. V2 = U2. V3 = U3 / «51 The covariant velocity derivatives are

given by

ik62v1(av“
V.j =3;E5;3 + --{F av“) - Fa'( k(ax7+ fivfi) + F; “(5x3 +r§bv“)

where summation over repeated indices is understood.

4

A list of the velocity derivatives Vijk is given in appendix A.

We now substitute the above in the tensorial form of the Navier-Stokes

equations:

iav1 iva _g_-’a aP kji

V(‘a—XF+raJ )= P —J-+ vg VJk. i=1. 2. 3.
i
'ali'+r;iva=
6X

Where p is the density. and v is the coefficient of kinematic
viscousity 'both. assumed constant. After some rearrangement. the

momentum equations can be written as

2
(1.1) (mute—1+ Uzi) u1 + .403)?" = - £§—+ v div (73'le ) - —U1
ax J5

ax2 1

g:

(1.2) V5(Ul;§T-+ 026:2) U2: J —pa:; + v div(v5'vU2)

 

2
.r 1 6 U26 ”3 103 1 a?

and the continuity equation

(1.4) infill) + 2«7592) = o
ax 6X2
where for a scalar function ¢(X1 X2)
=N_a£ _5‘£.B
axl axzB

and for a vector function 4’ = («91. 4’2. «93)

1 2
div o = 92T-+ 935- (vi = ¢1(X1. X2))
ax ax

Thus the gradient v and the divergence div. and later the Laplacian
are all two dimensional operators. This convention will be adopted
throughout the thesis to simplify notation.

We normalize the independent and dependent variables as follows:

1 J v 122.
X -axJ.U '5-11". P: 2 p
a
Substituting in equations (1.1) — (1.4) yields the dimensionless

momentum and continuity equations

3

(1.5) VAT-(£1153:- + u --a—)u1 + eu2 1

2
-JjJZP...d1 Jj .. 2.
6x1 v( vul) Jj' u

253—29112 - {T g + div(~/.Tvu2)

(1.6) «flulag—l + u

1 2
(1.7) (Hula-x6: + “26—63(2‘13 -. eulu3 = 4R 4- div (VJ-VH3) - 3- u:3
(1.8) diiju) = o
-a3 a?

where R = — is the Reynolds' number.
4p02 6X3

 

6
e = a/E is the curvature ratio. and now vfi'= 1 - exl.

Equations (1.5) - (1.8) are solved subject to the boundary conditions

(1.9) u = 0 1:1, 2, 3,

where F is the boundary of the unit disk
2 2
0—{x1+x2(l}

Note that when e = 0. equations (1.5) - (1.9) reduce to the equations
describing the motion in a straight pipe.

The velocity component u will henceforth be called the axial

3
velocity (or primary flow). and the transversal velocities ul. u2 will
be referred to as the secondary flow. Most of this work is devoted to

investigating efficient approximations of equations (1.5) - (1.9).

1.3 Literature review

In this section we give a brief account of some studies related to
our work. An extensive compilation of references can be found in the
review by Berger. Talbot. and Yao (1983).

We begin by giving an equivalent set of equations to (1.5) - (1.7)
which uses the Dean number D = 4Rv22; and the curvature ratio e as the

governing parameters. If in equation (1.5) - (1.7) we substitute

u
u3 = -§-and drop the prime. we obtain the following set of equations.

2
(1.10) Jflul-a-g—li» u2 522-) 111 + -;-u3 = — Jib-33+ div(\/.Ivu1) - Eul

J3 E + div(\/.Tvu

(111) «3(u15x_1+u2afiz)u2 =- 6x2 2)
(:2
(1.12) Jfiulaxl + u2 6x2) u:3 - eu1u3 = D + div (\[I- vu3) - JET-“3

(1 . 13) div(u~/J')

ll
0

where D = M is the Dean number as defined by McConalogue and
Srivastava (1968). Clearly equations (1.10) - (1.13) are completely
equivalent to equation (1.5) - (1.8) except for renaming the parameters
and rescaling the axial velocity “3'

Now if we set 6 = O in equations (1.10) - (1.13) we obtain a
conlnonly used set of approximte equations describing the motion in

"the limit of zero curvature".

_6__ 6 1 2 fl
(1.14) (u1+ u —) 111 + u = - + Au
2 6x2 2 3 6x1 1
_a_ .2. - - .92
(1.15) (111 (ix + u2 6x2) u2 .. 6x2 4- Au2
(1.16) (u1 33% 1121-) 113 -n+ Au3
1 “"2
(1.17) div u = 0
where A = i + i is the 2 dimensional Laplacian.
2 2
“"1 “"2

Equations (1.14) - (1.17) will henceforth be referred to as Dean’s

8

equations for the motion in the limit of zero curvature. The same
equations can be seen to be equivalent to the equations used by
McConalogue and Srivastava (1968). Greenspan (1973). Collins and Dennis
(1975). Dennis and Ng (1982) and Dennis (1980). except that the authors
of these works employed a polar coordinate system in the plane of the
cross section. and used a vorticity-stream function-axial velocity
formulation. The authors of the papers cited above used a variety of
numerical methods to solve Dean’s equations. and among them Collins and
Dennis (1975) seem to have achieved the most accurate set of results
covering the range 96 g D S 5000. The accuracy of their calculation
was later confirmed by the study of Dennis (1980) and that of Dennis
and Ng (1982).

Collins and Dennis (1975) also report very good agreement with
experimental results and results obtained by applying boundary layer
methods which they documented in their paper. .

We believe that there is enough evidence of the reliability of the
results of Collins and Dennis (1975). and therefore take their results
as a reference set of solutions against which we measure the accuracy
of our solutions. It should be noted that all the references mentioned
above deal with the flow in the limit of zero curvature. and that none
of them addresses the effect of curvature on the flow.

There are a few works that attempt to study the effect of
curvature. Truesdell and Adler (1970) used a finite difference method
to obtain solutions up to D = 3578 and e = .1. with a set of equations
first derived by Cuming (1952). The equations used by Truesdell and
Adler (1970) can be manipulated into equations (1.5) - (1.8) of this
work. The results of Truesdell and Adler (1970) however were not

sufficiently accurate beyond D = 1131. as the authors of that work

9

noted in their paper. Most of the works mentioned above used finite
difference methods which often requires a great number of mesh points
to achieve reasonable accuracy. which in turn puts tremendous bearing
on computer memory and execution time (Collins and Dennis (1975) used a
finite difference mesh with 2921 points and achieved the same accuracy
we did in this work with a finite element mesh of 489 nodes). Many of
the above authors also used a vorticity - stream function - axial
velocity formulation. and were confronted with the difficulty of
approximating the axial vorticity at the boundary of the region.
Finite difference approxinations also suffer the further difficulty
that for different cross sections. different formulations of the
problem must be used. For example. Collins and Dennis (1975) used a
polar mesh for a circular cross section and were forced to use the
equations of motion in polar coordinates. and later Collins and Dennis
(1976) had to convert the equations to rectangular form in order to
solve the flow problem in a triangular duct.

The switching of formulation often requires substantial
modifications of numerical techniques and computer code. We view the
problem of duct flow as essentially one problem regardless of the shape
of the cross section. and consequently believe that a unified approach
is needed to avoid redundancy. The situation becomes worse if one uses
a mesh that does not conform to the shape of the boundary. An example
of this is the work of Truesdell and Adler (1970) who used a
rectangular finite difference mesh to compute the flow in tubes of
circular and elliptic cross sections. and there they used a number of
fictitious boundary points to avoid the difficulty of approximating
derivatives near the boundary. Works on duct flows in cross sections

other than the circle include Cuming (1952). Truesdell and Adler

10

(1970). Joseph. Smith and Adler (1975). Collins and Dennis (1976).
1.4 Objectives and Organization

Our first goal in this work is to give an accurate mathematical
description of the problem. In chapter 2 we introduce appropriate
solution spaces and give an accurate variational formulation of
equations (1.5) - (1.8). We then consider the questions of existence
and uniqueness of the solution. The main objective of this work is to
build efficient numerical schemes to solve equations (1.5) - (1.8). and
since we believe that finite difference methods are rather cumbersome
and expensive for the reasons stated in section 1.3. we utilize finite
element methods to achieve our goal. Our solution method avoids all
the problems classically associated with finite difference methods. and
we believe it is far more efficient than methods used by earlier
investigators.

In Chapter 3 we describe a finite element scheme to reduce the
infinite dimensional problem to a finite dimensional nonlinear problem
whose solutions are described in Chapters 3 and 5. and there we also
give some details on the practical implementation of those schemes.
which we hope will cast some light on the advantages of finite element
approximations. and hence the computational efficiency of our solution.
In chapter 4 we present our error analysis. and give a bound on
Ilu - uhll. where u is the exact solution. uh is the approximate
solution. and II 9 II denotes the norm in 111(1))3

Numerical results and the physical discussion are given in Chapter
6. In Chapter 7. we apply a perturbation method to obtain solutions

for the flow in a tube of elliptic cross section. and use our scheme to

ll

generate a numerical solution of the same problem. This was done for
two reasons. On the one hand. it supplies further evidence that our
computation is accurate. and on the other hand it demonstrates the ease
with which finite elements are able to handle irregular boundaries. In
particular the same scheme and the same computer program can be used to
solve the problem for any shape of the cross section at the small
additional expense of generating a finite element mesh to discretise
the cross section.

From the physical point of view. it is our aim to show the effect
of curvature on the velocity field. and to show the significant
deviation of our solutions from those obtained in the limit of zero
curvature.

We also study the effect of curvature on the flow rate. and

discuss the adequacy of each of the parameters R and D to describe the

motion.

CHAPTER 2: EXISTENCE AND UNIQUENESS

2.1 Preliminaries and notation

In this chapter we use standard functional analytic notations.
For example if B is a Hilbert space. we denote the inner product by
(...). and distinguish the duality bracket between H and B” by the
superscripted bracket (...)*. thus for x e H and f e H“. (f.x)* = f(x).
The reason for this is that we do not always identify a Hilbert space
with its dual in the canonical way via the Riesz-Fisher theorem.
Throughout the chapter. 0 is the open unit disk. and all integrals are
taken over 0. We also use standard notations like 63(0) for Cu0
functions of compact support in 0. and 113(0) for real valued functions
on 0 which vanish in the sense of trace on the boundary of 0. and whose

generalized derivatives are also in L2(0). with norm defined by
2 2 an 2 du 2 1
Hull = Ilullo+ Ila-{Hy Hag-Ho iorueuocm

where II° II denotes the usual norm on L2(0). II° I ID will be used for
O

the norm in Lp(0). and to distinguish the norm on 113(0). we simply
denote the norm on the latter space by II- I I. We also use standard
notations like 003(0)3 to denote the space 63(0) x 03(0) x 03(0) and
similarly for 115(0)? Hé(0)3 etc. Product spaces are given the product
norm. Thus for u = (ul. u2. u3) e llama (vector functions will be

denoted by boldface letters)

MOD

2
Hull =

Mum’-

1:1

12

13

The continuity equation

6 iii 6 JJ‘
div(ui/j)=—(;1‘—l-—).+_(:_i2._)=o

suggests the introduction of the following space in which the above

equation is automatically satisfied.
Difinition 2.1.1:
For 0 S e < 1. define
V. = {u e 113(0)“: div (um = 0} ; 115(0)“

Observe that Vo

divergence-free vector fields which has been studied in detail in Temam

= {u e Hé(0)2: div u = 0) is the space of

(1984) and Girault and Raviart (1979). Note that Va is the kernel of
the linear operator Le: Hé(0)2 -» 19(0) defined by Len = div(m/.T). Note
also that Le is the composition div 0 Me where div is the divergence
operator and Me: lrl(1)(0)2 -b ll(l)(0)2 is the multiplication operator defined

by

g.

Mu
e

Lem 2.1.1:
V6 is a closed subspace of 06(0)? Thus V6 is a Hilbert space

with the norm II-II.

14

Proof: Observe the the divergence operator is continuous. and that Me

is a linear homeomorphism on llcl)(0)2 since J5. fl: 6 Cw(0). Thus L6 is
J

bounded and V6 = Ker Le [J
Definition 2.1.2:
2 2 .
Lo(0) = {p e L (0). [pdx = 0}
Learn 2.1.2:
The linear operator Le defined before lenlna 2.1.1 maps [1(1)(0)2 onto
2
LO(0).

Proof: Clearly for u e l-l(1)(0)2 we have

[div(u\/j)dx = - I {Tum ds = 0.
l"

l 2 1 2 2
since u 6. 110(0) . Thus Le maps l~l0(0) into L0(0).

By theorem 3.5 page 32 of Girault and Raviart (1979). the
divergence operator naps [1(1)(0)2 onto Lg(0). Hence given q 6. L(2)(0).
there exists v e 11(1)(0)2 such that div v = q. Set 11 = v/fi. Clearly
Leu = q B
Lem 2.1.3:

0° 2
V‘5 n 00(0) is dense in V6.

Proof: By theorem 3.8 page 36 of Girault and Raviart (1979).

15

V0 n C370)2 is dense in V0. The result now follows directly since ME

is a linear hemeomorphism on Hé(0)2. El
We now define the velocity solution space:

Definition 2.1.3:

x = ve x 113(0)

Observe that the continuity equation (1.8) does not involve the third
veloctiy component u:3 and thus X is the natural space for a solution u
of equations (1.5) - (1.8) to exist.

Since X is given the product norm. and 03(0) is dense in. 113(0).

the following lenlna follows directly from lenlnas 2.1.1 - 2.1.3.

Lenlna 2.1.4:

(a) X is a closed subspace of Hcl)(0)3. and is therefore a Hilbert space
(b) The operator Le maps il(l)(0):3 onto Lg(0). and X = ker Le“

(c) X fl 63(0)3 is dense in X. D

Remark: The two dimensional operator Le is extended in the obvious way

aiulm aiuzm

Leu - T... T . u = (111. u2. us) a l-l(1)(0)3

2.2: Variational Formulation of the Problem

To motivate the variational definition of the problem. let us

l6
assume that smooth solutions (u.p) exist for equations (1.5) - (1.8).
Multiplying equations (1.5) - (1.7) respectively by wl. w2. w3 e Hé(0).
integrating over 0. we obtain the following equation after integrating

the viscous terms by parts and adding the resulting equations.

3
2 1 a d
I vu ° vwwfi dx + e I 3(u1w1+ u3w3) dx + 2 IVE wJ(u1§-:-+ u -—-) ujdx

j=l 26x2

+ e [ (u3u3w1 - u1u3w3)dx = - [Vp ' wVfi dx + I F 0 w dx

0
where u = (ul. u2. u3). w (wl. w2. 3). F = [O J. and

4R
3

F o I =J:1FJ'J = 4RW3.

The above equation can be written more compactly in the form

ao(u.w) + B(u. u. w) = - (Vp. who + (F. w)0

. 3
where (.'.)0 denotes the inner product in L2(0)3. vu. vw = E vu ° vw
i=1

and
(2.1) ao(u. w) = [vu - vw JJ' dx + e2 [ (ET-(“1'1 + u3w3)dx

3
6 6
(2.2) B(u.v.w) = 2 f vfivj(u1 5;—-+ u2 5;;)vjdx + ef(u3V3w1 u1v3w3)dx
J=l 1
We now define the variational form of equations (1.5) - (1.8) and refer

to the problem by (P).

(P) Find u e x and p e L(2)(0) such that

17

ao(u.w) + B(u.u.w) = - (Vp. wVfi)“ + (F.w)* for every w e Hé(0)3.

In the above. (Vp. “(I)” is the duality bracket between [161(0):3
and Hé(0)3. thus Vp is viewed as a bounded functional on 113(0)3 defined
by (Vp. wVfi)” = - fp div (vvfi) dx. For the rest of this chapter

-1 3
F e HO (0) .

Proposition 2.2.1:
a0 is continuous. symmetric. bilinear. and coercive on B(l)(0)3.

Thus there exists a constant a > 0 such that

(2.3) acme) 2 a Ilul 12

Proof: The result follows immediately from the Poincare inequality and

the properties of V5. 0

Proposition 2.2.2:
B is well defined. trilinear and separately continuous on H5(0)3.

Thus the number B ) 0 defined by
B u.v.w . 1 3
(2.4) B = sup { u v ' . u. V. w e Hb(0) }

is finite.

Proof: The trilinearity of B is obvious once we prove that all the

integrals appearing in the definition of B are finite.

Some of the integrals defining B have the form

18

[fiwju dx(1$i$21$j$3)whereu v.chl)eH(0)

_.l.
iax1 ' i' j

By the Sobolev imbedding theorem. B(l)(0) is continuously imbedded in

L4(0). thus ui. w

6v
and since d—xl e L2(0). the generalized Holder inequality guarantees the
i

is L4(0).andl|u1||4 s c lluill. IIw_,II.1 s Cllell.

finiteness of the above integrals. and together with the imbedding

theorem and the boundedness of J5 gives

IIrwJu .6. a—"=1I dxsc Ilv II4 Ilu II4 II—lllosc llw II Ilu II IIv II
sCIIuII IIvII Ilwll.

In the above. the symbol C denotes different positive constants. The
other integrals involved in the definition of B have the form

[u 1v kadx. Using a similar argument. one can show that

IIungkIdx s c ||u|| IIvII IIvII.

The above inequalities imply that IB(u. VI. III S C Hall IIvII IIwII
and the proof is complete. 0
Now that we lmow that problem (P) is well defined. we give an

equivalent form of the problem. Consider the following problem
(Q) Find u e X such that

ao(u.w) + B(u. u. w) = (F. w)” for every w e X

19
Clearly if (u.p) is a solution of problem (P). then u is also a
solution of problem (Q) since div (u VJ.) = 0.

We now prove that the two problems are equivalent in the sense of

the fol lowing proposi tion.

Proposition 2.2.3:
Suppose u e X is a solution of problem (Q). Then there exists a
unique function p 5 [3(0) such that the pair (u.p) is a solution of

problem (P) .

Proof: Let us view Le as bounded linear operator from Bcl)(0)3 to
L(2)(0)*. where Lg(0)* is the dual of Lg(0) and the two spaces are
identified via the Riesz-Fischer theorem. Thus Le: l'l(1)(0)3 -b 13(0)” is
defined by

(Leu.q)* = I qLeu dx = I q div(u\/.I)dx. q 6 Lg(0)

By lelmiia 2.1.4.. L‘5 naps II(I)(0)3 onto L(2)(0)*. thus the range of Le

is closed. hence
(2.5) Range (LZ) = (ker 1.6)" = x”

* 2 ' 1 3 fl -1 3 + .
where Le: LO(0) 4 (110(0) ) = “0 (0) is the adjoint of Le and X is

the annihilator of x defined by x“ = {L e H31(n)3: l.(X) = 0}. Observe

that by definition

(8.31. n)” = new)" = Iq diva/fldx = - (vim/1')? q a L310)»: e 113(013

20

Now assume that u e X is a solution of (Q). Define an element
L e H61(0)3 by (.(w) = a0(u.w) + B(u.u.w) - an)“. w e 113(9)?

By assumption L(w) = 0 for w e X.thus L e X+ and by (2.5)
L 5. Range (LZ).

Thus there exists p e Lg(0) such that L = L:p i.e.

a0(u.w) + B(u.u.w) - (r.v)" = (sz.v)” = -(Vp.w\/.I)*
for every w e 115(0)?
We have proved the existence of the function p promised in the

statement. The uniqueness of p follows from the fact that Le maps

i-l(1)(l'2)3 onto Lgm) . D

2.3: Existence
In this section. we prove that problem (Q) has a solution.
Uniqueness is discussed in the next section. First we need the

following properties of B.

Proposition 2.3.1
(a) B(u.u.u) = 0 for every u e X.
(b) B(u.u.v) + B(v.v.u) = B(u—v.u—v.v) for every u e X. v e X.

Proof: Assume first that u e X n C3(0)3. In this case we have

3
a 6
B(u.u.u) = 2 I (In (u —+ —)u dx =
1:1 j 16x1 26x2 .1

21

3
2 “fl-u

l a a 2
_. __ + Ju __.)(u )dx = 1fu2d1v(m/j)dx =
2 1:1 16x1 26x2 1

.1.
2

IIMW

n w 3 n
Now for u e X. choose a sequence u e X n C062) such that u -> u
in X; This is possible by lemma 2.1.4. The trilinearity of B yields

the identity

B(u-un. u-un. u-up) = B(u.u.u) - B(up. up. up) + B(u. un-u. up)

+ B(ug-u. up.u) + B(up. u. up-u)

If in the above identity we let n 4Ȣ. use the fact that B(up.up. up) =

O and the separate continuity of B (proposition 2.2.2). we obtain

B(u.u.u) = O which proves (a). The proof of (b) is trivial.

Definition 2.3.1:

We define the auxilliary trilinear form

2.6 B =2 J3 dx
( ) 1(u' V W) j=l I "1(ulaxl+ “26x2)vj

Observe that
(2.7) B(u.v.w) = Bl(u.v.w) + e f(u3v3wl — ulv3w3)dx

The proof of the following proposition is similar to that of

proposition 2.3.1 and is therefore omitted.

Proposition 2.3.2:

For u. v. w e X

22

(a) Bl(u.v.v) -

l
O

(b) Bl(u.v.w)

-Bl(u.w.v) 0

Proposition 2.3.3:
If a sequence “n in X converges weakly to an element u e X. then

for every v e X we have

lim B(un.un.v) = B(u.u.v).
[H on

Proof t We prove that every subsequence of the numerical sequence
B(un. u“. v) has in turn a subsequence that converges to B(u.u.v). We
first prove the proposition for v e X n C;(0)3. By equation (2.7) and

proposition 2.3.2 (b). the following is true
n n n n n n n n
B(u .u .v): B1(u .u .v) + e f(u:3 u3 v1 - u1 u:3 v3)dx
n
=-Bl(un.nu.vu)+e[(3un 3 v1 -ulu3v3)dx

=_2 \U + dx-I- n n _ n n )dx
ng 1] uJ(ulax—l- “3632)” e[(u3u:'3v1 u1u3v3

Since un is weakly convergent in X. IIunII is bounded. thus by the

Sobolev compact imbedding theorem. un has a subsequence u“k which

converges strongly in L2(0)3.

l
4uJu1 in L (0) as

Thusu -»u inL2(0) ask-H”. henceuJ

“k ““14
J J “1

23

av
k -b 0°. and consequently flu?u1nka—vl-> {T u jui 5x1 in L1(0). since

1 1
and its derivatives are in I.” (0).

6x

v e C” 0(l'l) and hence v

J J

The above shows that

lim B(unk. unk.v) = -Bl(u.v.u) + e f(u3u3vl - u 111 3v

Rain

3"!"

= B1(u.u.v) + e f(u3u3v1 - ulu3v3)dx = B(u.u.v).

To prove the result for v e X. choose v1 e X n Cm 03(0) such that

IIv-vlll is "smll" and use the inequality

IB(un .u un.v) - B(u. u. v)I S IB(un .u n.v)- B(un. un.v1)l

+IB( un. uvn. 1) - B(u. u. v l)I+ [B(u.u.v1 ) - B(u.u.v)l D

Proposition 2.3.4: (apriori estimte)

A solution u of problem (Q) satisfies

a:

(2-8) Ilull S Ill-“ll la

where a is defined by (2.3) and HF“ = —— : w e X
llwll

Proof:

By assumption:

ao(u.w) + B(u.u.w) = (F.w)* for every w e X.

24

Choose w = u and recall that B(u.u.u) = O. we obtain

ao(u.u) = (F.u)*.

2
ll

Thus a Ilu S IIFII” IIuII and the result follows 0

Existence now follows from the theorem below. Its proof can be

found in Temam (1984) and Girault and Raviart (1979).

Theorem 2.3.5: (Existence)
Let H be a separable Hilbert space. Let ao be continuous.

symnetric. coercive bilinear form on H, and let B be a trilinear.

separately continuous form on H which satisfies
(a) B(u.u.u) = O for every u e H.

(b) lim B(un. u“. v) = B(u.u.v) for v e H and a sequence un converging
n-om

weakly to u in H.
Then for f 6 HT. the problem

ao(u.w) + B(u.u.w) = (f.w)* for every w e H

has a solution u in H. D
2.4: Uniqueness

x
Let the numbers a and B be defined by (2.3) and (2.4). and for F e X .

define

25

fl .
(2.9) IIFII” = sup {4%§;f%_ll: v e X. v # 0}
v

Theorem 2.4.1: (Uniqueness)
If
* 2
(2.10) IIFII (a/p
then the solution u of problem (Q) is unique.
Proof: Let u. v be solutions of problem (Q).
Then for every w e X
x
ao(u.w) + B(u.u.w) = (F.w)
x
a0(v.w) + B(v.v.w) = (F.w)
subtracting the last equation from the one above we have

ao(u - v. w) = B(v.v.w) - B(u.u.w)

Selecting w'= u - v in the above equation and using proposition 2.3.1.

we obtain
ao(u - v. u - v) = B(u - v. u - v. v)
By definition of a and B it follows that

2
aHu-flFsalh-vn Hfll

26

Since v is a solution. it satisfies the apriori estimate (2.8) thus

N
allu-v||2gEl-lfll IIu-vll2 i.e.

a

(a2 - burn") nu - m2 s

O

The hypothesis of the theorem now implies that u

u
5

Remark: In our problem. the functional F is the vector F = [O

on X by integration. thus
* 3
(my) = 2IFwdx=4RIw3dx
i=1 J J

Ilwdxl
and IIFII*=4Rsup{——3—:weX.wa£O

llwll

IIWII '
$Msup{—l-l-3—l-|2:wex.w#0 $41M?
W

0

Thus we have the following corollary which asserts the uniqueness of

the solution for "311311" Reynolds numbers.

2

 

Corollary: For R ( a . the solution of problem (Q) is unique.

Chapter 3: THE APPROXIMATE PROBLEM AND ITS SOLUTION

3.1: Motivation and Statement of the Problem.

In this chapter we introduce finite element (F.E.) analogues of
the continuous problems (P) and (Q). Our F.E. formulation resembles
the well-known Taylor-Hood approximation of the Navier-stokes
equations. References on F.E. approximation of the Navier Stokes
equations include the works of Hood and Taylor (1973). Girault and
Raviart (1979). Glowinski and Pironneau (1979). Bercovier and Pironneau
(1979). Le Tallec (1980). and Glowinski (1984).

We use Co piecewise biquadratic polynomials to approximate the
velocity components. and Co piecewise bilinear polynomials to
approximate the pressure. The reference elements and their respective
basis functions are shown in Figure 2.

Let the domain 0 be approximated by a union Oh of elementary
rectangles. Observe that the unit disk 0 cannot be covered exactly by
elementary rectangles. and thus Oh it 0. In practice however. one uses
isoparametric elements to accomodate the curved boundary. and since in
this thesis we do not study the effect of using isoparametric elements.
we shall disregard the distinction between 0 and Oh. and assume that
O = flh'

As is usually done in F 0 E - approximations. we are thinking of a

family ( of tessilations of 0. indexed by the parameter h which is

7131:
the maximum diameter of an element in Th. We also assume that h -) O.
The precise assumptions on the family {Th}h will be stated in the next
chapter.

Let Q1 and Q2 denote respectively the spaces of bilinear and

biquadratic polynomials in two variables. and define the finite element

27

28

v4
2
I

x(x-1)y(y-1)/4

. - - -1 2
(1'1) N2 (x 1)(x+1)y(y )/

 

x(x+1)y(y-1)/4

\1
$!
0‘
ULJ
z
u
I

-x(x+l)(y-l)(y+l)/2

2
p
I

x(x+1)y(y+l)/4

L

 

 

 

 

 

X N . -(x-1)(x+l)y(y+1)/2
8 9 a 6
N7 = x(x-1)y(y+1)/4
1 N8 = -x(x-1)(y-l)(y+l)/2
(-1,-1) 1L u - (x-1)(x+1)(y-1)(y+1)
1. 2 3 9

The reference element for the velocity and its basis functions.

 

 

 

 

 

Y
(1.1)

4 ._A. 3
M1 - (x-l)(y-l)/4
M2 = -(x+1)(y-1)/4
x M3 - (x+l)(y+l)/4
M4 - -(x-1)(y+l)/4

1‘: 2

(-1.-1)

 

The reference element for the pressure and its basis functions.

The reference elements

Figure 2

29

spaces HOh and 11h as follows:

(3.1) "(1”! = {uh e Co(fi): uh

 

e4';Q2.!.1h|r=0}

(3.2) 11., = {qh e c°(s‘i)= qh e 01. Iqh dx = 0}

 

In the above. e denotes a typical element. 1‘ is the boundary of 0. and
the integral in (3.2) is taken over 0.

Observe that "Oh is a subspace of Hcl)(0) and that "h is a subspace

of 13(0) (in fact is a subspace of 111(0)). The spaces (“(1)192 and

1 3

equation div (m5) = O is equivalent to the condition

are defined in the obvious way. Observe that the continuity

(3.3) I q div (us/fidx = O for every q e L§(0)

and in the approxinate problem we replace (3.3) by the discrete

condition that the approximate velocity uh satisfies

(3.4) fqh div (uhV/j) dx = O for every qh e "h

We now define the approximate solution space for the velocity as

follows:

(3.5) veh = (uh e (113,192: I qhdiv(uh\/.T)dx = O for every qh 6 Mb)

(3.6) . xh = v‘5h x “(luv

30

Remrk: veh is not a subspace of V6. and consequently Xh is not a
subspace of X.
We now state the approximte versions of problems (P) and (Q) of

chapter 2 .
Problem (Ph)

Find u.h e (H311

)3 and ph e uh such that
30(“h"h) * B(uh'uh’ h) = ‘(Vph' 'h‘mo + (F' b)”
for every wh 6 (H3193
I qhdiv (tuba/I) dx = o for every qh e uh
Problem (oh)
Find uh e xh such that

ao(uh. h) + B(uh.uh.wh) = (F.wh)*. for every 'h e Xh.

As in the continuous case. we show the equivalence of (Ph). and (Oh).
To achieve this we need a discrete analogue of lenlna 2.1.2 which

guarantees the existence and uniqueness of the pressure.

3.2: Brezzi-type condition

The proof of the following technical lemna can be found in

Bercovier and Pironneau (1979).

31

Lemma 3.2.1:

If every element e 5 Th has at least one vertex which is not on

the boundary of 0. then there exists a constant C > O.

 

1 independent of
h such that
(VPh- 'h)o
1 2
Ilwhl lo
for all ph 5 lb. 0
Lemma 3.2.2:

If the assumption of lemma 3.2.1 holds. then for 6 "small enough".

there exists a constant C2. independent of h such that for all ph e Uh

( . w 45)
(3.7) sup “a“ h 0 = “h e (H3192. h 3‘ o 2 c2 llvphllo
llwhllo

 

Proof: Recall that VGD= l — exl. lel S 1. thus

(”h' 'h‘mo = (Vph’ 'h)0 _ e (”11' "1'h’0
Huh”. _ Ilwhllo llwhllo

 

 

 

(vph. uh). I(vph. wol 2 (wk, uh). llvphl Io llxlwhl lo 2
... e "' e
Hth I0 :1th lo IIth I. ”th I.

 

(VPh' 'h)o

- . llvrhllo
Huh”.

32

Thus by lemma 3.2.1 we have.

 

 

(VPh' 'h‘mo 1 2
sup : wh e (HOh) . h x 0 2
llwhllo
( . V)
SUP Vph ho: whe(H(l)h)2, hi0 ‘5 ”VPhHOZ
llwhllo
(Cl - E.) HvPhHO
Now take e < Cl' and set Cé = Cl - e. D

3.3: The Equivalence of (Ph) and (Qh)
Lemma 3.2.2 now implies the equivalence of (Pb) and (Qh):

Define Leh: (Héh)3 at]: by (Lehwh.ph)” = I phdiv (whv5)dx.'
Then the adjoint Lzh: uh -» ((Héhfl" is defined by

(Lthh' h)” = ‘ (”ph' 'h‘mo

x 1 3
Condition 3.7 implies that Leh is one-to-one. and hence Lehmaps (th)
onto I:. Now the proof of the following equivalence theorem is exactly

like that of proposition 2.2.3.

Theorem 3.3.1:
Problems (Pb) and (Qh) are equivalent in the following sense: If

(uh. ph) e Xh x lh is a solution of (Pb). then uh is also a solution of

33
(Oh). Conversely if uh is a solution of (Qh). then there exists a

unique phe "h such that the pair (uh. ph) is a solution of (P h.) D

Remrks :

(a) The space Xh is the kernel of the linear operator Leh' and thus if
N is the number of interior mesh points for the velocities. and M(< N)
is the number of mesh points for the pressure. then lh has dimension

.(0 1h)3 has dimension 3N and Xh has dimension 3N - M + 1.

1
(b) Let {NJ}J=1. {119:1 1 be the basis functions for HOh and Eh

respectively. Define a mtrix D = (D1.Dz) by

(3-8) D = I NJ—dex

__i_“
(3.9) 132: NJ 5x2 fidx

lSiSM-l. ISJSN
Then condition (3.7) implies that the matrix D has full rank (M - 1)

3.4: The Solution of the Approximate Problem.

We turn now to the question of existence and uniqueness of the
approximate problem (Qh). Observe that the restriction of the
trilinear form B to the space Xh does not have the properties given in
proposition 2.3.1. because the space xh is not a subspace of X. and

thus the elements of Xh do not satisfy the continuity equation exactly.

34

This makes problem (Qh) more difficult to study since for example the
apriori estimte (2.8) is no longer valid for the solution of (Qh).
However. at the additional cost of restricting the size of the forcing
term IIFII” more than we need for uniqueness (Condition (2.10) in
theorem 2.4.1). we can show that problem (Qh) has a unique solution uh
which satisfies the same apriori estimate as the exact solution u.

The following condition is a standing assumption in this chapter

and the following one:

2
1‘ a
(3.10) IIFII < 33

where a and B are defined by equation (2.3) and (2.4).

2
Definition: For IIFII" < 3‘5. let

 

(3.11) r = (a-«eZ-ebnrn") /2p

Observe that r is the snaller zero of the quadratic equation

Br2 - ar +IIFII” = O. and that

(3.12) a - Br ) O

(3.13) a - 2Br > 0
Bl IFII"

(3.14) = ———2- ( 1
(ct-Br)

Observe also that when condition (3.10) is satisfied. the exact

35

solution u of the continous problem in unique. and by (2.8) it follows

 

that
HHP 2
(am) IhHs a =r-%}<r

Thus the exact solution lies in the ball
m={ue%mf=udlso

This suggests that we look for an approximate solution in Dr’ and

indeed there is a unique solution of (Qh) in Dr'
We begin with the following

Lemma 3.4.13

For uh e Xh. Iluhll S r. and for any L e XE. the problem
x
ab(vh. h) + B(uh’vh"h) = (L. h) for every 'h e Xh
has a unique solution vh.

Proof: The result follows from the Lax- Milgram theorem. since for all

vh e Xh we have

aoh’h' h) + B(uh'vhovh) 2 (“’5' luhl I) I lvhl '2 2 (a—Br) ”th '2

and a - Br ) O by (3.12). D

The above lemma allows us to make the following definition. Define a

36

mapping ¢= Dr diXh by ¢(uh) = vh. where vh is the unique solution of

(3.16) ao(vh"h) + B(“h'vh"h) = (F"h)* for every 'h e Xh

Lemma 3.4.2:

The mapping ¢ defined above maps Dr into Dr'

Proof: If we choose 'h = vh in (3.16) we obtain

:1th l2(a - m s ao(vh-vh) + B(u.,mhwh) = Mb)" .<. IIFI I“ Ilvhll.

N
IIFII

mus||vhllgm=r D

Remark: Observe that every solution of problem (Qh) is a fixed point
of 0.

Proposition 3.4.3:

¢ is a contraction on Dr'
Proof: For uh 5 Dr’ define Auh e 2(Xh.X;) by
*.
(Auh(vh). h) = 0(vh"h) + B(uh‘vh"h)'

'ao("h' h) + B(uh'vh"h)|:

 

*
IlAuh(vh)|| =sup IIvII wheXh.wh#0 2
h

Iao(vh.vh) + B(uhw .vhu a(vh. h) - IB(uhwhwhH

llvhll llvhll

 

 

37

m-pnqninmnzia-m)nmn

Since a - Br ) O (by 3.12). Auh is bounded away from zero. and is
therefore one-to-one.
Lem 3.4.1 states that Auh is also onto. The open mapping

theorem implies that Auh has a bounded inverse Au}:l and that

l
a-Br

 

min Ingflls

Observe that now by definition ¢(uh) = AuglF. We now show that O is a

contraction: let url‘. u: be in Dr’ and let

A1 = Au; (i = 1.2).

We have A-1 - A-1 A-}.

-1
2 1 = A 2(A1 - A thus by (3.17)

2)

l l

- -1
(3'18) IlAz-AIIIS—EIIA1_A2”
(“*BP)

1 2
It is easy to check that ”Al - A2” S B Iluh - uhll.

N.. Ilocufii - ¢<u;>ll = Ila-gr - A'irll s IIA'; - A'ill llFll”$
BHHP 2 1 2 1

----§-||uh - uhll = kIluh - uhll.

(a-Br)

where k < 1 is defined by (3.14) 0

Theorem 3.4.4: Problem (g9 has a unique solution uh in Dr'

annuqur—1

38

Proof: The result follows directly from proposition 3.4.3 and the

remark preceeding it. D

3.5: A Fixed Point Scheme.
Proposition 3.4.3 provides the following practical scheme for
finding the solution uh of (Oh). Choose an initial approximation

ug e Xi n Dr' and define a sequence (ufi) (n 2 1). converging to uh by

“:+1 = ¢(“:) 1-0- given “2' u:+l is the solution of the linear problem
n+1
uh exh-

n+1 l

(3.18) ao(uh . h) + B(ufimfi+ . h) = (F. h)” for every 'h e Xh.

Again by an argument similar to the proof of proposition 2.2.3 it is
easy to show that lemma 3.2.2 can be used to show that problem (3.18)

is equivalent to the following problem.

Given u: e Xh n Dr' find ufi+1 e (H0h)3‘ pfifl e Uh such that

n+1 n+1

(3.19) ao(u:+1.wh) + B(uh’uh , h) = -(Vph "hy530 + (F.wh)*

1 3
for every 'h e (HOh)

(3.20) I qhdiv (“2”l

Vfi)dx = 0 for every qh e Eh
Observe that the solution of (3.19). (3.20) can be obtained by solving
a linear system of algebraic equations which is constructed as follows:

let (N )N be the basis functions for H; and let {M )M-1 be the
1 i=1 h J 1:].

basis functions for Uh. then the space (H311)3 has the basis

”-21

 

39

(3.21) ((N1.O.O). (0.N1.0). (0.0.Ni)}: 1

Clearly. equations (3.19) is satisfied for all wh e (HOh)3 if and only
if it is satisfied for all the basis functions (3.21). The same

applies to equation (3.20) and the basis functions M Thus if we

1'
substitute the basis functions (3.21) in (3.19). we obtain 3 sets of
linear systems each containing N equations. Simiarly (3.20) gives rise
to a system of 11-1 equations. The combined system of 3N + M - 1
equations can then be solved for the nodal values of the velocities and
the pressure. The matrix of the system described above will be
referred to as the grand fluid mtrix. As an illustration. let us

describe explicity the equations that correspond to using wh = (0.0.111)

in equations (3.19).

Let ur'l = (“31:1 . “12:1. ugl). where for example
n+1_g N dltUn+l-( )T
“3h ‘ 1:1 “33 J' a“ e 3 ‘ “31' ° ° " “3N

be the nodal values of u3;1.'1he substitution wh = (0.0.111).

1 g i g N. yields the mtrix equation

(3.22) A3+1 Un+1 = c

h

where G is the column vector whose it component is (F.(O.0.N1))*. and

the (i.j) entry of A3” is

[VJ-vNi-VN

J 1 ax2

2 an 6N
+ (:75- ' “Elwin; + ‘5 Ni”? 51 + “g 4) dx'

40

Similar equations can be derived by substituting (N1.0.0). (0,111.0) in
(3.19) and NJ in (3.20). Observe that equation (3.22) involves the
nodal values of “13:1 only. and thus the grand fluid matrix can be
rearranged into a block diagonal mtrix. Observe also that this is no
coincidence. indeed the triliniear form B was defined in such a way to
guarantee this nice feature of the grand fluid matrix. Other
definitions of B are possible. and for some of these definitions. the
existence proof given in chapter 2 can be greatly simplified. however.
such definitions lead to more dense matrices than the one described
above. The reader is referred to Taylor and Hughes (1981) for more

details on the computer implementation of F.E. method in flow problems.
3.6: The discrete Stokes Problem.
The continuous stokes problem associated with the flow under
investigation has the form
(SP) Find u e li(1)((l)3 and p e L3H!) such that
a0(u.w) = -(Vp.w\’j)* + (F.w)* for every w 5 113(0)?

div(m/.T) = 0

Again problem (SP) can be shown to be equivalent to the problem
($) Find u e X such that

a°(u.w) = (F.w)* for every w e X.

Problem (SQ) has a unique solution by the Lax-Milgram theorem.

The discrete form of problems (SP) and (SQ) will have great importance

41

in clnpters 4 and 5. and hence deserve the brief description we give

below.

If we use the same F.E. approximations of u. p as was done in
section 3.1. one obtains the discrete versions of problems (SP) and
(s0):

(SPh) Find uh e (35193 and ph 45. uh such that

ao(u.h. h) = - (vph. thO + (F.wh)* for every wh 5 (H3193.

I qhd1V('h\/.T)dx = 0 for every qh e "11'

(80h) Find uh e Xh such that

}

x
ao(uh.wh) = (F.wh) for every wh 5 Kb

Equations (SPh) lead to the following system of linear equations for

the nodal values [11.02.U3 and P of ulh' u2h. “3h and (3h respectively.

 

 

 

 

 

 

’A o (1)1)T o‘ ' u1 ‘ ' c:1 ‘
1 “r
(3.23) o A (112) o 02 = (:2
' 131 132 o o P o
_o o .o A _ 03 . . G3 .

 

where the ith components of Cl. G2. G3 are given by

N
(cl)1 (F. (Ni.0.0)) .

N
(62)1 (F. (0.N..0))

42
)1
(c3)! = (F. (0.0.N1)

D1 and D2 are defined by (3.8) and (3.9). and the (i.J) extries of A
and A1 are

2
6
A13 = [(VN1 ° VNJ J3"?—

fiNiNJNX'

1

Aij

= [in1 - VNJJI dx.

Observe that D = (D1. D2) has full rank. and that A and A1 are
syuInetric. positive definite. and by numbering the mesh nodes
appropriately A and A1 are banded with the band width considerably
smaller then the size of these mtrices. The practicalities involved

in solving (3.23) are outlined in appendix B.

CHAPTER 4: ERROR ANALYSIS
In this chapter we give an error bound on Iluruhll. where
throughout the chapter. u will denote the exact solution. and uh will
denote the solution of the approximate problems (Ph) and (Qh)'
Thus u e X satisfies the equations
(4.1) ao(u.w) + B(u.u.w) = (p. div(w\/J'))o + arm)”
for every w e Hé(0)3.
(4.2) ao(u.w) + B(u.u.w) = (F.w)* for every w e X.
and uh e Xh satisfies the equations
(43) a1w1+B< )=-( mun)"
° ouh'h uh’uh’h wh‘h o ’h
1 3
for every 'h e (th)

(4.4) ao(uh. h) + B(uh'uh"h) = (F.wh)* for every 'h e Xh.

2..
”I
that satisfies (4.1) and (4.2). and uh is the unique solution of (4.3)

It will be assumed that IIFII” < thus u.is the unique element of X

and (4.4) that satisfies

(4-5) lluhll S r.

2
Observe that under the assumption that IIFII“ < g5; the exact solution
43

44

u also satisfies

(4.6) Ilull S r.

We also assume that the hypotheses of lemma 3.2.2 are satisfied. thus 6
is small enough to allow condition (3.7) to hold. We quote condition
(3.7) here for the ease of reference:

There exists a constant Cé ) 0. independent of h such that

(”911' 'h‘mo, 1 3

(4'7) 3UP '1'! 5 (H011) - 1') fl 0 2 02 I'Vphllo
Huh”.

for all ph e Mb.

In the sequel we will continue using the letter C to denote (possibly)
different positive constants. however when a specific condition such as
(4.7) is used. the constant C2 will be used to call the reader‘s
attention to the specific condition used. we assume that our family
(Th)h of tessilations of 0 is regular in the sense of Ciarlet (1978).

thus if he denotes the diameter of element e e T . and pe denoted the

diameter of the largest circle that can be inscribed in e. then

(i) There exists a constant a ) 0 such that for every e e Th and

every Th

he
3“"

(ii) The quantity h = max { he : e 5 Th} approaches zero.

45

Besides the regularity of the {Th} . we also assume that {Th)h
h

satisfies the inverse assumption:

(iii) There exists a constant v > 0 such that for every e e Th and

every Th

:3"?!
1A
:2

Observe that conditions (i) - (iii) are by no means restrictive in
practice. We begin by stating the following lenma whose proof can be
found in Ciarlet (1978)

Len-11a 4.1.

Under assumptions (i) - (iii). there exists a constant C3 > 0.

independent of h such that

Il'hl '0

(4.8) inf - wh 6 (H392. wh ,. o 2 can :1

Huh”.

Lem 4.2

For 11 e X. there exists a constant C > 0. independent of h such

that
(4.9) . inf {IIu-vhll:vh e Xh} S
CianIu-v Il+lIIu-v II =v e111?)
h h h o h H0h

. 1 3
Proof- Fix vh e (HOh) . and let (2h. wh) e th "h be the solution of

the discrete stokes problem

46

(4.10) ao(zh. h) = —(wh.wh\/j)o + ao(v .wh) for every wh e (HOh)3
(4.11) (vqh.zh\/J-)o = O for every qh e “h

The boundedness of a

0 yields

o‘vh' 1h"):

a
(4.12) IIvh-Zhll 2 C sup { II II w 6. H3603. w i O} 2
w

 

a(v - .w)
Csup{0 h Zl1h:whe(l-l‘l)l,1)3. hi0}:
llvhll

(w.w\/.T)
Csup{ h h O:whe(ll(1)h)3. hi0}
Ilwhll

where the last equality in (4.12) is by (4.10). (4.12) yields

(W-V‘Ij)
h 1‘ ozwhe(H(l)h)2. h¢0}2

 

(4.13) Ilvh-zhll 2 c sup{
llwhll

 

(W . w ff) H" II
c sup { h h 0 : wh 5 (11392} inf {——h—g : 'h e (115192}
. Ilwhllo Ilwhll

Now (4.13). (4.7) and (4.8) give
(4.14) IIvh-zhll 2 c czc3 Ilwhllo . h

We rename C C2 C3 to C. Using (4.14). the ellipticity of a0. and

(4.10) respectively we have

47

llv - 112 a(v-.v-1
(4.15) hIIvh-ZhHSC h 2" go 0 hzh hzh =
llwhllo Ilvrrhllo

 

 

C (W - (VhTZhh/jh)

 

Ilvrhllo
(4.15). (4.11). and the fact that u e X imply

(W - (vb-u)‘/j)0

 

hIIv - II S C
h zh
||WhHO

Thus

(W 0 (VII-uh/jh) .

 

(4.16) hIIvh-zhll S C sup { "h e uh} S

”th Io

 

{(' . (vb-u)‘/‘T)O
C sup :

w e L2(D)2} = C Il(vh-U)V5]|o S C I'Vh'ullo
Ilvllo

where the last inequality in (4.16) is by the boundedness of ¢3L By

the triangle inequality and (4.16) we have
(3
Ilu-zhll s IIu-vhll + ”vb-2.,” s llu—vhll + ..- llu-vhllo
and this implies (4.9) and completes the proof. 0
Theorem 4.3:

Under the assumptions stated before lemma (4.1). there exists a

constant C ) 0 such that

48

llu-uhll S Ch2

Proof:

Fix vh e xh. and let wh = uh - vh.
Consider the expression

E = ao(wh.wh) + B(wh.uh.wh)
By (2.3). (2.4) and (4.5) we have
(..m t 2 (a - Br)llwh| :2
By (4.4) and the definition of B we have
(4.18) 1: = ao(u.h.wh) — a0(vh.wh) + B(wh.uh,wh) =
-B(uh,u.h.wh) + (F.wh)" - ao(vh.wh) + B(wh.uh.wh) =

(F.wh)" - B(vh.uh.wh) - ao(vh.wh)
(4.18) and (4.1) yield
(4.19) = -(p.div(whY53)o + ao(u.wh) + B(u.u.wh)

‘ 80(vh’wh) "' B("h"’h"'h)

49

Using the fact that wh e Xh. a simple mnipulation of (4.19) gives
(4.20) E = -(p—qh.div(whvfi))o + a0(u-vh. h) + B(u—vh.uh. h)

where qh e "h is arbitrary. It now follows from (4.20). (4.5). (4.6)

and the triangle inequality that
(4'21) E S Cll'hll Hp—thlo + Cll'hll Hu-vh” + Brll'hll [nu-V11” +

llu-uhll] .<. c llwhll IIp-thlomllwhll IIu-vhll + 2n llwhll IIu-vhll
2
well...”

Now (4.17) and (4.21) yield

(«z-2m ”th12 s c Ilwhll [llp-thlo + llu—vhlll i.e.
(4.22) llwhll 4 152m [llp-thlo + Ilu—thIJ

Observe that a - 2Br > 0 by (3.13). Now the triange inequality and

(4.22) give

IIu-uhll s llu—vhll + IIvh-uhll s cmp-qhno + IIu-thIJ

for every vh e Xh. and every qh 6 11h.
Thus

50

(4.23) IIu-uhll s c 1nf{||p-thl=<1h «2 uh} + c inf {IIu-vhll= vh e xh}

(4.23) and (4.9) give
(424) Ilu-uhll s c inf {IIp-thlo= qh e 15,}
1 , 1 3
+ C inf {I'll-Vb” 4' Kl Iu-vhl '0' vh 5 (Hob) }

Now the result follows from (4.24) if we choose for vh the interpolant
1 3

Uhu of u in (11011) and for qh the interpolant th in uh. By the

general interpolation theory in Sobolev spaces (see e.g. Ciarlet

(1978). Oden and Carey (1983)) we have

2 3 2
IIu-Hhull == 0(h ). Ilu - Hhullo = 0(h ). IIp-Uhpllo = 0(h ). and the

proof is complete. 0

Remark: The results quoted above require that u 15 [13(0), p e “2(0).
Although our results in chapter '2 do not include a study of the
regularity of solution (u.p). we believe that the solution is actually
in COUZ). This is intuitively obvious since the domin D is bounded.

80 e C”. and the data (boundary and forcing) are all C”.

CHAPTER 5: LEAST SQUARES APPROACH

5.1: Introduction

In section 3.5 we described an iterative scheme to solve the
nonlinear problems (Ph) and (Oh). Each step of the iteration required
the solution of a large system of linear equations. but unfortunately
the coefficient natrix used to solve for ufi+l was dependent on the
previous approximtion n: Thus at each step of the fixed point scheme
(3.19) one must form and solve a different set of linear equations.
This turned out to be rather expensive since we were interested in
computing the solution over a wide range of values of R and e.

In this chapter. we give an alternative scheme for solving the
nonlinear equations (Ph) and (Qh). The method described here is a
straightforward extension of the techniques described in Glowinski
(1984). (Mr discussion will be brief since this method of solution is

well described in the above cited reference.

5.2: A least squares formulation

The assumptions and notations of chapter 3 will be maintained
here. Probelm (Qh) can be converted into a minimization problem as
follows:
For a fixed uh e Xh. define an element Eh e Xh to be the solution of
the discrete Stokes problem

Eh e Xh is such that for every "h e X'h we have

(5.1) ao(fh.nh) = ao(uh-nh) + B(uh-uh-nh) - mun)”

Problem (5.1) will be referred to as the state equation. and Eh will be

51

52

termed the state variable that corresponds to “h' Observe that the
right hand side of (5.1) is a linear functional in rib. and thus (5.1)
is a perfectly well-defined (discrete) Stokes problem whose solution Eh
exists and is unique by the elliplicity of a0. We now define a
functional Jh: Xh -b Xh by

(5.2) Jh(uh) =% 0(Eh' Eh)

where Eh = Eh(uh) is the solution of the state equation (5.1). The
proof of the following proposition is obvious if one assumes that
problems (Ph) and (Qh) have solutions which is the case for example if
IIFI I” < 112/413.

Proposition 5.2. 1:

An element u.h e Xh is a (global) minimizer of Jh if and only if uh

is a solution of problems (Pb) and (Qh)' D

As a corollary of the above proposition. Jh has a unique global

minimizer uh 6 Dr which is also the unique solution of (Pb) and (Qh) in

2
Dr ( here we again assume that IIFII” ( 5415 ). Now instead of solving

(Oh). we try to solve the following minimization problem:

(uh) Find uh e Xh such that
Jhwh) S Jh(vh) for every vh e Xh.

5.3: A conjugate gradient scheme.

Following Glowinski (1984). we use the following conjugate

53
gradient scheme to solve the minimization problem (Mh). Let ug e Xh be

an initial approximation of uh. and let 23 5 X5 be the solution of the

"linear variational equation"

(5.3) ab(z:.nh) = < Jg(u:).nh > for every "h e Xh.

o o n+1 m+1 n+1
Set 'h - zh. For m 2 0. compute uh . zh . h as follows

Step 1 (Descent)

Compute

(5.4) in = Argmin who; - Mfi). x e IR}
and set

(5.5) u:+1 =u:-7\m w:

Step 2 (New descent direction)

Define z:+l to be the solution of the linear variational equation

m+1

(5.6) ao(zh . nh) = < Jg(u:+l). "h > for every "h e Xh

 

. <4“- 4+1- 4)
(5.7) Set em” = 0 m m

30(211' 2h)
(5.8) and wig.” = 2T1 + «m1 w:
Replace m by m+1 and go to (5.4). a

In equations (5.3) and (5.6). Jh(vh) denotes the Frechet derivative of
Jh(vh + t"1.) ‘ Jh(vh)

 

Jh at 'h e Xh defined by < Jh(vh)' "h > = limo t

54

(1')h e Xh). It can be easily shown that

(5 9) < Jh(vh).nh > = ao(Eh.nh) + B(nh.v .fh) + B(vh.nh.§h)
where Eh = §h(vh) is the state variable corresponding to vh.
Observe that the linear variational equation (5.6) is a discrete

m+1

Stokes problem for the unknown function zh Observe that computing

2:” actually requires the solutions of two Stokes problems: first of
all. one must compute the state variable Eh(u.h) in order to use (5.9)
to compute the right hand side of (5.6). and then solve the Stokes
problem (5.6) for 2:“ once ( J};(uhm). "n > is known. The one
dimensional minimization problem (5.4) deserves a brief conment since
it is in fact the most expensive step in the above scheme. Given “E'

w:. define a real valued function

g(x) = Jh(u: - hwfi) A e m
and to find the minimizer Am in (5.4). we solve equation
(5.10) g'(x) = 0
Observe that
(5.11) g'(x) = - < J'(u: - Awfi). wh >

In our implementation of the above scheme. we used the modified false
position method to solve equation (5.10). Observe that each evaluation
of g'()\) (through (5.11)) requires the solution of one Stokes problem

to find the state variable Eh corresponding to u: - Awg. which is

55

needed in evaluating the right hand side of (5.9). Consequently the

solutions of several Stokes problems must be computed in order to find
Am“

It must be noted that although the conjugate gradient scheme
requires the solution of many Stokes problems in each iteration. it is
actually a very efficient scheme since only the forcing term needs to
be computed each time. but the bilinear form ao(. . .) is the same for
all the the Stokes problems needed to implement the above scheme. This
amounts to computing the matrix in equations (3.23) only once
throughout the entire computation. and once this has been done. solving
the Stokes problem only requires a back substitution to find the nodal
values of the velocity and pressure.

Solving the matrix equation (3.23) associated with ths discrete

Stokes problem is outlined in Appendix B.

CHAPTER 6: NUMERICAL RESULTS

6.1 Computational procedure

As we mentioned in section 1.4. the main objective of this work is
to find accurate approximations of equations (1.5) - (1.8). To achieve
this goal. we implemented the schemes described in chapters 3 and 5 to
compute the solution uh of problems (Pb) and (Nb) respectively. Since
an initial approximation is needed in both cases. the approximate
solution “h was computed for a sequence of values of R (keeping 5
fixed). starting at the initial value R = 25 and increments AR = 25.
The initial approximation for R = 25 was chosen to be zero. and after
the solution for a certain value of R has converged. R was incremented
and the solution for the previous value of R was used as an initial
approximation for the next value of R. This procedure was repeated
until the final desired value of R was reached.

The above procedure was carried out for e = .1. .2. . . .. .5 and
for values of R ranging from R = 25 to R = 1118.03.

The criterion

lama) - u“ ml _
(6.1) max 5“ 11‘ < 10 3

lug‘hii) I

 

was used as an indication that the sequence uhn has converged to the
solution uh.

In (6.1) the maximum is taken over 1 S j S 3 (the velocity
components) and over all the nodes i of the mesh used. Observe that

(6.1) implies that

56

57

ll'“- “II
“h “h w<10-3

(6.2)
Hugl I.

Note however that (6.1) is a much more stringent condition than
(6.2). The same computer programs were used to solve Dean’s equations
(1.13) - (1.17) for the range 96 S D S 3500. The following slightly
modified version of the fixed point scheme (see section 3.5) was found
to accelerate the rate of convergence considerably

Given an initial approximation ug. define

D

= ¢(uh)I

MT.—

.3

or...

n+1 n

uh = NJ: + (l-hhlh

where O S A < l is an averaging parameter whose optimal value is
determined by numerical ewerimentation. All the computations were

carried out in double precision on a VAX - 11/750 computer.

6.2 Some comparisons

First of all it must be mentioned that the two schemes gave
identical results for all the cases tested. but as expected. the least
squares scheme proved to be much faster than the fixed point scheme of
Chapter 3.

This must not be surprising since the fixed point scheme requires
a huge number of matrix factorizations. The relative speed of the two
schemes depends on the mesh size. for example. to compute the solutions

of Dean’s equations for 96 S D S 1000 (AD = 100) using MESH 24 (see

58

Figure 3). the fixed point scheme requires 39:04 minutes of CPU time.
while the least squares scheme requires 57:45 minutes. Observe that in
this case the fixed point scheme is faster than least-squares.
However. to compute the flow at D = 2100 using MESH 96. the
least-squares scheme requires 1:16:04 hours vs. 5:45:57 hours for the
fixed point scheme.

We now turn to practical accuracy considerations.

In order to determine an adequate mesh size that gives good
accuracy. three different mesh sizes were used to solve equations
(1.14) - (1.17). Figure 3 shows the three meshes used. Table 1 shows
the values of the maximum axial velocity “3m at different Dean
numbers. while table 2 shows the values of the axial velocity u3(0) at

the center of the cross section

59

4

' . MESH 24
89 velocity nodes,33 pressure nodes,24 elements.

 

 

 

 

 

 

 

\\
Y\\

. MESH 96
321 velocuty nodes.113 pressure nodes,96 elements.

 

 

 

 

MESH I48
489 velocity nodes,171 pressure nodes,148 elements

Mesh configurations

Figure 3

1000

1000

6C)

MESH 24 MESH 96
23.30 23.35
84.26 83.77
150.92 141.59
-—— 238.64

Table 1: Variation of u with mesh size
3max

MESH 24 MESH 96
22.45 22.44
64.97 63.84

103.45 99.54
--—- 158.07

Table 2: Variation of u3(0) with mesh size

MESH 148

 

83.72
141.34

237.48

MESH 148

63.81
99.32

157.08

61

It is clear from tables 1 and 2 that the results obtained by MESH
148 and MESH 96 are in excellent agreement for D S 2000. thus
indicating that either mesh is adequate for this range of the Dean
number.

In tables 3 and 4. we compared the same quantities (um. u3(0))
obtained from the finest mesh (148) with the results of previous
investigators. The tables indicate that our results are perfectly
consistent with those of previous studies for D S 2000. The agreement
is less than perfect for D > 2000 however. For this reason we restrict
the results in section 6.3 to the range D S 2000. and restrict the
values of R and e accordingly. thus for e = .1 we wake the restriction
R S 1118.03. and for e = .5. R S 500.

We believe that if MESH 148 is refined all around the boundary
(thus increasing the total number of velocity nodes to about 700

points). the flow in the entire laminar region can be accurately

compu ted .

605.72

1000

605.72
1000

2000

Present Study

83.72
96.53
141.34
237.48
317.72

354.72

Table 3:

Present Study

63.81
72.10
99.32
157.08
206.57

230.31

Table 4: Comparison of u3(0)

Collins & Dennis

(1975)
83.69
96.53

141.30
236.50

351.40

Comparison of 11:3

Collins 8; Dennis

(1975)

63. 70

99.00

154.70

224 . 70

Dennis
(1980)

83.67

 

141 . 10

 

314.80

Dennis

(1980)

63.78

99.16

203. 50

63

6.3 Results and discussion

The numerical results are summarized in Figures 4 - 11. Figures
4 - 7 show the contour lines of the axial velocity (u3) and the stream
function (W) of the secondary flow at R = 250 and R = 500 in the range
1 1

The stream function 4: is defined by

6x

(6.3) 9‘”- = - JJ'uz.
1

gig; = «5.1.
Clearly figures 4 - 7 show that the axial flow in diminished and the
secondary currents become stronger as the curvature increaces.

In order to study the effect of curvature on the flow when the
Dean number is fixed. we computed the solutions for D = 1000. D = 2000
in the range 0 S e S g; and the results are shown in Figures 8 - 11.
Observe that in Figures 8(a) - 11(a) (for e = O). the results are in
excellent qualitative agreement with those of Cbllins and Dennis (1975)
and Dennis (1980).

Notice however that in this case the secondary flow is diffused as
the curvature increases (see figures 9.11). The reason for this is
that for a fixed value of D. the Reynolds number R decreases as 5
increases (Recall that D = 4R5). We thus believe that from the
physical point of view. the pair (R.e) is a better set of parameters to
study the motion than (D.e). because when D is fixed. one cannot study
the effect of curvature independently of the effect of the pressure
gradient on the flow.

Table 5 shows the values of the flow ratio Q/Q0 at R = 250.
R = 500 as 5 increased from l—-

10
straight pipe at the same Reynolds number.

1
to 5. where Q0 is the flow rate in the

64

In table 6. we give the values of the flow ratio Q/QD for
D = 1000. D = 2000 and the same range of values of e. where now QD is
the flow rate in the limit of zero curvature at the same Dean number.

Tables 5 and 6 show that when R is fixed. the flow ratio Q/QO
decreases rapidly with the curvature ratio. and that the variation of
Q/QD with e is minimal. Thus for a fixed Dean number. curvature has
virtually no effect on the flow rate for practical values of 8. One
must keep in mind however that for a fixed Dean number D. the curvature
ratio 6 cannot be increased without decreasing the pressure gradient.
and these two competing factors seem to stablize the resistence of the
tube (which should increase with increasing values of R and e). and
this is what we believe causes the flow ratio Q/QD to be almost

GODS tant .

0:250:30

013.29%:

Table 5:

Table 6:

65

R = 250
0.748179
0.685238
0.645136
0.615091

0.591105

D = 1000
0.978992
0.959563
0.942704
0.928089

0.915433

The flow ratio Q/Q0

The flow ratio Q/QD

R = 500
0.648668
0.585044
0.545185
0.516832

0.494718

D = 2000
0.970511
0.950805
0.933770
0.919073

0.906423

66

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

68

 

 

.2

(b) R=500. Epsilonul

.1

(o) R-SOO. Epsilon=

 

1‘ 1m

 

(d) R2500. Epsilon=.4

 

 

 

 

 

5

(e) R=500. Epsilon=.

locity 0t R=500

0!

of oxi

Contour lines

Figure 6

69

 

 

 

(b) R3500. Epsilon-.2

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(o) 0=1000. Epsilon-.3 (R-322.75)

7O

 

 

(a) 0-1000. Epsilon-
(c) 031000. Epsilon-.2 (82395.28)

 

 

 

 

  

 

 

(f) D=1000. Epsilon-.5 (R=250)

-.4 (R-279.51)

(e) D-lOOO. Epsilon

 

 

D=1000

sty 0t

Figure 8

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

I / 1)), // L
l ‘l l

- I .
2.... g

 

 

 

 

 

 

 

 

film i

(e=)0 2000 Epsil =-.4(R- =5

‘ i llll “Hill/7?.

 

59. 02)

i I velocity at D=2000

Figure 10

8
l

I” i )0 2000. :psi L0 a "mm 0= 29100 59” i0! f =1 8:11:03 if

[/i ////\\\\ ///// W.

*3. i, f l l “ lllllllr 20.213.45.55) 111
(c) 0=2000. Eps.lon=.2 (R= —.790 56) (d) 0- 20 .sip

//./._// (T /l\\\\ ///// ( //\\\\

 

.7 Jill

)D= 2000. EpSl siol =.5 (R= =500)

73

 

 

 

 

 

 

(f) 0'2000. E sin

.02)

 

 

 

 

 

 

(o) 0-2000. Epsilon

 

OWUWER'7: A HJUUNNHHONEIXMTHXI

7.1 Introduction and formulation

The primry aim of this chapter is to give another example to
demonstrate the fact that our solution procedure applies to duct flows
in cross sections other than the circle.

Here we present a perturbation solution for the flow in a toroidal
pipe of elliptic cross section. Results and comparisons with numerical
solutions are given in the next section.

Consider the flow in a toroidal pipe whose cross section is an
ellipse with axes 2a. 2b where one of the axes makes an angle a with
the normal I of the center line of the tube. The geometry is shown in

Figure 12.

 

.— o -0-

 

 

 

d—r z .3,
Geometry and coordinate system for elliptic cross section

Figure 12
74

75

As was done in chapter 1. it can be shown that the dimensionless

momentum equations are

2 2
(7.1) (ugx-+v%)u+%—=-ng+Au-:—E(Bg;-6%)u-gE(Bu-5v)
65w2 6 562

 

(7.2) (ug-x-+v%y-)v- =‘%+AU':‘(B%‘537)V+§'(BU‘5V)

W:

(7.3) (ug;+ ‘85”? - 3651“» - 6v) =
and the continuity equation

(7.4) g—jufi) + gy—(w/E) = 0

Where u.v are the dimensionless velocity components in the directions

X1. X2 and w is the dimensionless axial velocity. In equations (7.1) -

(7.4). x = Xl/b. y = X2/b are the normalized independent variables. B =

-b3 a?

-——-is
4p02 6X3

the Reynolds number. where gig-is the pressure gradient. The boundary

cos a. 6 = sin a. e = b/E. G = (1- efix + 66y)2 and R =

 

conditions are

where F is the boundary of the elliptic region c x2 + y2 S 1. and c =
b2/a2 is the aspect ratio.
Observe that is we set a = b. a = 0 (thus c = 1. B = 1. 5 = 0).

equations (7.1) - (7.4) reduce to equations (1.5) - (1.8).

76

In addition to the assumptions of Chapter 1. we assume that e. < <

l and R = 0(1). Define the stream function 4' by

H

(7.5) u =

$123

Qty:
63"

an"
a]

substituting (7.5) into (7.3) yields the equation

2
(7.6) Aw=ffigffi+fzm%-5gy)w+§_'-%-§_!(pgy—+ag—;)w

6 , 6w 6w
whereifig—=§L€x--%37

Eliminating p from (7.1) and (7.2) by cross differentiation. and

keeping terms only up to C(52) we obtain

225 6

a 'A 9— %)w - E (a a - 537nm

+6 +6

2 1
(7.7) A .p = - p
Jc'a(y.x) ( 6’

Our discussion for obtaining a perturbation solution of (7.6) and
(7.7) will be brief since our formulation and solution procedure is
very similar to that of Srivastava (1980).

Theaxial velocity w and the stream function 4: are perturbed about

6:0
2
(7.8) w=wo+ewl+ew2+...
2
(7.9) w=ewl+e¢2+...

Substituting (7.8), (7.9) into (7.6), (7.7) we obtain the following obtain the

77

of linear equations upon comparing the coefficients of different powers

of e
(7.10) AVG = - 4R
2 a a 2
(7.11) A 4.1 = (5517+ 5 5;) wo
awn-'0) a a
(7.12) Awl = W + (55x- - 65)w0 - 4R(Bx - 6y)
2 a”'1"“"1) a a

(7.13) A (.2 = 777.3?)— + 2mg],- + ag—x-Mwowl) - 2(pa—x- - 5 %)(Aw1)

awl. 0) mm)

(7.14) Aw2'= (Bx - 5?) m)— * "T‘a y.x)
6N .w )

a a a a 2 '
we: ' bay—"1" 'omay‘ * 55;:Nr 4W" ‘ 53') + wo-
Equations (7.10) - (7.14) are solved subject to the boundary conditions

J = 0 on the boundary

cx +y =1

The exact solutions of equations (7.10) - (7.14) can be computed in
turn since in fact all the functions wJ. 403 are polynomials in x.y.

The solution for (7.10) can be easily seen to be

78

2

(7.15) wo = KR(cx2 + y - l). K

1+c

Now (7.15) and (7.11) give the equation

Azwl = 4 cK2R26x(cx2 + y2 - 1) + 4K2R2By(cx2 + y2 - 1)

whose solution is

(7.16) 11 = R26P1(x.y.c) + 22522(x.y.c)

where for example

2 2 2 2 2
P1 = X(p + p x + p y )(cx + y - 1)
oo 20 02

and the coefficients p . p . p are obtained by substituting (7.16)
00 20 02

in (7.11) and comparing the coefficients. This leads to a 3X3 matrix

equation for the coefficients p . p , p .
00 2O 02

Solutions of equations (7.12) - (7.14) can be obtained similarly.
Many more details can be found in Srivastava (1980).

It must be observed that. although the procedure for solving
(7.10) - (7.14) is very simple in principle, the process of comparing
coefficients and inverting matrices to obtain the coefficients of the
(polynomial) solutions becomes impossible to do by hand as one proceeds
to compute wl. $2. w2. A rather lengthy computer program was written

to carry out the calculation.

79

7.2 Results and discussion

'e begin by giving the results for an elliptic region where c = %.

thus the cross section is

2
E

2
4+y $1

Figures 13 - 15 show the first order stream lines 471/ R2 for
w w
a .. 0. a - '4" a = 2 respectively.

It must be mentioned here that the results of Srivastava (1980)
can all be easily reproduced if we choose a = 0 in our formulation.

All the comparisons given below are for the flow rate Q = [w
(where the integral is taken over the elliptic region). From (7.8) we

have
Q = ]w + efw + (:sz
0 1 2'

It is easy to see that w1 is an odd degree polynomial and therefore

does not contribute to the flow rate. thus

Q=Qo+52fw.

where 00 = fwo is the flow rate in the straight pipe.

It is easy to show that

w2 = RD1(x.y) + R3E1(x.y) + R5F1(x.y)

where D1. E1. F1 are polynomials and R is the Reynolds numbers.

Now the flow ratio Q/Q0 is given by

80

(7.17) = 1 + e2(D + R213 + R417)

Q__
QO

where D = ga-IDI. E = gs-[El. F = [F1 are constants that depend on

o°|"’

B. 6 and c.

It is worth mentioning here that the results given by Srivastava
(1980) for the flow ratio coefficients D.E.F are completely erroneous.
In table 7. we quoted Srivastava (1980) for the coefficients D.E.F for
different values of c. and in table 8 we give the correct values of the

same coefficients.

O

.°.°.°.°
8

1.0201

0

9.000
‘1
(II

1.25

D
-1.7617875
-1.6960690
-0.0052929

0.0201415
0.0208331
0.0213993

0.0223002

Table 7: Cbefficients in Equation (7.17)

81

E
-0.0037665
-0.0019303
-0.0010727
-0.0006594
-0.0006334
-0.0006085

-0.0003923

(Srivastava(1980))

D E
0.419048 -0.003770
0.133333 -0.001935
0.053724 -0.001077
0.022655 —0.000663
0.020833 -0.000637
0.019102 -0.000612
0.004678 —0.000395

Table 8: Cbefficients in Equation (7.17)

.(Present Study)

5665566

4.5155552.

.00001032
.00000212
.00000083
.00000038
.00000036
.00000033

.00000017

.00000386

.00000188

.00000039

.00000037

.00000035

.00000017

82

Finally we employed the method of chapter 5 to generate numerical

solutions of the problem with c = % %

and the flow ratio Q/Q0 was computed for each R in the range

a: forR=6.R=20.R=50.
-1— g e g .45.
10

The flow ratio 0/00 was then calculated for the same values of R
and 5 using the asymptotic formula (7.17) and the two sets of results
are sunmarized in Figures 16 - 18. As is obvious from Figure 16. the
results are in excellent agreement at R = 6. but the asymptotic
solution deviates considerably from the numerical solution as R
increases. (Recall that the perturbation solution is valid under the
assumptions that 6. << 1. R = 0(1))

Finally observe that at R = 6. the flow rate actually increases

with increasing values of e. In fact. for the given data (c = l

71"“:

‘11-), equation (7.17) takes the form

Q/Q0 = 1 + e2(l.809925 - 0.033573R2

- 0.000055R4)

and the expression in brackets is positive for R < m = 7.058.
The numerically computed values of Q/Qo reflect the same phenomenon
(see Figure 16).

It must be mentioned that the same phenomenon (that Q/Q0 increases
with e at low Reynolds numbers) was reported by Larrain and Bonilla
(1970) for toroidal pipes of circular cross section. and was later
confirmed by Wang (1981) who discovered that the same phenomenon occurs
for helically coiled tubes at small values of e and low Reynolds

numbers .

><

 

 

 

 

 

 

 

 

\Wy. .
/fln why

mnnmma stmm mom Q n 00 Awmnncnomnwo: mopcnwoav

 

 

 

 

 

 

 

mwmcnm Hm

84

 

 

 

 

mnnmma venom mom 9 u bmo Avmnncnomnwo: mochwoav

wwmcnm we

85

 

 

 

 

Stream lines for a = 900 (Perturbation solution)

Figure 15

86

 

 

 

 

 

0.00" q d «I “I q
owe 98 93 can owe 9w... ewe own. 98 96 98

Dunno... om. ocgmfiww oz 310$ 3).—.5

wwocnm Hm

meta"

gnu" 2:31.31 place:

manages 8.5.3

87

0.001
9.3.L
u.v?
0.3L
0.3..

   

gnu" 2.5.013. 3.5.2.

mar—o"

 

Pun

 

an... 0m. OCZMZMWW 02 3.02 3).—.5

wwosnm Hg

98 93 9.... 9mm 9% 9m. ammo 9mm 9»: 9mm 98

12333.0: sofas:

88

 

I,
ll
’
I,
I"
-"'l

«0:0

90 4 . 1 .
98 9mm 93 9:. 9m. 9m 9m. 9%.. 98 98 mg

 

 

find 0.... ocmfi>fiuw oz 3.02 «$.20

mwocnm Hm

gnu” Zea-zoo. 3.5—o:

u. 1.3.33.0: 8.5.3

Appendix A

The covariant velocity components Vi

 

 

 

 

 

jk
v1 =62‘Ul
.11 a(x1)2
vl =52”1
.22 a(x2)2
Vl = - K J §!1-- xZUI
.33 5x1
azu2
v211 = a(x1)2
v2 .2322.
22 a(x2)2
auz
V733="‘/j§'
.3 - 62 [fi]-2x__s_[_vi]
11 ' 1 2 1
a(X) Jj Jj ax if
v3 . .2 [.03
22 a(x2)2 J7
a U3
V?” = - K315; [F]

Appendix B

In this appendix we outline the method we used to

solve the matrix equation

 

 

 

 

r 1 T a r a 1
(3.23) A 0 (D) 0 U1 = r61
1 2 “r
0 A (D ) 0 U2 02
I)1 132 o o P o
.0 o o A. 3’3. .63.

 

 

4

The. reader is reminded that the matrices A and
A1 are symmetric, positive definite, and banded, thus
they were computed, stored compactly and then factored in
place using LINPACK routine DPBFA which carries out the
Choleski decomposition of a symmetric positve definite
banded matrix.

Observe that (3.23) is a block diagonal system.thus
US can be computed separately by

(3.1) u =A’1c

where [lea is found by using the LINPACK routine DPBSL
which carries out the back substitution step since the
Choleski decomposition of A is now available.

The discrete velocity vectors U1 , 02. and the
discrete pressure P were computed as follows :

system (3.23) is equivalent to the equations

(B-2) AU + DTP = c.

(3.3) DU = o

90

where n = (111.112). A = [A o J. c = [cl]. u = [01].
1

0 A

observe that

(11.4) p .-.- (11A'111Tr1 D A"1 c.-
and that once P is known, U = U) can be computed
u2

directly from (3.2) since

(13.5) U A-1(G - DTP).

The matrix nm'Hf' was computed by applying the

back substitution routine as many times as fiT has
columns, and then premultiplying the result A-IDT by D.
Observe that computing DA'IDT requires many back subs-
titutions for A, and then a matrix multiplication, but
the effort taken to accopmlish this is well justified
since several thousand solutions of (3.23) are
required.

Since DAihfl‘ is symmetric and positive definite

(it is a full matrix however) , its Choleski decompos-
ition was found by applying the LINPACK routine DPPFA,
which factors lyfdfir in place, and stores it compactly

since it is symmetric.

92

Now each time we solve equation (3.23). we proceed as follows:

(1) Find 03 by (3.1)

(11) Find P by (11.4)

(111) Find u = [111] by (3.5)

. U2

 

 

LIST 01" REFERENCES

LIST OF REFERENCES

[l] Bercovier,M. and Pironneau,0. (1979)

Error estimates for finite element method solution of the
Stokes problem in primitive variables. Numer. Math.
33 : 211-224.

[2] Berger,S.A. , Talbot,L. and Yao,L.S. (1983) Flow in curved
pipes. Ann. Rev. fluid Mech. 15 : 461-512.

[3] Carey,G.F. and Oden,J.T. (1986) Finite elements. The Texas
finite element series. Vols. 4,6. Prentice-Hall.

[4] Ciarlet,P. (1978) The finite element method for elliptic
equations. Studies in mathematics and its applications.
Vol.4. North-Holland.

[5] Collins,W.M. and Dennis,S.C.R. (1975) The steady motion of
a viscous fluid in a curved tube. Q.J. Mech. Appl. Math.
28 : 133-156.

[6] Collins, W.M. and Dennis, S.C.R. (1976) Steady flow in a
curved tube of triangular cross section. Proc. R. Soc.
London. Ser.A. 352 : 189-211.

[7] Cuming,H.G. (1952) The secondary flow in curved pipes.
Aeronaut. Res.Counc. Rep. Mem. No. 2880.

[8] Dennis,S.C.R. (1980) Calculation of the steady flow
through a curved tube using a new finite-difference method.
J. Fluid Mech. 99 :449-467.

[9] Dennis,S.C.R. and Ng,M. (1982) Dual solutions for steady
laminar flow through a curved tube. Q.J. Mech. Appl. Math.
35 : 305-324.

93

 

 

 

94

[10] Girault,V. and Raviart,P.A. (1979) Finite element

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[18]

approximation of the Navier-Stokes equations. Springer,
Lecture notes in mathematics No. 749.

Glowinski,R. (1984) Numerical methods for nonlinear
variational problems. Springer series in computational
physics.

Glowinski,R. and Pironneau O. (1979) On a mixed finite
element approximation of the Stokes problem.

(I) Convergence of the approximate solutions.

Numer. Math. 33 : 397-424.

Greenspan,A.D. (1973) Secondary flow in a curved tube.
J. Fluid Mech. 57 : 167-176.

Hood,P. and Taylor,C. (1973) A numerical solution of the
Navier-Stokes equations using the finite element
technique. Computers and fluids. 1 : 73-100.

Joseph,B. , Smith,E.P. and Adler,R.J. (1975) Numerical
treatment of laminar flow in helically coiled tubes of
square cross-sections. Part 1. Stationary helically
coiled tubes. AIChE J. 21 : 965-974.

Larrain,J. and Bonilla,C.F. (1970) Theoretical analysis
of pressure drop in the laminar flow of fluid in a coiled
pipe. Trans. Soc. Rheol.l4 : 135-147.

Le Tallec,P. (1980) A mixed finite element approximation

of the Navier-Stokes equations. Numer. Math. 35 : 381-404.

McConalogue,D.J. and Srivastava, R.S. (1968) Motion of
fluid in a curved tube. Proc. R. Soc. London Ser.A.

307 : 37-53.

 

 

[19]

[20]

[21]

[22]

[23]

95

Srivastava,R.S. (1980) On the motion of a fluid in a
curved pipe of elliptical cross-section.

Z. Angew. Math. Phys. 31 : 297-303.

Taylor,C. and Hughes,T.G. (1981) Finite element
programming of the Navier-Stokes equations. Pineridge
Press.

Temam,R. (1984) Navier-Stokes equations : Theory and
numerical analysis. Studies in mathematics and its
applications. Vol.2. North-Holland.

Truesdell,L.C. and Adler,R.J. (1970) Numerical treatment
of fully developed laminar flow in helically coiled tubes.
AIChE J. 16 : 1010-1015.

Wang,C.-Y. (1981) On the low-Reynolds—number flow in a

helical pipe. J. Fluid Mech. 108 : 185-194.