
Citation 
 Permanent Link:
 http://digital.auraria.edu/AA00003811/00001
Material Information
 Title:
 Approximating the incompressible Navier Stokes equations using a two level finite element method
 Creator:
 Nesliturk, Ali Ihsan
 Place of Publication:
 Denver, Colo.
 Publisher:
 University of Colorado Denver
 Publication Date:
 1999
 Language:
 English
 Physical Description:
 xii, 134 leaves : illustrations ; 28 cm
Thesis/Dissertation Information
 Degree:
 Doctorate ( Doctor of Philosophy)
 Degree Grantor:
 University of Colorado Denver
 Degree Divisions:
 Department of Mathematical and Statistical Sciences, CU Denver
 Degree Disciplines:
 Applied Mathematics
 Committee Chair:
 Franca, Leopoldo P.
 Committee Members:
 Mandel, Jan
Knyazev, Andrew Russell, Thomas F. Trapp, John
Subjects
 Subjects / Keywords:
 NavierStokes equations ( lcsh )
Finite element method ( lcsh ) Finite element method ( fast ) NavierStokes equations ( fast )
 Genre:
 bibliography ( marcgt )
theses ( marcgt ) nonfiction ( marcgt )
Notes
 Bibliography:
 Includes bibliographical references (leaves 128134).
 General Note:
 Department of Mathematical and Statistical Sciences
 Statement of Responsibility:
 by Ali Ihsan Nesliturk.
Record Information
 Source Institution:
 University of Colorado Denver
 Holding Location:
 Auraria Library
 Rights Management:
 All applicable rights reserved by the source institution and holding location.
 Resource Identifier:
 44076366 ( OCLC )
ocm44076366
 Classification:
 LD1190.L622 1999d .N47 ( lcc )

Full Text 
APPROXIMATING THE INCOMPRESSIBLE
NAVIER STOKES EQUATIONS USING A TWO
LEVEL FINITE ELEMENT METHOD
by
Ali Ihsan Nesliturk
B.S., Middle East Technical University, 1992
M.S., University of Washington, 1996
A thesis submitted to the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
Applied Mathematics
1999
This thesis for the Doctor of Philosophy
degree by
Ali Ihsan Nesliturk
has been approved
by
Andrew Knyazev
ii/i? M
Date
Nesliturk, Ali Ihsan (Ph.D., Applied Mathematics)
Approximating the incompressible Navier Stokes equations using a two level
finite element method
Thesis directed by Professor Leopoldo P.Franca
ABSTRACT
We consider the Galerkin finite element method for the advection
diffusion equation and the incompressible NavierStokes equations in two di
mensions, where the finitedimensional space(s) employed consist of piecewise
polynomials enriched with residualfree bubble (RFB) functions. The stability
features of the residualfree bubble functions for the NavierStokes equations
are analyzed in this thesis. We show that the enrichment of the velocity space
by bubble functions stabilizes the numerical method for any value of the viscos
ity parameter for triangular elements and for values of the viscosity parameter
in the vanishing limit case for quadrangular elements. To find the bubble part
of the solution, a twolevel finite element method (TLFEM) is described and
its application to the advectiondiffusion and the NavierStokes equation is
displayed. Numerical solutions employing the TLFEM are presented for three
benchmark problems. We compare the numerical solutions using the TLFEM
with the numerical solutions using a stabilized method.
m
This abstract accurately represents the content of the candidates thesis. I
recommend its publication.
Signed
Leopoldo P.Franca
IV
DEDICATION
To my family: Avse, Mustafa, Husameddin, Nuran, Halime.
ACKNOWLEDGMENTS
Foremost, I would like to thank my advisor, Leo Franca, for his en
couragement in preparation of this thesis, his support, guidance and time spent
in discussions throughout my studies at the University of Colorado at Denver.
I would also like to extend my thanks to Martin Stvnes for his useful ideas, and
to Andrew Knyazev, Jan Mandel, Thomas Russell, John Trapp for their efforts
as members of my dissertation committee, to Antonini Macedo for providing
several C++ routines, to Frederic Valentin for providing the visualization tool,
and to Saulo Oliveira for review of the dissertation.
CONTENTS
Figures ............................................................... ix
Tables................................................................ xii
Chapters
1. Introduction......................................................... 1
2. The Residualfree Bubble Method and
a Twolevel Finite Element Method............................... 10
2.1 ResidualFree Bubble Method....................................... 10
2.2 A Two Level Finite Element Method................................. 13
3. The Residualfree Bubble Method (RFB)
for AdvectionDiffusion Equations.............................. 18
3.1 The Problem Statement and RFBs.................................... 18
3.2 Stability Estimates............................................... 21
3.3 Error Estimates................................................... 28
3.4 A TwoLevel Finite Element Method................................. 31
3.5 Numerical Results ................................................ 34
4. The Stability Analysis of the RFB method for
the NavierStokes Equations ................................... 45
4.1 The RFB Method for the Incompressible
NavierStokes Equations....................................... 47
5. Computational Algorithm............................................. 72
5.1 Computational Algorithm on Global Mesh............................ 73
5.2 Computational Algorithm on Submesh................................ 83
vii
G. Numerical Results
89
6.1 Cavity flow...................................................... 90
6.2 Backward facing step flow........................................ 92
6.3 Flow past a cylinder............................................ 107
7. Conclusion........................................................ 123
Notation Index....................................................... 126
References........................................................... 129
viii
FIGURES
Figure
1.1 The Galerkin solution of the NavierStokes equation by equal
order linear interpolants ............................... 2
1.2 The polynomial bubbles change their height as Reynolds number
goes up but not its shape ............................... 9
1.3 The residualfree bubbles change their height and shape as Rey
nolds number goes up......................................... 9
2.1 A twolevel mesh................................................ 15
3.1 The solution of the problem (3.1) for various values of f....... 19
3.2 The definition of the line segment (x~,x+)...................... 23
3.3 The definition of ha............................................ 24
3.4 Statement of the rotating flow held problem..................... 35
3.5 Advection in a rotating how held................................ 36
3.6 Statement of the thermal boundary layer problem................. 39
3.7 Temperature cut at y = 0.25 39
3.8 Thermal boundary layer problem.................................. 40
3.9 Problem statement of an example with analytical solution. . 42
3.10 L2 and H\ errors for k, = 1: TLFEM ( ) and GLS ( ) .... 43
3.11 L2 and H\ errors for k = 0.01: TLFEM () and GLS ()... 44
4.1 Types of inflow boundary........................................ 56
IX
4.2 The point (x*,y*) is the rightmost point in A'.................. 56
4.3 x'j x is always nonnegative inside the triangle............... 60
4.4 The point (x,y) is upwind to (xj,y) in the streamline direction 61
6.1 The statement of the driven cavity flow problem ................ 93
6.2 The problem meshes tested: 400,1600 and 6400 elements .... 94
6.3 Pressure contours for Re = 400 ................................ 95
6.4 Pressure elevation for Re = 400 ............................... 96
6.5 The pressure along vaxis at x = 0.5 for Re = 400 97
6.6 Velocity held for Re = 400 98
6.7 Streamlines for Re = 400 99
6.8 Pressure contours for Re = 5200 100
6.9 Pressure elevation for Re = 5200 101
6.10 The pressure along yaxis at x = 0.5 for Re = 5200 ........... 102
6.11 Velocity held for Re = 5200 .................................. 103
6.12 Streamlines for Re = 5200 .................................... 104
6.13 Submeshes tested: Uniform vs. Non uniform .................... 105
6.14 The effect of submesh to the shape of bubbles................. 105
6.15 Pressure improvement by non uniform submesh................... 106
6.16 The statement of the backwardfacing step problem............. 108
6.17 Successively refined meshes for the step problem.............. 109
6.18 Backward step how: Pressure contours for Re = 150............. 110
6.19 Backward step how: Pressure elevation for Re = 150 Ill
6.20 Backward step how: Velocity held for Re = 150................. 112
6.21 Streamlines detailed in the circulation region................ 113
x
6.22 Backward step flow: Horizontal velocity for Re = 150............ 114
6.23 The statement of the flow problem past a cylinder............... 116
6.24 The problem meshes tested: 450,1800,7200 elements............... 117
6.25 The cylinder problem: Pressure contours for Re = 26............. 118
6.26 Pressure detail behind the cylinder for Re = 26................. 119
6.27 The cylinder problem: Velocity field for Re = 26................ 120
6.28 Streamlines detail in the circulation area for Re = 26.......... 121
6.29 The horizontal velocity for Re = 26............................. 122
xi
TABLES
Table
5.1 Entries of the contribution matrices
1. Introduction
The modeling of realworld flows like air flow, water flow, etc., to
predict physical quantities such as pressure, velocity, temperature, etc., is an
important subject in science and engineering. There are two main approaches
to deal with such problems:
Experimental approach
Computational approach
Although the experimental approach simulates the main characteristic features
of the flow, it is usually expensive, if not impossible to set up the experimen
t. Therefore the computational approach has been gaining more popularity,
especially due to the dramatic advances in speed of computation and storage.
It starts by defining the flow as a set of mathematical equations describing
the physics of the flow. Then numerical techniques are applied to obtain a
reasonable simulation of the physical phenomenon.
The NavierStokes equations model Newtonian fluids, such as air and
water. Their applications range from aerodynamics to nuclear reactors and
therefore they have received significant attention. However, their numerical
1
simulation is not a simple task and many numerical methods in the literature
have been applied to these equations with limited results. The finite element
methods are among the most successful and robust numerical methods in the
literature applied to flow problems [47, 39, 19] ( see also [55] for an introduction
to the FEM). Throughout this work, we present and compare a variant of the
mixed finite element method for the NavierStokes equations ( a comprehensive
treatment of the mixed finite element methods can be found in [8]).
Figure 1.1. The Galerkin solution of the NavierStokes equation by equal
order linear interpolants
Applications of the Galerkin finite element method to incompressible
NavierStokes equations were introduced in the early 1970s [50, 21]. It was soon
recognized that the use of equalorder interpolations for both velocity and pres
sure variables, which is the most desirable choice from the implementational
2
point of view, generates spurious pressure modes (see Fig 1.1). To design stable
numerical methods, there are two classes of discrete function spaces introduced
in the mixed finite element literature:
Circumventing BabuskaBrezzi condition [16, 22, 31, 41, 43, 46, 48, 7,
14, 42]
Satisfying BabuskaBrezzi condition [2, 5, 1]
The finite element methods circumventing the BabuskaBrezzi con
dition achieve stability by adding meshdependent perturbation terms to the
formulation. These terms enhance the coercivity of the formulation and en
ables the usage of the velocitypressure pairs known unstable because they fail
to satisfy the BabuskaBrezzi condition. However, the amount of perturbation
or the value of the stability constants is not known a priori and needs to be
tuned up by means of the error analysis and/or experiments [33].
The finite element function spaces satisfying the BabuskaBrezzi con
dition have been well studied and various authors suggested elements and pairs
of functions spaces that pass the BabuskaBrezzi condition. Let us consider
an example, the MINI element of Arnold et.al. [1]. The MINI element is
based on the choice of equalorder continuous piecewise linear interpolations
for velocity and pressure with the velocity variable enriched by bubble func
tions, polynomials vanishing on the element boundary. The bubble functions
3
can be eliminated at the element level because they are zero on each element
boundary. That elimination leaves behind a matrix formulation based on the
reduced space of equalorder linear interpolants. It has been pointed out that
stabilized PetrovGalerkin formulations of streamlinediffusion type can be de
rived for several problems by adding bubble functions to the velocity space of
piecewise linears with a stabilization parameter that depends on the shape of
the bubble functions [57, 9, 40, 3, 10, 27, 58, 61]. In diffusion dominant cases,
the MINI element formulation of the problem yields acceptable results. When
the inertia term is dominant, however, the method starts to produce degraded
results. Therefore the polynomial bubble functions are not sufficient for such
problems.
The problems caused by the advectiondominancv in the NavierStokes
model are the most difficult component of the whole problem. As it is pointed
out above, even the appropriate pair of approximation spaces may yield poor
results in such situations. It should be noted here that the difficulty is similar
to the advectiondiffusion equation when the Peclet number is high. In the last
several decades, a great amount of effort was devoted to the investigation of this
problem and various numerical methods developed to overcome this problem
can be found in the literature [62, 54], The most popular methods, in chrono
logical order, can be classified into five main categories: The classical artificial
4
diffusion (or viscosity) method [47], the streamline diffusion method [16], the
discontinuous Galerkin method [47], the methods using the Shiskintype, a pri
ori refined meshes [62] and the first order systems least squares method. The
classical artificial diffusion ( or viscosity) method modifies the diffusion term
in the problem to obtain nonoscillating solutions but it introduces a consider
able amount of extra diffusion in both streamline and crosswind directions and
produces less accurate solutions. It turns out that to reduce the oscillations
in the standard Galerkin method, it is appropriate to add a diffusion term
in the direction of streamlines only. Such method is usually called either the
streamline diffusion finite element method (SDFEM) or the streamline upwind
PetrovGalerkin (SUPG) method. This method introduces less crosswind dif
fusion than the classical artificial diffusion, however it may still be inaccurate.
Another way of dealing with the difficulty arising in the dominance of advection
term is the discontinuous Galerkin method. It has similar stability and con
vergence properties as the SUPG method. In practice, it performs somewhat
better than the SUPG. However, it is implementationally more complicated
and expensive as it requires to solve larger systems of linear equations. The
methods using the Shiskintype, a priori refined meshes usually works well if
the data is smooth [62]. The first order systems least squares (FOSLS) method
is an alternative to these problems, but it requires solving for more unknowns
5
than in the original equations [17, 18, 24],
In this work, to deal with the problems caused by the dominancy
of the inertial term and by the incompatibility of spaces due to the mixed
method nature of these equations, we consider the residualfree bubble method,
a more robust and accurate method, using a special type of bubble functions:
The residualfree bubble functions. The RFB functions minimize the residual
inside each element by satisfying the strong form of the problem locally. In
Brezzi et.al [6], it has been pointed out that adjusting the shape of the bubble
functions is equivalent to tuning up the stabilization parameter in the SUPG
like stabilized methods and they suggest a criteria for choosing the bubble
function that also works in advection dominated cases. It turns out that the
RFB functions also obeys the criteria given in Brezzi et.al. [6]: They are able
to adjust themselves to produce large gradients, which is crucial to gain the
stability according to [6], as well as their heights when the Reynolds number
is high (see Fig 1.21.3). The residualfree bubbles were first introduced by
Brezzi and Russo [15]. It has been proven in the limiting case ( i.e. the
coefficient of the highest order term in differential equation goes to zero) that
the addition of RFBs to the approximation space adds the stability to the
discrete problem using triangular and quadrangular elements of the advection
diffusion equation and to the discrete problem using triangular elements only
6
of the NavierStokes equations [15, 11, 64]. In the subsequent chapters, we
extend these results to quadrangular elements for the NavierStokes equations
and prove how the addition of the residualfree bubble functions to the space
of piecewise linears adds stability to the finite element discretization using
triangular elements, not only in the limiting case but in any critical case in
which the inertial term is dominant. We also note that the RFBs can be
employed in adaptive algorithms [65, 63], in obtaining other new stabilized
methods [35] and in deriving several numerical tricks [35]. See [26] on the
limitations of bubble functions. Before solving the matrix equation coming out
of the formulation, the residualfree bubbles should be computed inside each
element by a finite element method described in each element separately. We
do this by the twolevel finite element method. The twolevel finite element
method was first introduced by Franca and Macedo in [32] for the Helmholtz
equation and extended to the advectiondiffusion equation in Franca et.al. [34]
( it is also a part of this work and the details can be found in the chapter
devoted to the advectiondiffusion equation). To be able to apply the twolevel
method to NavierStokes equations, we introduce a vectorial version of the
method. Since the RFB computations are independent from one element to
another, these computations can be efficiently carried out in parallel machines.
Computational issues are not pursued in detail herein, instead we point out
i
the potential of the RFB method as it produces more accurate results than the
wellknown stabilized methods.
We organize the thesis as follows: In chapter 2, we give an overview
of the residualfree bubble functions (RFB) method together with a twolevel
finite element method (TLFEM) which we use to find the bubble part of the
solution in each element. In chapter 3, we review the RFB method applied to
the advectiondiffusion equation and present the stability results obtained for
high Peclet numbers. Numerical results obtained through the TLFEM is also
presented at the end of this chapter. In chapter 4, we investigate the stability
features of the RFB method for the NavierStokes equations for a wide range
of viscosity parameters. The computational algorithm we employ to solve the
RFB formulation of the NavierStokes equations will be described in chapter 5.
In this chapter, it is also notable that how the vectorial form of the TLFEM
can be incorporated into the algorithm to solve a set of nonlinear differential
equations. To illustrate the potential of the RFB method, we present some
numerical results in chapter 6.
8
Figure 1.2. The polynomial bubbles change their height as Reynolds number
goes up but not its shape
Figure 1.3. The residualfree bubbles change their height and shape as Rey
nolds number goes up
9
2. The Residualfree Bubble Method and
a Twolevel Finite Element Method
In this section, we describe the residualfree bubble method and a
twolevel finite element technique. In the residualfree bubble method, we en
rich the finite element space by an appropriate set of functions, the socalled
residualfree bubble functions to obtain a more stable and accurate discretiza
tion. Residualfree bubble functions are approximated by a twolevel finite
element method, which we outline in the second part of this section.
2.1 ResidualFree Bubble Method
where ff is a domain in 1R2 with boundary dQ, L is a linear differential operator,
u is the unknown scalar or vector valued function and / is a given source
function. We assume that L is such that this problem is well posed.
The system (2.1) can be transformed into a variational form by multi
plying the residual by a set of test functions and integrating by parts. To specify
a Galerkin finite element method for (2.1) out of the variational problem, we
Consider the boundaryvalue problem
0 on di},
/ in
(2.1)
10
begin by partitioning il into elements I\ in a standard way (e.g., no overlap
ping, no vertex on the edge of a neighboring element, etc.), then we choose a
finitedimensional space \\ which is related to the choice of partition and which
satisfies V), C V', where \ is the space of functions in which we seek a solution
of the continuous variational problem. We define h = inax^{diam(/v)} as the
mesh diameter. Then the Galerkin method reads: Find Uh \ 'h such that
where a(, ) is a bilinear form of the variational problem of the system (2.1).
For Galerkin methods that use bubbles, each Vh Vh is the sum of
a standard piecewise (isoparametric) polynomial and a bubble function that
is selected later to be residualfree. For the polynomial component of V/,, we
are mainly interested in the space of piecewise linears, for simplicity. Thus we
write
where v\ is the polynomial component of v/, and Vb is its bubble component.
Similarly we write
a{uh,vh) = if,vh) VufcGVfc,
(2.2)
Vh = Ui + vb Vvh
(2.3)
Uh = U1 + ub
(2.4)
where we recall that Uh is the solution of (2.2).
11
We require the bubble functions to vanish on the boundary dK of
each element K. For the particular case of residualfree bubbles, we further
define the bubble component vb of each vh to satisfy the original differential
equation on the interior of each K, i.e,
Lvh L(vx + vb) = f in I<.
This simply implies that for each I\ we have
Lvb Lv\ + / in A', (2.5)
vb = 0 on dK.
For general bubble functions including the residualfree ones, the van
ishing of the bubble functions on each dK allows us to use the static conden
sation procedure: in (2.2) take vb = Vb on K and Vh = 0 elsewhere to obtain
a(ux + ub, vb)K = (/, vb)K (2.6)
where the subscript K indicates that integration is restricted to the element
K. Our choice above of residualfree bubbles ensures that equation (2.6) is
satisfied automatically. Thus, after solving (2.5) strongly (or approximately)
it is enough to test the discrete problem by piecewise linear functions and the
numerical method that we implement reads: Find Uh = u\ + ub such that
a(ux,vi) + a(ub,vi) = (f,vx) for all vx V'i, (2.7)
12
where \ \ is our space of piecewise polynomials and is given in terms of ux
and / by (2.5). Thus the enrichment of the discrete finite element space by
residualfree bubbles results in a modification of the variational formulation of
the standard piecewise polynomial Galerkin method by the addition of a(ub,
to the bilinear form. The stabilizing effect of the additional term a(ub,Vi) will
be extensively investigated in the following chapters.
2.2 A Two Level Finite Element Method
To find the residualfree bubble part of the solution, %, we need to
solve (2.5). This can be done, using the fact that the linear operator L is linear,
by decomposing (2.5) into its constituting components. Thus we solve for the
bubble shape functions ipi and ff, with i varying from one to the number of
element nodes (rien):
>< 1 II (2.8)
fi = 0 on dl\ . (2.9)
where the s are the local basis functions for Ui and
Lip/ = / m A' (2.10)
V! = 0 on dK . (2.11)
13
Thus, if
(2.12)
then
ub =
+ ff ,
(2.13)
and substituting (2.12) and (2.13) into (2.7) we get the matrix formulation
where the c* s are the finite element approximations to the solution at the
nodes. It turns out that the first step of this procedure may be as complex as
trying to solve the original PDE problem. Therefore, we consider the following
strategy [32], At the preprocessing stage of a finite element code, we consider
approximating the differential equations for the bubble shape functions and
ff, (equations (2.8)(2.11)) bv another finite element method. In other words,
for each element we consider a submesh (a mesh defined for each element,
see Fig 2.1) where we will solve the equations (2.8)(2.11) by a suitable finite
element method. This can be done using a nonstandard method: the Galerkin
least squares method (GLS), for example. We will need to repeat this procedure
for each element and for each basis function. Once these problems are solved
for the residual basis functions, then the Galerkin method as given by (2.14)
is used to produce the approximation in the original mesh.
14
Submesh
Figure 2.1. A twolevel mesh
Let us now concentrate in the submesh problems for each residual
free basis function. We construct a stable nonstandard approximation to the
PDE problems given in (2.8)(2.11) over an element K. Let us denote by I\*
(see Fig 2.1) an arbitrary element in the sub mesh with diameter h* and by ip*
the basis function for a piecewise polynomial interpolation in the sub mesh, so
that our unknown bubble basis functions
p>\ and ip^, respectively, i.e.,
Pi' = X>fV; <2'15)
l
15
(2.16)
Here l runs over all unknown interior nodes in the submesh, say through N*.
The index i refers to the specific bubble function that we are trying to compute.
(Recall that we need to solve for a number of bubbles equals to the number
of basis functions used to span the piecewise polynomial defined in the global
element I\.)
To be more precise, we present the GLS formulation of the method
explicitly. The matrix formulation of the problem for the GLS method reads:
For each i (from 1 to nen) find c\l\ l = 1, 2,.., N* such that
for m = 1,2,.., iV*. We remark that the value of the stabilization parameter t
depends on the PDE problem and its discretization. Note also that the driving
force on the right hand side is the piecewise polynomial basis function ipi for
the element K.
Y {{L^k, ip*m + tLiP*Jk.}
Similarly we account for the basis function associated with /, iff :
Find c\*\ l = 1,2,..,N* such that
16
for m = 1,2,N* with the same stability parameter.
Once the constants and c\^ are found, we substitute them in
(2.15) and (2.1C) to get the approximate residual basis functions and .
These approximations are then used in the global problem
Y Ci[a(^i, ipj) + a(ipi, 0j)] = (/, ipj) a((ff, ipj) (2.17)
i=l
where tfi are replaced by cpf* and ipf is replaced by >pff. The global problem
now can be solved. The bubble part of the solution may also be recovered by
ub = Ycwi + Vf,
i
with in lieu of and Y in lieu of
17
3. The Residualfree Bubble Method (RFB)
for AdvectionDiffusion Equations
3.1 The Problem Statement and RFBs
We now consider the convectiondiffusion equation. Convection
diffusion problems have many practical applications. They also arise from the
linearization of the NavierStokes equations, therefore it is especially important
to devise effective numerical methods for their solution.
We start with a onedimensional example to present the difficulties of
solving the advectiondiffusion equation. We denote the diffnsivity parameter
by e in ID. Consider the following boundary value problem:
eu" + u = 1 on (0,1) (3.1)
u(0) = 0
u{1) = 0
The problem (3.1) has the solution
lx _1
e e
u(x) = x; (3.2)
1 e~
If e, the coefficient multiplying the highest order (second order) derivative term
in the equation, is big enough, the solution will be smooth and then standard
18
numerical methods will yield good results. However, as c is approaching to
0, the effect of the second order derivative term is only visible at the very
downwind end of the domain resulting in a sharp boundary layer (see Fig 3.1).
If this is the case, the standard numerical methods are only able to produce
reasonable results on a very hue mesh which is too expensive, if not beyond
our computational power at all (think about a 3D counterpart).
Figure 3.1. The solution of the problem (3.1) for various values of f.
To overcome this difficulty, several nonstandard numerical methods
have been proposed [16, 15, 29, 43, 46]. In this chapter, we develop a finite
19
element method within the Galerkin framework to obtain a stable and accu
rate discretization of the advectiondiffusion equation for a wide range of the
diffusivity parameter, which is the source of the singularity. Numerical results
confirming the theoretical findings will be presented at the end of the chapter.
The advectiondiffusion operator, in higher dimensions, can be written
as
L = k A + aV, (3.3)
where k is the diffusivity parameter, assumed to be a small positive constant
and a 7^ 0 is the velocity field. We assume for simplicity that a is constant
on H.   stands for the Euclidean norm in 1R2. We obtain the discrete problem
modified by RFB functions by applying ideas discussed in the previous chapter:
Find Uh = u\ + ub such that
a{uh,vi) = {f,vi) Vtq 6 V'i , (3.4)
where
a(uh,v j) = a(uuvi) + a{ub,vi)
= k(Vui, Vui) + (aViti, r>i) + Ac(Vuh, Vni) + (aVufc, Vi).
(3.5)
Thus the last two terms in (3.5) are the contribution of the RFB method.
20
Now we establish stability estimates for the bilinear form in a variety
of cases relative to the elements we employ and to the values of the diffusivity
parameter k.
3.2 Stability Estimates
Let Uh E Vh be arbitrary. Then Uh = u\ + Uf,. We decompose
further as Ub = ub + ub, where u is defined on each K by
Lub
ul = 0
Lu\ in I\,
on dK,
(3.6)
and u[ is defined on each element K by
Lu( = f in K,
uJb = 0 on dK.
(3.7)
Now set Uh = Ui + ub. From (3.5), the numerical method can be written as
a(uh,ui) = a(ub, ii\) + {f,ui) for all vx \\.
(3.8)
Thus, as ub part of bubble function is only related to the consistency, we prove
that the bilinear form a(, ) is coercive over the reduced space V/* x V'i where
\'h = {7/i : Vh E Vh} under the assumption that is the space of piecewise
linears or bilinears. Once the coercivity inequality is obtained, the generalized
LaxMilgram theorem [56] will be the main tool to assure the existence and
uniqueness of the discrete problem (3.4).
21
Case I: k > 0 and Triangular Element
Consider the case that k > 0 and \ \ is the space of piecewise linears
on a triangular discretization. Let Vh G W he arbitrary, and consider a(v}uv\).
We write  0.a' for the L2(K) norm, and  \ for the standard Hl(K) semi
norm. From (3.5), we get
a{uh,ui) = ^ ((Vhft, Vui)K + {aVuh,ui)K)
where we integrated (aVui, i)a, a(Vu, Vi)a' and (aU\)k by parts
and used the facts that Vi is a conforming finite element, space, = 0 on each
dK, a is constant and Ku\ = 0 on each element I\ since U\ is linear.
For each triangle I\, let dK~ = {x G dK : a n(x) < 0} be its
inflow boundary and dK+ = {x G dK : a n(x) > 0} its outflow boundary,
where n is the outwardpointing unit normal to dK. Assume that a n(x) is
bounded away from zero; then, excluding vertices, for each element K we have
dK = dl\~ U dK+, that is, no edge parallel to the flow exists.
K
HuiIi,a' + (aVwi,i)/f + Vu^A' + (a'Vu, uy)K)
K
K
Let be the reduced solution of (3.6), i.e., for each K,
(3.10)
22
Figure 3.2. The definition of the line segment (x ,x+).
Let (x,x+) be a generic line segment that is parallel to the flow and lies in a
single I\ with x~ 6 0K~ and x+ G dK+ (see Fig 3.2). Clearly
.~.o
u
1 rx
(x) = Ui(x ) ni(x) = / aVuids for x G (x ,x+). (3.11)
hd Jx~
Observe, by [51] and integral version of the CauchvScliwartz inequality, that
lim
KO
/ Ub= ub
J K J K
(3.12)
In (3.9), we replace u by ub, due to (3.12) and use (3.11):
o(fi*,Ui) T (llll.K ~ (pVllJi)
(3.13)
E ( Klilb + (aV,;)V f ( f *)'* ) (314)
K J x~
23
K
53 ( K\U]\hK H
(aVmi)2a' K //(A')
53 (kIu'Ii.a + ^lla'Villo,A')
(3.15)
(3.16)
where s is the arc length parameter, ha is defined as the length of the longest
segment, parallel to the flow and contained in K (see Fig 3.3) and n(K) is the
area of the triangle K. This is the desired coercivitv inequality for the case
when k ^ 0 and Vj is the space of piecewise linears on triangles I\.
Thus the numerical method is stabilized by the RFB method like
the socalled StreamlineUpwind Petrov Galerkin (SUPG) method. However,
in the present case, the stability parameter explicitly comes out of the RFB
method being ha/(3 a). At this point, one may wish to compare the perfor
mance of the new stabilization parameter ha/(3 a) with the classical definition
of the stabilization parameter h^/(2 a) in stabilized discretizations like SUPG
(/*a is given as the longest side for triangles and as in Harari and Hughes [44]
for quadrangles being diagonal of a rectangle, for example). In [15], Brezzi
24
and Russo perform such test and report qualitatively similar results for both;
the scheme using the parameter coming from the bubble procedure presents
some oscillations, vet produces sharper layers. However, when the classical def
inition of the stabilization parameter is employed replacing Hk bv ha, Franca
and Valentin have obtained quite good results for advectiondiffusionreaction
equation in a number of situations for a wide range of parameters [30].
Case II: 0 < k
The coercivity result can be extended to the case where 0 < k 1
and V'i is the space of piecewise linears on triangles K. Rewriting (3.9) and
using the result obtained in Case 1, we get
a{uh,ui) =  ('Ufc,aYu1)A)
K
=  (*,aVwi)K + (ub ul,aVui)K)
K
= (KluilIk + + (6 aVuOtf^ (3.17)
To obtain coercivity, we need to prove that the third term in this sum is dom
inated by the second term. This can be achieved by the virtue of the following
lemma.
Lemma 3.1 Let K be a fixed triangle in our triangulation. Suppose that no
edge in the triangulation is aligned with the direction of the flow and that
25
ub ub is defined on K by
/
L{ul   <) = kA ub in A',
ul 4 = 0 on dK~, (3.18)
<  4 = U l{x~) Ui(x+) on dK+.
\
Then there exists a fixed positive constant C, which is independent of k, ha
and a, such that
lflggo,ic
M2hlJ2
a3/2
K
I a V w 111 o, a' 
Proof: See [34] or Lemma 4.1 for a proof in a more general case.
Now we can finally obtain the desired inequality (see [34]).
Theorem 3.2 Assume that no edge in the triangulation is aligned with the
direction of the flow and that
k < /iaa min < 1,
1
64C2
for all triangles I\ in the triangulation. Then
a(uh,u\) > MmiIIk + 7Tbjlla'VwilliU )
K \ M /
Proof: Lemma 3.1 shows that
\(u0bulaVUl)K\
Kll2h]J2 K , ,
+ ] NaVuxUi*
a3/2
(3.19)
(3.20)
26
By (3.17) and (3.20),
K \ ii /
where we have used (4.31) for the last inequality and considered the cases
Case III: k > 0 and Quadrangular Element
Now we further present a similar result extended to the case where
k 0 and V'i is the space of piecewise bilinear functions on quadrangles K. In
this case, the coercivity is obtained in a different norm which is not equivalent
to the SUPG norm. This result can be found in [11], however it is stated
below for the sake of completeness.
Theorem 3.3 Assume that k
tion. Then
1 < 64C2 and 1 > 64C2 separately.
K
27
Proof: See [11].
3.3 Error Estimates
In this section, we study the convergence features of the Residualfree
bubble method. In section 3.2, we have obtained the coercivitv inequality
under some assumptions on the orientation of the mesh and the relation be
tween k, a and ha. The error analysis presented below relies on these stability
estimates. Convergence rates found below are confirmed by the numerical re
sults presented in the next section (see subsection An example with analytical
solution^ and figures therein).
Theorem 3.4 Let the solution uh be the solution of the RFB method and
suppose u Hs({}) f]HQ(Q), 1 < s < 2. If
K
then we have the following convergence rates:
K v 11 / K
28
Proof: Let tq be the linear interpolant of u, e/, = u\ ui and ?/ = ui u.
The proof follows the lines of the error analysis of the SUPG method and it is
given as follows:
Â£ (lVe,.HS.K + 7jblla^HL
< D{ui uh,eh)
= B(u\ u, ch) + B{u uk, eh)
= (a Vij,eh) + K{Vr],Veh)
= (r/, a Veh) + (Vr;, Veh)
Ii7/l 1,^.' Ila ^7e/io,A' + Vr/o,A Ve/(o,A
K K
= ^{KlY112 IMo,a (ha'y)1/2 a Veh0lA
K
+ ^K,'2Vr,,K K1/2Vek,h
<Â£
K 7
U,l
2 + bll a m2
0,A' + 2 l'a ve
h No,K
5IVdI!?,k + ^Vc2k),
where we have used the wellknown inequality
, a2 b2
ab <1
 2 2
for any real number a and b.
Now choose 7 = Noting that 7 is bounded below away from
zero, i.e. 7 > Co >> 0 for some constant C0, we can subsume the 2nd and 4th
29
terms into the right hand side of the inequality above. Thus we obtain
K
24 la
(321)
The first term in (3.21) can be estimated by interpolation theory [4]. For the
second term in (3.21), we remark that the coercivity inequality holds under
the assumption that
Now the right hand side of (3.21) can be estimated and applying the triangle
A similar independent result can be found in Brezzi et.al. [12]. Their
proof is very different than the result presented herein. Their proof has the
advantage of allowing that edges of elements may be aligned with the flow, how
ever it has the disadvantage of assuming that the source function / is piecewise
constant throughout the domain. A more general a priori error analysis can
be found in Brezzi et.al [13].
and thus, letting (3 = min {l, g^_}, we get
(3.22)
inequality, and the theorem follows.
30
3.4 A TwoLevel Finite Element Method
We now describe the application of the twolevel finite element method
to the advectiondiffusion problem using the Galerkin Least Squares (GLS)
method to solve each bubble problem. We remark again that to be able to
solve the global problem, we need to know what the bubble basis functions
are. Thus corresponding to equations (2.8)(2.11) we need to solve in each
element I\ of the global mesh:
a Vifi kALpi
a V?pi + kAipi in K
Â¥>*
0 on dK
(3.23)
and
a Vtpf kAlpj / in I\
/ = 0 on dK
(3.24)
for i = 1, ..,nen. However, solving these equations may be as difficult as solv
ing the original PDE problem. Instead, we find the approximate solutions of
these problems using a stabilized method like the Galerkin Least Squares in a
submesh and use these solutions in place of the exact bubble basis functions
as the discretization of the global problem is still in the Galerkin framework
with piecewise linears.
To construct the GLS approximation of the differential problems
31
(3.23)(3.24), we denote by I\* an arbitrary element in the snbmesh with di
ameter h* and by tpf the basis function for a piecewise linear interpolation in
the submesh, so that our unknown bubble basis functions and p>f can be
approximated by
i (3.25)
T = ETW (3.26)
i
Here l runs over all unknown interior nodes in the submesh, say through N*.
The index i refers to the specific bubble function that we are trying to compute.
(Recall that we need to solve for a number of bubbles equals to the number
of basis functions used to span the piecewise polynomial defined in the global
element A'.)
We then formulate the GLS method in matrix formulation as follows:
For each i (from 1 to nen), find c[l\ l = 1,2,N* such that
Â£ [(a VC, VC) + (VC, VC) + (a Wf.ra VC,)]
l
= (a + kA^k, ii*m + ra Vi/d ) (3.27)
and find such that
Â£ cfKa VC, C) + (VC, VC) + (a VC, ra VC)]
l
= (/,C + ra'V^), (3.28)
32
for l,m = 1,2In our computations, we use the stability parameter r
given as in [29]:
T(x,PeK(x)) = Â£(PeK(x)),
PeK{x) =
a(x) h*
6 k(x)
Â£{PeK(x)) = <
Pex(x), if 0 < Pex(x) < 1,
1.0, if Pe/v(x) > 1,
1/2
ia(x)i = m*)!5
Jl
(3.29)
(3.30)
(3.31)
(3.32)
Once the constants and cj^s are found, we substitute them into
(3.25) and (3.26), respectively, to get the approximate residual basis function
and tphj These approximations are then assembled to the global problem
^ry[(a V^i, Vj) + (V^, V^) + (a V^,^)] =  (a Vtp/,^)
2=1
(3.33)
in place of the exact residualfree bubble basis functions
The global problem now can be solved. The bubble part of the solution may
be recovered by
ub = + (pf,
l
with in lieu of
(3.34)
33
3.5 Numerical Results
In this section we report three series of experiments for convection
diffusion problems with the twolevel finite element method introduced in the
previous section. We compare our results with those obtained by the GLS
method used in the entire domain.
Advection in a rotating flow field
This problem is defined on a unit square of coordinates 0 < x, y < 1,
where the flow velocity components are given by
a\ = (y 0.5), a2 = x 0.5
Along the external boundary u = 0 and along the internal boundary (OA),
u = (cos(47t?/ 7r) + 1), 0 < y < 0.5
(see Fig 3.4 for problem statement). The diffusivitv is k = 106. A uniform
mesh with 30x30 elements is employed with bilinear elements Q1. Twolevel
results are obtained by subdividing each element into 7x7 subelements: each
subelement is generated by a linear mapping from a uniform subdivision with
7x7 subelements in the reference square domain. In the following results, we
have not added the bubble part of the solution to the bilinear part in the
twolevel method.
34
Elevation plots are shown in Fig 3.5 for GLS and TLFEM methods.
In this problem, both methods perform well and the results are almost indis
tinguishable.
y
u=o
Figure 3.4. Statement of the rotating flow field problem.
Thermal Boundary Layer Problem
Let us consider a rectangular domain of sides 1.0 and 0.5, subjected
35
I
I
i
i
i
0 0
i
!
0 0
Figure 3.5. Advection in a rotating flow field.
GLS
TLFEM
36
to the following boundary conditions (see Fig 3.6 for problem statement):
u 1, when < x = 0, if 0 < y < 0.5,
1 y = o.5, if 0 < x < 1.0,
u 0, when y = 0, 0 < x < 1.0,
u = 2 y, when x = 1, 0 < y < 0.5.
The flow components are
ai = 2y, a2 = 0 in Q
and the diffusivity is k = 7xl04. This problem may be viewed as the sim
ulation of the development of a thermal boundary layer on a fully developed
flow between two parallel plates, where the top plate is moving with velocity
equal to one and the bottom plate is fixed. Taking the top plate velocity as the
characteristic flow velocity, we have a Peclet number Pe = 714 for this flow .
We employ a rectangular mesh consisting of 21 equally spaced nodes
in the xdirection, 6 nodes uniformly distributed in the interval 0 < y < 0.1
and 6 nodes equally spaced on 0.1 < y < 0.5. We also subdivide each element
into 20x20 subelements: each subelement is generated by a bilinear mapping
from a uniform subdivision with 20x20 subelements in the reference square
domain.
In this problem, both GLS and TLFEM perform similarly except in
37
the neighborhood of the outflow boundary layer (see Figures 3.73.8). TLFEM
is a little oscillatory, yet GLS is overdiffusive therein.
An example with analytical solution
In the convergence analysis, we consider a simple problem on a unit
square that can be solved analytically, so that we can compare the GLS and
TLFEM results with the exact solution. Consider the unit square subjected to
the following boundary conditions (see Fig 3.9 for problem statement) :
/
z = 0, if 0 < y < 1.0,
when y io, if0
(3.35)
x = 1, if 0 < y < 1.0,
\
u = sin(7ra:), y = 0, 0 < x < 1.0.
The flow components are
oq =0, 02 = 1.0 in Q .
Using separation of variables, the exact solution is given by :
1
emiy em2y
) sin(7rx)
ll =
where
1 \/l + 4k27T2
38
Figure 3.6. Statement of the thermal boundary layer problem.
1.08' 1.08' /I
0.91' 0.91"
0.74 0.74"
0.57' 1 0.57"
0.4 H 0.4 H
0 0.3 0.6 0.9 0 0.3 0.6 0.9
Figure 3.7. Temperature cut at y 0.25
39
GLS
tlfem
Figure
s.8. Thermal boundary layer problem.
40
and
1 + \/l + 4k27T2
m2 = V
Ik
We consider two different values for the diffnsivity, namely, k = 0.01
and k = 1.0. We compute the L2 and F/x errors using 25point Gaussian
quadrature on each global element except for the rows of elements between
0.9 < y < 1, where there is an outflow boundary layer for k = 0.01. The
TLFEM approximation is obtained by subdividing each global element into
5x5 subelements. The results are shown in Figures 3.103.11. We obtain second
order convergence in the L2 norm and first order in the Hi norm when k = 1,
but for k = 0.01, the convergence rate slightly decreases for both and the GLS
method is slightly more accurate than TLFEM.
41
y
1
u=o
u=o
u=o
Advection=(0,l)
U=sin(7ix)
x
Figure 3.9. Problem statement, of an example with analytical solution.
42
L2 errors
Hi errors
Figure 3.10. L2 and Hi errors for /c = 1: TLFEM () and GLS ( )
43
2
10
2
10
3
Li2 errors
Hi errors
Figure 3.11. L2 and Hi errors for k 0.01: TLFEM ( ) and GLS ( )
4. The Stability Analysis of the RFB method for
the NavierStokes Equations
The incompressible NavierStokes equations are among the most im
portant models in Computational Fluid Dynamics. There are many applica
tions of these equations in many scientific and industrial areas. The Navier
Stokes equations are mainly used to solve problems in heat transfer, aerody
namics of vehicles, aerodynamics inside motors and problems originating from
marine and meteorological sciences. Therefore these equations are of interest
to many researchers and a large number of work can be found in the literature
(see, for example, [37, 38, 62, 59, 64, 49, 20] and references therein).
On the other hand, the incompressible NavierStokes equations are
difficult to solve especially for problems of practical interest where the Reynold
s number is high. Its mathematical study is very hard and incomplete. The
uniqueness of the solution in 2D or 3D are assured under very restrictive condi
tions if not impossible at all. The additional difficulties arising in the analysis
of NavierStokes equations can be expressed in three items:
The presence of the nonlinear term
The incompressibility constraint
45
The presence of boundary layers at high Reynolds number
In the residualfree bubble analysis of the NavierStokes equations,
we consider the linearized version of the steadystate problem thus not dealing
with the nonlinearity in derivation of the coereivity inequalities. However, we
work with the full nonlinear version of the problem in numerical tests. There
the nonlinearity will be resolved by the Newtons method.
The second difficulty, the incompressibility constraint, is related with
the choice of approximation spaces for velocity and pressure variables. One
cannot choose any combination of finite element spaces, especially finite di
mensional spaces of equal order which make the implementation easy. The
difficulty can be solved by enriching the velocity space by an appropriate set of
functions, bubble functions, for example. As polynomial bubble functions (or
even any bubble function; a continuous function zero on the element boundaries
and nonzero in the interior of the element) fulfill the requirements imposed
by the incompressibility constraint, choosing residualfree bubbles will equally
well satisfy this constraint.
At high Reynolds number, the Galerkin formulation of piecewise poly
nomial spaces exhibit pressure oscillations. The mechanism is similar to the
advectiondiffusion equation when the Peclet number is high. The deficiency
in the Galerkin method will be avoided by using the methods developed for
46
the advectiondiffusion equation presented in the previous chapter.
Now we derive the stability estimates for the discrete problem for
the critical values of the Reynolds number where boundary layers are present.
Triangular and quadrangular meshes are treated separately.
4.1 The RFB Method for the Incompressible
NavierStokes Equations
Consider the linearized incompressible NavierStokes equations given
by:
(Vu)a i/Au + Vp = f in Q,
V u = 0 in Q,
u = 0 G o
where Q is a domain in 1R2 with boundary dQ, u is the unknown velocity field,
p is the unknown scalar pressure function and f is a given source function, u is
the viscosity parameter and a is the advection vector assumed to be constant
throughout the domain.
The continuous weak formulation of the problem can be stated as :
Find u e V = {H^(tt))nsd and peP = C(Q) D L2{Q) such that
B{u,p,\,q) = F{v,q) V{v,
where nsd is the space dimension of the problem and the bilinear and the linear
47
forms are given by
B(u,p\v,q) = ((Vu)a,v)+ 2i/(Vu, Vv) (V v,p) + (Vu,q),
F{v,q) = (f,v).
To introduce a finite element method, we begin by partitioning the do
main into triangles or quadrangles in a standard wav ( no overlapping, matching
vertices across inter element boundaries, etc..) and define the conforming finite
dimensional spaces that we wish to work with:
V, = {v Â£  v ^Â£ (Pi(K) B(I\))nsd, K Â£ T/J,
Ph = {p Â£ C(Q) D I p fte Pi (A'), K Â£ T/J,
where B(K) is the space of the residualfree bubble functions. For a given
K Â£ Th, any function in B(K) is zero outside I\ and satisfy the following
system:
{(V(ui + u6))a ;/A(ux + u6) + Vpi = f in K ,
(4.3)
ip, = 0 on dK.
For the problem (4.2), the standard Galerkin method enriched by
residualfree bubble functions now reads: Find u* = Ui + u;, Â£ \\ and px Â£ Px
such that
B(uh,pu\h,qx) = (f, v/j) Vjvfc,^} Â£ Vh X Pi , (4.4)
48
where
B{uh,P\\vh,qi) = ((Vu/l)a,vfc)+ 2i/(Vu/>, Vv,*) (V v^,/q)
+ (V Ufc,^).
The problem (4.4) can be decomposed as follows:
B{nh,pu vl50) = (f.Vi) V{vJ e Vi , (4.5)
B{u.h,pi,0,qi) = (f,0) VjgJePi, (4.6)
B(uh,pi; v6,0) = (f, v6) V{vj} G Vj, (4.7)
where the subscript K indicates that integration is restricted to the element K.
Note that the problem (4.7) is automatically satisfied by our choice of bubble
functions; it is zero on the element boundary and satisfies the momentum
equation in (4.3). This special choice of functions reduces the bilinear form
(4.4) to the following simpler form: Find u* = Ui + ip, Â£ V), and pi Â£ Pi such
that
B(uh,pi, Vi,qx) = (f, Vj) V{vj, qi} G V\ x Fj , (4.8)
where
B(uh,pi;vuqi) = ((Vu/l)a,v1) + ^(Vuh, Vvj) (V Vi,px) + (V uh,qi).
(4.9)
49
We now turn to the stability analysis and prove the coercivity inequal
ities for the bilinear form (4.9) in a variety of cases relative to the elements we
employ and to the critical values of the viscosity parameter u.
Let rq, Â£ Vh be arbitrary. Then u/( = iq + iq. We decompose iq,
further as iq, = u + u(, where u is defined on each K by
Au = Lui Vpi in A',
(4.10)
< = o
on OK,
where L is the vector form of advectiondiffusion equation, that is,
Au = (Vu)a t'Au .
From (4.8), the numerical method can be written as
B(uh,Pi; vi,9i) = B{ufb,pivi,qi) + (f, Vi) for all {v^rq} Â£ Vi x Pu
(4.11)
where we set iq = Ui + u.
Let the reduced space be \\ = {vq, : vft Â£ Vh, Vh = Vi + v, v[ =
0}. As the tijÂ£ part of the solution is related with the consistency of the
method, we only need to prove that the bilinear form B(.,.;.,.) is coercive
over Vh x P x Vq x Plt under the assumptions that Vi and Pi are the space
of piecewise linears on triangles I\, that the triangulation is regular, and this
regularity is uniform as h > 0. We shall also make assumptions later on the
50
orientation of the mesh and the relationship between v, a and ha, where ha is
defined as in the previous chapter and of the same order with the diameter of
the element I\, hx ( That is, there is a constant Cx such that hx Cha).
Case I: i/ > 0 and Triangular Element
Let u/i Vh be arbitrary, and consider D(uh,p\',ui,qi). We write
II o,d fr the (L2(D))nsd norm on any domain D, and  iiA for the standard
(Hl(K))nsd seminorm. From (4.9),
#(u/i,Pi;ui,pi) = (^(Vuft, Vu!)A^ + ((Vu/j)a, U!)A (V ui,p!)
K
+ (V ufc,pi))
= X! (,y(Vui, Vu^x + ((Vuj)a, ui)A + i/(Vu,Vui)k
K
+ ((Vu^)a, u^x + (V ug,pi))
= (HVaiH + ((Vu1)a,u1)A (u, (Vu^a + Ypi)/r) ,
K
(4.12)
where we integrated by parts and used the facts that u = 0 on each dK, a is
constant and Aui = 0 on each K since Ui is linear.
For each triangle K, let dK~ = {x G dK : a n(x) < 0} be its
inflow boundary and dK+ = {x dK : a n(x) > 0} its outflow boundary
as before, where n is the outwardpointing unit normal to dK. Assume that
a n(x) is bounded away from zero; then, excluding vertices, for each K we
51
have DI\ = OK U 0K+.
Let u be the reduced solution of (4.10), i.e., for each K,
(Vu)a = (Vui)aVpi in A',
on OK
(4.13)
u = 0
Let (x,x + ) be a generic line segment that lies in a single K in the streamline
direction with x_ E dK~ and x+ E dh+ (see Fig 3.2). Integrating along the
line segment (x~,x+), we obviously get
u(x) = j~j (^J ({^ni)a + '^Pi)ds^j for x e (x ,x+)
= n((Vu!)a + Vpi)a x x for x E (x_,x+). (4.14)
a
On the other hand, using that the volume of a pyramid is given by
Height x Area of base
\ olume =:
we observe that
dx
ha n(K)
3
where x E (x ,x+). Hence, recalling the relation (3.12) and that ((Vu1)a +
Vpi)a is constant in an element, the third term in (4.12) can be written, with
52
an explicit statement of stabilization parameter, as
EW'(Vu>) + vpOk
K
EW'(Vu>) + v',')
K
K
/\ x x , (Vu[)a + V/q
K
K
V ^((Vu^a + V/q)2A
L' a
x
X
((Vui)a + Vpi)2ft
^^rH(Vui)a +v^Ho.x
K 13131
(4.15)
We finally make nse of this result and note that a is constant throughout the
domain which implies that ((Vui)a, Ui) = 0 (as Ux = 0 on dQ), to acquire the
desired coercivity result from (4.12);
Case II: 0 < is
For a small positive value of is, the viscous term should be taken
into account in the analysis. This makes the analysis harder compared to the
previous case. Regarding the bilinear form in the present situation, we basically
need to bound the additional terms in the bilinear form appropriately to get
the coercivity inequality.
K
53
Thus, using the result (4.15) from the previous section, the bilinear
form in the present situation can be written in the following simplified form:
B(uh,pi;uupi) (4.17)
= E (HuiI?., + ^IKVuOa + VPi^ + (u u, (VUl)a + VPl) .
We need to show that the third term in (4.18) is dominated by the
second term. This can be accomplished by realizing the following lemma.
We begin by rewriting the differential equation that u satisfies on each
element K: By (4.10) and (4.14),observe that
ms ug) o o <2 ^1 1 II in A',
1 O C> ug)(x) = 0 on dK~
(ug <)(x+) = E((VUl)a +Vpi) 1/, ; x+ X  on dK+
For notational convenience we set z = u uÂ£.
Lemma 4.1 Let K be a fixed triangle in our triangulation. Suppose that no
edge in the triangulation is aligned with the streamline direction and that z is
defined on I\ by
Lz = I/Aug in A',
z = 0 on <9 A', (4.18)
z = 4>(x+) on <9A'+,
54
where
0(x)
1
((Vu,)a + Vpi)A x x
(4.19)
a
Then there exists a fixed positive constant C3, which is independent of v, ha
and a, such that
Proof:
Define the function T on K by
/
L<4> = 0 in K,
< 4> = 0 on dK~,
= 4> 011 dK+,
so that is a smooth boundarvlayerlike extension of
w = z . Then
If dK~ is a single side of K, then r'Au = 0 on I\, and then w = 0 all over
the domain I\ implying 2 = but in general this is not the case. Therefore
we assume that dK~ consists of two edges of K (see Figure 4.1).
Lw = J'Auj) in K,
w
0
a
a
Singlesided inflow boundary Twosided inflow boundary
Figure 4.1. Types of inflow boundary.
Let us denote the cartesian coordinates in two dimensions by (x, y) .
Since a is constant on K, without loss of generality we assume that a = oiei
for some a\ > 0, where ei is the unit vector in the direction of the positive x
axis.
We first bound w in terms of the crosswind derivative of u. Then de
velop an argument to bound the crosswind derivative of u in terms of stream
line derivative of this function.
Choose (x*,y*) dK+ such that x < x* for all (x,y) G I\ (see
y
a
x
Figure 4.2. The point (x*,y*) is the rightmost point in I\.
56
Fig 4.2). Multiply the ith component of Lw = by the ith component
of (x* x)w(x, y) and integrate over I\, for i G {1, 2}, to get
/ (x* x) (Lwi)(x, y) Wi(x, y) dx dy
J K
f ~ (l c)
= I u{x* x)\Vwl(x,y)\2 + ^{x* x)w2{x,y) dx dy
= f [v{x* ~ x)\Ywi(x,y)\2 + (al/2)w2{x,y)] dxdy, (4.20)
J K
where we integrated by parts and used wl = 0 on dK. Next, integrating the
lefthand side of (4.20) by parts and using the inequality
a b <
a2 + b:
we get
/V *) v) *(*, y) dxdy
7, J x)\Vwi(x,y)\2 + ^w2{x,y) dxdy
<
+
' K
is(x* x)\Vulb(x,y)\
2^
ax
I Vu6(ar,?/):
dx dy.(4.21)
Combining (4.20) and (4.21), we have
\ fJjwHx>y)dxdV
ll \u(x* ~ x)\Vwi(x,y)\
1 J K 1
If u(x* x)\VKb(x^y)\2 + ~\V{1lb{x,y)\2
2 Jk L Oi J
2Iw2(x,y) dxdy
dx dy.
57
Consequently, we infer that
w2(x, y) dx dy
<
<
v(x* 
uha +
x) +
2^\
i )
2 v2
a i )
Vu
o
i,b
\^ulb{x^y)\2 dx dy
(x,y)\2dxdy.
Now
\Wi
0 ,K
< 21 +1Ir
a
1/2
1
1
(4.22)
for each i G {1,2}. Thus
IMk* < 2v"2 j lVugo,A: (423)
We now bound Vuo,A' hi terms of its streamline derivative. Recall that
we assume that dK~ consists of two edges of K and the edges of each I\ are
bounded away from being parallel to a = aie^ This means that we can write
the equations of the two edges E\ and E2 in dl\~ as i' = rrijy + kj for each
j G {1,2}, where the rrij and kj are constants that satisfy \m,j\ < C\ for some
fixed constant C\.
If (x, y) lies in the part of K that is downwind to a point (x~, y) on
the side Ej, j G {1,2}, then, letting A, = ((Vui)a + Vpi),a, from (4.14)
cl
58
we have
ulb(x,y) = Aj x x
= A, (x x )
= Ai(x {TTlflJ + kj)),
since a = Gqei and a\ > 0. Then it is easy to see that
du0lb(x,y)
dy
\\imj\ < CiAi = Ci
9ui,b(x, y)
dx
(This inequality bounds the crosswind derivative of uj) in terms of the stream
line derivative of this function.) Hence, for i 6 {1,2},
Vu6o,A < (m6)io,A' + 11 I o,A'
< ((1 + Ci)/ai)aVu6o,A'
= ((1 + Ci)/ai)((Vui)a + Vpi)jo,A',
where we have used (4.13). Invoking this inequality in (4.23), we obtain
(v^hT
\ A2
Next we bound 11^o,at As an application of maximum principle
for the advectiondiffusion operator, this will be accomplished by designing a
suitable barrier function for T: A function such that if it satisfies <3> < 6 on
w lo, AT 5: 2vz (1 + Cl)
V
+ d
(Vvi)a + Vpi
0,A'
(4.24)
59
OK and A<1> < LG on K then  < 9 on all of A; see [60]. At the end, an
upper bound for the barrier function will also be the desired bound for .
We shall do this for the case where dK+ consists of two edges of K;
the case where 0I\ + is a single edge is similar but easier.
Figure 4.3. xj~ x is always nonnegative inside the triangle.
Let the two edges in 0K+ be E\ and E2, where the equation of Ej is
xL = m' y + kj for j = {1,2}, and the rn'] and k'} are constants. Extend each
edge Ej to form a complete line (which we still call Ej, see Fig 4.3). Given
{x,y) Â£ I\, define x+ = x~j(x, y) by the requirement that (x~j ,y) lie on Ej (see
Fig 4.4).
60
y
a
x
Figure 4.4. The point (x,y) is upwind to (x*,y) in the streamline direction
We set M and N such that
M, = max {Id), I}
dK+ 1 J
where 0 is given by (4.19) and
N = 1 + {m[)2 +
Define the function Oj(x,y) by
9j{x,y) =Me'(lrl)W > 0, for j e {1,2}.
Observe that
9
M on Ej,
< 0j < M on dK+ \ Ej,
(4.25)
6j < M in K,
61
For each j, a calculation shows that
a? ( l + (m')2\
= Jr1)
e, > o,
since N >1 + (m')2. Now set G(x,y) = G\(x,y) + G2{x,y). Then
LG = LGi + LG2 >0 = 1$
on K, and if (x,y) Â£ 9A'+, then
0(x+) = 0!(x+) + 02(x+) > M > $(x+) = 0(x+),
because, by construction, G3 = M on E3] also G > 0 = $ on dK~. That is, we
have shown that G is a barrier function for $ on K. The maximum principle
now implies that $ < G on K. We make use of the wellknown inequality
(a + b)2 < 2 (a2 + b2) for all a,b e 1R,
62
to estimate the norm of (I> bv the barrier function:
Iloilo,A' < ll^llo.All^+^llo^
<
2 ^ 1/2
Mieax{x+x)/(vN) + y^eai{xtx)l{uN)^ I
n (x"^" ,?/)
/ 2(Mfe2a^x)l^N)\ dxdT
<
'x+Â£<9A'+ J (x~ ,y)
[ [{ y) 2 (Mfe2ai^xm,/N)) dxdT
Jx.+GdI\+ J(x~ ,v) '
+
1/2
'x+6<9A' + J (x~ ,y)
C (x^ >2/)
/ 2 (A/?e2ai w+fcl*)/(^JV)) dxdr
J (x~ .V)
<
'x+5A'+ J (x~ ,y)
+
[ [( V) 2 (M?e2ai{m2y+k2~x)/{l/N)) dxdT}
Jx+edK+ J(x~,y) I
1/2
< < I A/i f(le2ai(mi"+fclI~)/(,/Af))) rfr
ix+edK+ a\ ' '
L
+ / AA
'x+Â£SA+ 1
g2al(m2!/+fc2x (IY
1/2
< .w
^dr}1'2
x+eaA"+ i J
(4.26)
63
where we have used the fact that for j 6 {1,2},
1 e2ai(x+x) < x for x e gK
and
e2a1(xtx+) < j for x+ e QK+
Since ((Vv^ai ^Pi)i\K is constant on each I\, it follows from the regularity
of the triangulation that
Mi = sup$ < ha ((Vu1)a + Vpi)i\K \/a\
dK+
< hK ((Vui)a + Vpi)i\K /ai.
(4.27)
(4.28)
Thus the bound obtained for above can be estimated further:

iIo,it <
2Mi'/Nv1^2hlJ2
1/2
<
2\TNvll2h%2 (VVl)a + Vpi)
i if
3/2
<
2v^CV/2/4/2((Vv1)a + VPl)
i 0,if
3/2
(4.29)
using the fact that
ha ((VVl)a + VPl)t\K  < C2 ((VVl)a + VPl)
i o,if
(4.30)
64
for some constant C2 depending on the minimum angle condition in the trian
gulation. To see the relation (4.30), observe that
ll((Vv1)a + VpiJiH* = ((VVl)a + VpO.U ( / dQ)1/2 = ((Vv1)a + C^lhK
Since hx = C ha for a non degenerate triangulation, by the triangle inequality,
from (4.23) and (4.29), where C3 = 2\/2 (1 + C\) + 4 y/N C2 C.
It is crucial that the bounding coefficient above is a multiple of the
viscosity constant. This fact actually enables 11s to embed the last term in
(4.18) into the third one in the same equation which is known to be coercive.
The next theorem uses this result to establish the coercivity inequality.
Theorem 4.2 Assume that no edge in the triangulation is aligned with the
streamline direction and that
v < ha\a\ min jl, ^ (431)
for all triangles I\ in the triangulation. Then
BiUh^puUuPi) > ( Hul 1,K +
K
12 la
(Vu1)a + Vp13 A
65
Proof: Lemma 4.1 shows that
1/2 7 1/2
(u u, (Vui)a + Vpi)A < C3 ( r 3y2f jn I 11(^L11)a + ^iPillo.K
By (4.18) and (4.32),
a3/2
(4.32)
h
B{uh,pi\uupi) = y'(^ui?^ + pT(Vu1)a + Vp
k V 3N
+ (u Ub, (Vuj)a + Vpi)*)
2
11 lO.K
> (^iuiii,/f+ ^fii(Vui)a+v?
Hlo,/f
C%  , h rwT I (Vui)a + Vpio>A
a3/2 a
>
(^luilIk + ^jll(Vui)a + vPilli,/f) >
where we have used (4.31) for the last inequality and considered the cases
1 < 64Cf and 1 > 64Cf separately.
The stability inequality derived in Case I is essentially in the limiting
case (is = 0) of Theorem 4.2. The large meshPeclet number hypothesis (4.31)
66
is not a serious practical restriction as long as the no edge in the triangulation
is aligned with the streamline direction. If the assumption (4.31) is violated
for some I\, then the mesh is locally fine relative to the diffusion parameter u
and the use of a stabilization technique (such as bubbles) is not crucial. As a
corollary, in the piecewise linear case, the coercivity inequalities derived above
implies that our RFB method is less stable than the usual SDFEM/SUPG
method as the stabilization parameters for each method are given by ha/(3 a)
and hK/{2 a), respectively, where ha/(3a) < hK/(2 a) holds for any ele
ment. Nevertheless this apparent loss of relative stability appears to be minor.
Furthermore, this turns out to Ire an advantage for the RFB method in prob
lems where a boundary layer exists as the SUPG method produces overdiffusive
results on and around the layers: See the numerical results in [11, 15, 30] and
numerical results in section 6.
Case III: u > 0 and Quadrangular Element
A similar coercivity result, extended to the case where u > 0 and \\
is the space of piecewise bilinear functions on quadrangles I\, can be derived
in a norm which is not equivalent to the SUPG norm. We follow the lines of
the analysis given in Brezzi ct.al. [11].
For quadrangular case, from (4.9), the bilinear form can be simplified
to the following form:
67
= ^(i/(Vuh, Vui)ff + ((Vufc)a,Ui)Ar (V U!,/?!)
I<
+ (V uh,Pi))
= (^(Vu!, Vui)^ + ((Vuja, + KVu, Vui)A
K
+ ((Vu)a, ui)A + (V ug.pi))
= Y1 (^lVuilk (u> (Vui)a + ^pi)r<) ,
K
(4.33)
where we integrated by parts and used the facts that u = 0 on each dK, a
is constant, Ui = 0 on dQ, and the term (Aui,ujj) is negligible for u
with respect to the other terms. Since we are essentially in the limiting case,
we can replace u by and substitute u by its solution given as in (4.14).
Then a componentwise analysis of the term containing the bubble function in
(4.33) provides a lower bound for the bilinear term:
K(Vui)a + VpO/c
K
~ ^2 K,(viii)a + vpi)K
K
= ^ (ja (/ ((Vui)a + VPi)ds^J (Vui)a + Vpi^
68
X
((Vu[)a + Vpi)ds ) ((Vuja + V/Ji)
dK
jif lx (a'Vul +Pl,x)ds (aVi/j +pijX)
^/xX(aVM? +Pi,y)ds (a Vuf + piJ
dK
E
' K
4 lap
+
' K
4 lap
aV
aV
(aVuJ + P\,x)ds
\ 2
dK
2
(aVtii + pi,y)ds
dK
E
/
+ Pi,*)ete
jiw L [C
f [ (aVuj + Pi,y)ds
al 4X
an dT
an c?r
z{*wL[Â£ (a'v!+^
+ [ [ (aViij + pi,y)ds
4 lal 4x
an dT
a n dr
(4.34)
The final expression is positive because an > 0 at the outflow bound
ary. Thus, from the analysis, the effect of enriching the velocity space by the
69
residualfree bubbles is to add stability through these positive squared terms.
Now we summarize the results in the following theorem.
Theorem 4.3 Assume that
n ha a
for all quadrangles K in the discretization, then we have
B{uh,pi,ui,pi) > y^uili.
K
+
eU
4 lap
K
<9K+
(aVitj + pi,x)ds
+
&[ [ (aVu? + pi,y)ds
al JdK+ Jx
anr/r
a n dr
The stability results for the first two cases were obtained under the
assumption that no edge of elements in the discretization is aligned with the
streamline direction. However, this assumption seems to be only a technical
difficulty. We will present several numerical experiments showing that the RFB
method produces stable results even in the lack of this assumption.
We finally point out that the coercivitv is induced by the component
Uh = u\ + of Uh\ the component u[ plays no part here. That is, stability is
70
induced by the Lspline component Uh of the trial function uh The remaining
component u(, which appears on the righthand side of (4.11), affects only the
consistency of the method. Thus our argument does not rely on that the chosen
bubble is residualfree; it will still be valid, yet deteriorating the accuracy, if
u[ is replaced by any other choice that satisfies u[ = 0 on each OK and is
independent of Uh
71
5. Computational Algorithm
To verify the stability results obtained, we test a series of benchmark
problems in the next section. In this section we describe our computational
algorithm to solve these problems.
The solution method consists of two main bodies: One at the global
mesh level and and the other at the local mesh level. On the global mesh
level, the problem will be in the Galerkin framework while the bubble part of
the solution is still unknown. Therefore we define a submesh in each element
where we can compute the bubble functions before solving the problem on the
global mesh. This will be accomplished by using a nonstandard finite element
method. An appropriate solver is invoked to solve the problem after assembling
the bubble part of the solution into the discretization at the global level.
In the first section, we present the solution techniques we used at the
global level. Implementational issues related to the submesh will be discussed
in the second section. Note that the algorithm we describe here takes the full
nonlinear version of the NavierStokes model into account.
72
5.1 Computational Algorithm on Global Mesh
revisit the problem we wish to solve:
/ (Vu)u 2h>V e(u) + Vp = f in Q,
V u = 0 in Q, (5.1)
u = 0 on dQ,
under the same assumptions as before. The weak formulation of the problem
(5.1) reads: Find u Â£ V = Hq(Q) and p Â£ P = C(fi) D Lq(Q) such that
v,q) = F{v,q) V {v,
where
B{u,p;v,q) = ((Vu)u, v) + 2i/(e(u), e(v)) (V v,p) + (V u, q),
F(v,q) = (f,v).
We specify a finite element method by choosing the following finite dimensional
spaces on a fixed regular triangulation or quadrangulation of the domain Q:
Vh = {vÂ£(//01(f2)rdvAÂ£(F1(A')JB(/Ord,/v erh},
P\ = {peC(Q)nLl(n) IpUgP^A'), A'Â£ T},
73
where B(K) is the space of residualfree bubble functions as before. For the
implementational simplicity, we use the relation
Ve(u) =
Au +
1
2
dxi
fl(Vu)
dX2
(5.3)
(5.4)
to rewrite the system (5.1) in an equivalent form as
(Vu)u i'Au 1 Vp = f in Q,
< V u = 0 in Q, (55)
u = 0 on d(2,
by means of the continuity equation. We use the momentum equation in
(5.5) to define the residualfree bubble functions. Thus, for a given K T/,,
ip, G B(K) is zero outside K and satisfy the following system
/
(V(ux + U;,))(ui + ub) uA(ui + Uft) + Vpi = f in K,
(5.6)
ub = 0 on dK.
We note that the system (5.6) is used to find the bubble part of the solution at
the preprocessing stage which will be extensively described in the next section.
The special choice of bubble functions above (see the system (5.6))
enable us to employ the static condensation procedure. Then we infer that
74
the Galerkin method modified by residualfree bubbles reads : Find {}
\ \ x Pi such that
B{Uh,Ph',Vi,qi) = (f,vi) V{v!,9i}
(5.7)
where
B{uh,ph;vi,qi) = ((VuA)uh, Vj) + 2i/(e(ufc), e(v!))
(V VUP\) + (V uh,qi),
where u^, is given in terms of ui through the relation(5.6).
Now we concentrate on the nonlinearity emanating from the advec
tion term in (5.7). A linearization procedure employing the NewtonRaphson
method will be used to resolve the nonlinearity [28].
As the superscript is responsible for the number of the Newton itera
tions, let us write approximations for velocity and pressure part of the solution
explicitly in terms of the basis functions:
+
(5.8)
i 1,... ,nsd
(5.9)
75
where u\ and p\ represents the nodal approximation values to the solution for
the corresponding components, nen is the number of velocity nodes, nenp is
the number of pressure nodes, i/ds are the bilinear shape functions and
residualfree bubble basis functions that we will precisely define later. To be
able to linearize the discrete problem (5.7), we further decompose the nodal
values of the solution ul^n+l and p?+1 as
where u1^1 and p? are the nodal approximation values at the previous Newton
iteration step, u and p are the correction values of the solution at the current
Newton iteration step. Set Un = (un,pn) and U = (u,p). Denote the system
of all nonlinear equations we obtain when we accumulate all terms in (5.7) in
one side by G. Let U = (u,p) be a zero of G ( or equivalently a solution of the
nonlinear discrete problem) and denote by Un = (un,pn) an approximation to
the solution. By Taylors expansion around Un,
0 = G(U) G(Un) + G'{Un)U
where G' is the Jacobian matrix. The correction vector U is obtained by solving
G'{Un)U = G(Un).
(5.11)
76
Hence the Newtons method is equivalent to solving (5.11) for U and
updating the solution by
Un+l = Un + U, (5.12)
where U is guessed initially.
When the Newtons method is applied to the problem (5.7), it is
equivalent to making the following approximation
((v<+1K+1,Vl) ((v
+[((Vu^)uh,v1) + ((Vu/i)uJJ, Vl)]
obtaining the linearized form of the discrete problem (5.7) in new unknowns:
Find {uh,ph} e Vh x Ph such that
Bh{uh,Ph',vi,qi) = Fh(vi,qi) V{v1} qx] Vi x Pi, (5.13)
where
Bh(uh,ph,vuqi) = ((VuE)u/,, vj) + ((Vu^uH.Vi) + 2i/(e(u/l),e(v1))
(V vy,pi) + (V uh,qi),
Fh(vu(h) = (f,v1)((V
+ (Vv1>P?)(VuS>9l).
77
Before solving the system (5.11), we make some reasonable approx
imations for the first two terms of the bilinear form in (5.13) for the sake of
implementational simplicity:
((VuaK.vj) ~ ((Vu^u^vj),
((Vu"K,Vl) ~ ((Vu"K,Vl).
Now we write the discrete problem we solve in practice in the final
form: Find {uh,Ph} Â£ f/j x A such that
fihh(u/i, Ph;vi,<7i) = Fhh(vi, qi) V{vj, qi} Â£ Fi x , (5.14)
where
Bhh(uh,ph\vi,qi) = ((VOua.V!) + ((VuftX.vj) + 2i/(e(u/l),e(v1))
(V + (V (uh),gi),
Fhh{vuq\) = (f,Vl)((Vu1)2i/(Â£),Â£(v1))
+ (Vv1,p^)(Vu,g1)
To obtain the entries in the matrix form explicitly, we write the com
ponents of the incremental solution as a linear combination of the basis func
tions (We keep the same notation used in (5.8)(5.9) ):
78
Ui + U b
(5.15)
=
"b = J2 KVa*i
a=l,... ,nen i1,... ,n.srf
+ p* pa
a=l,...,nenp z=l,...,nsd
+
z=l,... ,nsd
ui = Y1
a=1,... ,nen i=l,... ,nsd
Pi = (516)
a=l,... ,nenp
Now we insert U&, Ui and p\ into (5.14) so that the discrete variational
problem (5.14) can be written in the following matrix form. Denoting the
number of equations by neq,
**neaxnea
X Ur,
= b
n o/i v 1
(5.17)
where A and b are given by
N\ elLin + NVelBub + K + KVB : NPreBub + KPB G
A =
Gl + YelBub : PreBub
(5.18)
79
[F] [il/] u" [N} pn
(5.19)
[g] ura [R]pn
where each submatrix consist of the contribution matrices where explicit
statements of their entries can be found in Table 5.1:
M = [iVrZ/1] + [M'Bl] + [A'] + [K\'B]
N = [NPD1] + [KPB] [G]
Q = [G^ + lVelBub]
P = [PreDub\
Note that the matrices NVelBub, KVB, NPreBub, KPB, V'elBub,
PreBub, NVB1 and NPB1 are absent in the Galerkin formulation of the
problem when piecewise linears are used only. These matrices are actually
contributed to the formulation by the residualfree bubble method. These
contributions are responsible for the stability of the scheme.
Following the assembly phase to form the global linear system of
equations, we solve the system by invoking an iterative solver or a direct solver.
We choose the GMRES of Saad and Schultz [66] as an iterative method which
is able to solve the nonsymmetric systems. The number of Krylov vectors, the
number of passes in GMRES and a tolerance of how well one solves each linear
80
(Arl elLin)pq = + ((Vu^e,,^)
{NY el Dub) pq = ((V<^ej)u7,^ei) + ((Vuf)^,^)
(F )pq = 2j/(e(^e,),e(^*ei))
(KVB)pq = 2i/(e(v^ej),e(^el))
(NPreDub)ph = T/j = h...,nsdK(^lHgj)U^'lla^) +
(KPB)* = Ej=i,...,ns),e(V'iel))
(Gpi) = (^,V ($&))
(Glq) = (V(^e_j),^a)
[V elDub)aq = (V(^ej),^a)
(PreBub)ai ~ ^J=l,...,nsd(^7 '
(F)P = (/, $&)
(NVLl)pq = ((V^e^u^^e,)
81
(Nvm)pq = ((V^eQu",^)
(VPBl)pb
p = nsd (a. 1) + i i Â£ {1,... nsd}
+ 1 CC Â£ II J Â£ {!, nsd}
Table 5.1. Entries of the contribution matrices
system, are chosen by the user. GMRES is preconditioned by the diagonal. A
skyline direct solver which makes efficient use of the memory is used otherwise.
At the top of the specified solver, we put the NewtonRaphson loop.
In practice, the nonlinearity is solved by a combination of the Newtons and
QuasiNewtons Methods. As a result of this combination of methods, we point
out the changes in the implementation as follows : The Main Matrix is saved
in element blocks and computed in the first iteration of this loop, and only
when necessary. In other words, we compute the residual at every iteration
and if after solving the linear system the solution produces a smaller residual,
then the same matrix is kept and a new iteration is performed with the same
matrix (a quasiNewton iteration, rather than a Newton iteration). If, however,
the residual is increased after a QuasiNewton step, then the QuasiNewton
solution is dropped and a Newton step is performed (computing the matrices
82
in element, blocks). This procedure avoids the computation of an updated
matrix at each iteration. This NewtonQuasiNewton loop is finished if either
a user specified number of iterations is reached or if the tolerance, also specified
by the user, is reached.
5.2 Computational Algorithm on Submesh
To compute the contribution of residualfree bubbles to the problem
on the global mesh, we need to calculate the bubble functions in each element.
This can be done by solving (5.6). However, the equation (5.6) is not linear due
to the presence of the nonlinear advection term. In global problem, we have
dealt with nonlinearity by a linearization process, however, we approximate
here the nonlinear term (VuÂ£+1)uÂ£+1 by (VuÂ£+1)u" Then the system (5.6)
can be written as
/
(VuÂ£+1)u? ^Au"+1 = f (Vu"+1)u" + !/Au"+i Vp7l+l in K ,
 ufe = 0 on dK,
(5.20)
where the superscript again stands for the number of Newton iteration. The
equation (5.20) can be viewed as a scalar advectiondiffusion equation in each
component.
Finding the solution of the problem (5.20) may be as hard as solving
the strong form of the original problem. Therefore, we specify another finite
83
element method for (5.20) in the interior of each element, approximate the
bubble functions by this method and then use these approximations in the
problem at global mesh level in place of the exact solutions of the bubble
functions.
Thus, to find u"+1 which is dependent on u+1 through (5.20), we
solve for the bubble shape functions y^s, Pas an(l Vfi wifh a varying from
one to the number of element velocity nodes (nen), a varying from one to the
number of element pressure nodes (nenp) and i varying from 1 to nsd:
where L = u V vA, the ipla's and the z/^s are the local linear or bilinear
basis functions for Uj and p\, respectively, and
(5.21)
(5.22)
/
Lp)et = /*e, in AT,
(5.23)
\
Thus, if
84
p?+1 = Ert+va,
a
then
<+I = E E <'"+1rts.
a=l,... ,nen i= 1,... ,nsd
+ E "*+l E ^
a=l,nenp z=l,...,nsrf
+ Y1 P/*
2=1,... ,725d
and substituting into (5.14) we get the matrix formulation (5.17).
Now, consider a submesh (a mesh defined for each element) where
we solve the equations (5.21)(5.23) by a suitable nonstandard method finite
element method, the Galerkinleast squares method (GLS), for example. We
will need to repeat this procedure for each element and for each basis function.
Let us now describe the application of this strategy to the vectorial
advectiondiffusion problem (5.20) using the GLS method. Thus corresponding
to (5.215.23) we need to solve in each element K of the original mesh:
 u Vv^e, vAip^e, = u? + uAa = 1,.., nen in K
iplaei = 0 on dK
/
u V (p\ei uA (p^e, = f^e* a = 1,.., nenp in K
y iflei = 0 on OK
85
a = 1,nen
in I\
I
u V(flfet i/A(plfei = ftet
95ye, = 0 on DI\
To construct the GLS approximations of these differential problems,
let us denote by I\* an arbitrary element in the submesh with diameter h* and
by ^* the basis function for a piecewise linear interpolation in the submesh,
as in section (3.4), so that our unknown bubble basis functions ^s, <^~s and
Here l runs over all unknown interior nodes in the submesh, say through N*.
The index (a,i) and (a,i) refers to the specific bubble function that we are
trying to compute. (Recall that we need to solve for a number of bubbles
equals to the number of basis functions used to span the piecewise polynomial
defined in the global element A .)
The matrix form of the problem obtained through the GLS method
reads: For each a (from 1 to nen), for each a (from 1 to nenp) and for i =
(5.24)
86
1,... nsd, find l = 1,2,N* and c\a'l\ / = 1.2,V', such that
E^IK + ('Wl.vft)
l
+ W V* V&Pt, r  nAC,))]
= (< ' C + r(J Vt/>,* tAfi))
Â£>!a,i)[(ur VrPIPJ + W.Vft)
l
+W v$" vArf, t( u?  ^a?o]
dlpa,K
= (
<5x,
fC + r  ^a^o,
for in = 1, 2,IV*. We use the stability parameter r given by [28] as:
T(x, Pe^.(x)) =
h*
2u"(x)
Â£(PeA.(x)),
(5.25)
PeK*(x)
6u
(5.26)
((PeA.(x)) =
1
PeA'.(x), if 0 < PeA.(x) < 1,
1.0, if PeA.(x) > 1,
(5.27)
nsd
1/2
K(*)2= > u^(x)
i=l
(5.28)
We similarly account for the basis functions associated with f, tpj
87
Find cpl\ l = 1,2, ..,N* such that
for i G {1,nsd}.
them in (5.24) to get the approximate bubble basis functions pl^h > Pp Pp
These approximations are then plugged into the global problem with the help
of the representing identities in (5.16) in place of the exact residual functions
p^s, cp~s and p^s. The global problem may now be solved for the piecewise
linear part of the solution. We can also recover the bubble part of the solution
approximately, restricted to a global element K, by
+
(5.29)
i=l,... ,nsd
with p1^' in lieu of pla and pp in lieu of p\.
88

