Citation
A mixed finite volume element method for accurate computation of fluid velocities in porous media

Material Information

Title:
A mixed finite volume element method for accurate computation of fluid velocities in porous media
Creator:
Jones, Jim E
Publication Date:
Language:
English
Physical Description:
vii, 101 leaves : illustrations ; 29 cm

Subjects

Subjects / Keywords:
Finite dynamics ( lcsh )
Finite element method ( lcsh )
Multigrid methods (Numerical analysis) ( lcsh )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

Notes

Bibliography:
Includes bibliographical references (leaves 100-101).
General Note:
Submitted in partial fulfillment of the requirements for the degree, Doctor of Philosophy, Applied Mathematics.
General Note:
Department of Mathematical and Statistical Sciences
Statement of Responsibility:
by Jim E. Jones.

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:
35327241 ( OCLC )
ocm35327241
Classification:
LD1190.L622 1995d .J66 ( lcc )

Downloads

This item has the following downloads:


Full Text
A MIXED FINITE VOLUME ELEMENT METHOD
FOR ACCURATE COMPUTATION OF FLUID
VELOCITIES IN POROUS MEDIA
by
Jim E. Jones
B.S., New Mexico Tech, 1986
M.S., New Mexico Tech, 1989
A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
Applied Mathematics
1995


This thesis for the Doctor of Philosophy
degree by
Jim E. Jones
has been approved
Stephen F. McCormick
Thomas A. Manteuffel
John W. Ruge
Thomas F. Russell
John Trapp
Date


Jones, Jim E. (Ph.D., Applied Mathematics)
A Mixed Finite Volume Element Method for Accurate Computation of
Fluid Velocities in Porous Media
Thesis directed by Professor Stephen F. McCormick
ABSTRACT
A key ingredient in the simulation of flow in porous media is the accurate
determination of the velocities that drive the flow. Large-scale irregularities
of the geology, such as faults, fractures, and layers, suggest the use of irreg-
ular grids in the simulation. In this study, the approach to this problem was
to apply the finite volume element methodology, developed by McCormick,
in conjunction with mixed methods, developed by Raviart and Thomas. The
resulting mixed finite volume element discretization scheme developed here
can be applied in a clear and straightforward way to irregular grids and
is appealing because of its local conservation properties and its direct and
accurate representation of physical intercell flux terms. Several multilevel
algorithms are developed that provide efficient methods for solving the set of
equations that this discretization produces. This thesis includes numerical
results from a variety of test problems, from Poissons equation to prob-
lems with anisotropic, discontinuous, and tensor diffusion coefficients. These


results show that this approach has the potential to generate accurate ap-
proximate fluid velocities and that the multilevel methods can provide fast
solvers.
This abstract accurately represents the content of the
candidates thesis. I recommend its publication.
Signed
Stephen F. McCormick
IV


Contents
Chapter
1 Uniform Grids .................................................1
1.1 The Mixed Finite Volume Element Discretization............1
1.2 An FVE Based Multigrid Algorithm .......................13
1.2 Computational Results ...................................25
1.2.1 Comparison to Finite Differences for Poisson Problem ... 25
1.2.2 Discontinuous and Anisotropic Diffusion Coefficients .... 30
1.2.3 An Alternative Multilevel Algorithm..............35
1.2.4 Full Tensor Diffusion Coefficients ..............42
Chapter
2 General Quadrilateral Grids..................................46
2.1 The Mixed FVE Discretization ............................46
2.2 Multilevel Solvers ......................................55
2.2.1 FVE Based Multigrid..............................55
2.2.2 A Two-Level Approach ............................60
2.3 Computational Results ...................................62
2.3.1 Multigrid Performance............................62
v


2.3.2 Two-Level Performance................................67
2.3.3 Accuracy of the Mixed FVE Discretization ............70
Chapter
3 Technical Derivations ...........................................77
3.1 Uniform Grids and Tensor Permeability.......................77
3.2 Derivation of FVE Stencils for General Quadrilaterals .....80
3.2.1 Scalar Diffusion Coefficient ....................80
3.2.2 Tensor Diffusion Coefficient ....................89
3.3 Using Fine Grid Integrals to Calculate Coarse Grid Integrals .... 91
Chapter
4 Summary and Future Work .........................................97
Bibliography.......................................................100
vi


ACKNOWLEDGEMENTS
This thesis presents the results from a five-year period of research, and
it would not be complete without a word of special thanks to those whose
assistance enabled me to complete it.
First of all, I wish to thank Professor Steve McCormick for giving me the
opportunity to work on a host of interesting research projects. I also thank
Professor McCormick for the wide freedom that I was given as a graduate
student to do independent research, and his patience as I labored to bring
this research to its conclusion.
I also wish acknowledge Professor Thomas Russell for his assistance with
this work. Without his expertise in the held of reservoir simulation, this
research could not have taken place. I also wish to thank Professor John
Ruge and Professor Steve Schaffer. While they were not directly involved in
the day-to-day work on this project, I learned most of what I know about
multigrid methods through collaborative research with these two men.
I would also like to thank Dr. Pat Roache, Ecodynamics Research, and
Sandia National Laboratory for their financial support at various times dur-
ing this research.
Finally, I wish to acknowledge the countless sacrifices made by my family
- my wife Dorothee, my daughter Frederika, and my yet un-named second
child that allowed me to complete this thesis. Most of all I thank my wife
for her unfailing love and support during this time.
vii


1 Uniform Grids
1.1 The Mixed Finite Volume Element Dis-
cretization
In this first section, we introduce the new mixed finite volume element dis-
cretization technique in a familiar context: solving the diffusion equation on
a uniform, square grid. However, it should be clear from the presentation
that the discretization developed in this section, and the multigrid algorithm
constructed in the next section, can easily be modified to apply to any tensor
product grid.
We begin by considering the following partial differential equation defined
on some square coordinate-aligned domain 0 in 1Z2:
f- V A(x)V(x) =/(x), X E o,
[ A(x)V(x) rj = g(x), x G <90. ' '
In the context of reservoir simulation, this is the pressure equation for in-
compressible single-phase flow, where is the pressure in the reservoir 0
and the boundary condition specifies the flux on <90. Because one of our
1


goals for the new discretization is accurate approximations of flow velocities,
we will begin by reformulating this equation as a first-order system of equa-
tions, where velocity appears explicitly. This is done by introducing the flow
velocity variables via the definition
v = -AV(f>, (1.2)
then rewriting the partial differential equation in (1.1) as
V v = /. (1.3)
In the context of reservoir simulation, definition (1.2) is Darcys law and
equation (1.3) is the mass conservation law. In reservoir simulation, this
same approach of treating flow velocity explicitly has been used in mixed
Unite-element methods with considerable success [11],[12],[18].
To continue the derivation of the discretization, assume that A is a diag-
onal tensor,
Here, ax and ay may be functions of position. Section 3.1, contains the details
of the discretization process for the case where A is a full tensor. Let u and
v denote the components of the velocity,
(it, vf = v.
Then the second-order partial differential equation (1.1) may be written in
2


the form of a first-order system in 0:
ax~1u + d(f)jdx = 0 (u equation)
< ay~xv + d(f)/dy = 0 (v equation)
du/dx + dv/dy = / (p equation),
with boundary conditions,
u(x) fi'(x)5 x ^^west U ^^east
u(x) H(x)*> x ^^south ^ '^nortlr
The labels u,v, and p for the three equations in equation (1.4) are introduced
simply for convenience.
To discretize this system, we follow the finite volume element (FVE)
principles developed in [7],[13],[14], The two major components of any FVE
discretization scheme are the choice of control volumes over which to inte-
grate the continuous equation and the choice of finite element spaces for the
unknowns.
First, to introduce the control volumes for the FVE discretization, con-
sider a uniform square mesh Qh with mesh size h that covers 0. We introduce
three sets of control volumes, one for each of the three equations in (1.4).
Examples of the control volumes are shown in Figure 1.1. We denote by hi
the set of all volumes U that will be used to discretize the u equation in (1.4).
Similarly, we use the notation V and V for the sets of volumes V and P for
the v and p equations, respectively. Note that the control volumes P corre-
spond to the grid blocks defined by Qh, the control volumes U straddle the
vertical grid edges of Qh, and the control volumes V straddle the horizontal
grid edges of Qh.
3


Figure 1.1: Uniform Grid
4


Next we will introduce the finite elements for the FVE discretization. For
our finite element space, we use the lowest order Raviart-Thomas elements
[16],[19] on the the control volumes V:
uh is linear in x and constant in y on each P £ P,
vh is linear in y and constant in x on each P £ P,
Note that uh is continuous across vertical edges of the grid and that vh is
continuous across horizontal edges of the grid. The location of the nodes for
each of the unknowns with their indexing is also shown in Figure 1.1. The
flux boundary condition in equation (1.1) specifies the values of uh on the
west and east boundaries of 0 and the values of vh on the north and south
boundaries of 0. Noting some of the characteristics of the finite element
spaces, the pressure space contains all functions that are piecewise constant
on Qh. A typical basis function for uh is shown in Figure 1.2.
Note that the normal component of this basis function is continuous across
all grid edges, is zero on all edges except the edge shared by the two P
volumes which are the support of the function where the normal component
equals a nonzero constant. With these properties we can guarantee that the
computed flow velocity will also have continuous normal components across
grid edges. This will also be the case for general quadrilateral grids as we
will see in Section 2.1. The true physical solution also has this property -
continuous normal component of velocities but not every numerical scheme
for approximating it does, as [17] pointed out.
5


Figure 1.2: u Basis Function
We are now ready to discretize the equations. Start by taking the u
equation in (1.4) and integrating it over each U £ U. This will result in a
discrete version of Darcys law. As an example, let U8J be the volume in hi
that is centered at the interior uh node (i,j), as shown in Figure 1.3.
We then have
/ (ax ru + dcfi/dx^j dxdy = 0. (1.5)
Consider the two terms in the integral separately. For the first term in
equation (1.5), we split the integral into two parts, one for each half of U4-j.
This yields,
/ (ax dxdy = / (cix 1u^ dxdy + / (cix 1u^ dxdy, (1.6)
t/ X-J 7 7 t/ X-T ^ t/ X-T ^
6


Figure 1.3: UVolume
7


where U* and U[? denote the left and right halves of the control volume.

Now because the finite element approximation uh is linear in x and constant
in y on each of the halves, by direct calculation it follows that
J t dxdy = a^1^)
A similar expression is obtained for the integral over the right half of U4-j.
Note that here we have assumed a^-1 to equal the constant value a:E_1(/) on
the left half of the control volume. If it is instead piecewise constant on some
finer mesh, as will often be the case in the section on computational results,
then the integral can be split into parts on this finer mesh and evaluated.
The details of this piecewise integration can be found in Section 3.3, where
we show, in the more general case of quadrilateral grids, how to use integrals
calculated on fine grids to generate integrals on coarser grids. Now for the
second term in equation (1.5), we have
' U.
(dcf)/dx) dxdy = / / (d/dx) dx dy.
Then apply the divergence theorem to get
f (f [d(h,y) (0,y)) dy.
(1.7)
Then, as the finite element approximation x = 0 and x = h, it follows that
((h,y) ~ (0,y))dy = h
hJ
,j
(1.8)
8


Putting these observations together, the discrete equation that equation (1.5)
gives rise to is thus
Y \ax 1(?) (3+ Ui+itj) + ax 1(l) (3ufj +
{^,3 ~ $-1,3) = O'
If ax~x is constant, the equation simplifies to
-1
X
uUi + K,
"I- u
*+i j
Note that this expression is quite similar to the expression one would get
using standard cell-centered finite differences to approximate Darcys law:
hi
(1.9)
The difference between the two expressions is that, with the mixed FVE
method, tridiagonal mass matrix multiplies the velocities rather than a triv-
ially invertible diagonal matrix. One way of viewing equation (1.9), is as
a midpoint rule approximation to the integral in equation (1.6). In Sec-
tion 1.3, we will show that this can result in less accurate approximations for
the velocities. More importantly, as we will see when we consider irregular
grids, the mixed FVE method has the property of generating discrete ver-
sions of Darcys law for geometries where it is unclear how to use standard
cell-centered finite differences.
Integrating the v equation in (1.4) over an interior V volume yields a
similar discrete expression involving nodal values of vh and 9


To complete the discretization process, integrate the p equation in equa-
tion 1.4 over the volume in V centered at the interior denote this volume by P4-j. See Figure 1.4.
Integration is simply a matter of applying the divergence theorem. From
the definition of u and u, we have
/ (du/dx + dv/dy) dxdy = / V-vdxdy.
Transforming the volume integral to a surface integral yields the following:
fPij V vdxdy = fap.'. v -yds =
= fgpe udy fdpw udy + fgpn vdx fgps vdx}
where dPp- denotes the east edge of the control volume, and similarly for the
west, north, and south edges. Then, since the finite element approximations
uh and vh are constant on edges, we have
fgpe udy fgpw udy + fgpn vdx fgps vdx
h(u
V.J l,J l,J
h o,h h h \
t + l,3 U+3 ^ Fj + 1 Fjh
The discrete version of the p equation in (1.4) is therefore just
h(ui+i,j ~ uh + vti+1 vt) = h2kr
Here, on the right hand side, we have assumed that / is a piecewise constant
function on V, with /8J denoting the value of / at the way of looking at this is: we have used the midpoint approximation to the
integral; if we had assumed that we had more information about /, we could
have used a more accurate approximation for its integral. For example,
10


Figure 1.4: P Volume
11


knowing / as a function of x and y, we could have used a Gauss quadrature
rule to approximate the integral.
In summary, the discretization has produced discrete versions of Darcys
law for each volume U £ hi and each volume V £ V. These equations relate
the pressure drop between cells to a linear combination of velocities. The dis-
cretization has also produced discrete conservation equations for each volume
P G V. For each P volume, we have that the discrete flux out of (or into) the
volume balances with the strength of the source (or sink) in the volume. This
conservative property is an important feature of FVE discretization schemes
in general. The mathematical modeling of many physical systems (including
flow in porous media) starts with a conservation law stated in integral form.
In going from the continuous mathematical model to the discrete mathe-
matical model, preserving some analogy of this conservation law is desirable.
For time dependent applications, this discrete conservation property is often
necessary to generate acceptable solutions.
The discretization technique described in this section can easily be gener-
alized to handle problems in three dimensions. Also, it can be easily modified
to handle the Dirichlet problem, where the flux boundary condition in equa-
tion (1.1) is replaced by
In [10] this mixed FVE discretization is applied to Poissons equation with
Dirichlet boundary condition on a three dimensional domain. This reference
12


contains the details of the necessary modifications.
1.2 An FVE Based Multigrid Algorithm
We begin with a very brief overview of multigrid fundamentals, assuming that
the reader is already somewhat familiar with the held. Good references are
the introductory tutorial [6] and the advanced guide [3]. Another useful
reference is [14], which presents multigrid (and multilevel adaptive) methods
for use in conjunction with FVE discretization techniques.
Consider the problem of finding the solution to the discrete set of equa-
tions
Lhuh = fh. (1.10)
We think of equation (1.10) as coming from the discretization of some contin-
uous partial differential equation on a uniform grid with mesh size h, although
multigrid methods can be applied to discrete problems with no continuous
origin (e.g., [4]). If the number of unknowns is too large to solve the equation
directly, one might try simple Gauss-Seidel relaxation to find an approximate
solution, uh, to (1.10). One typically finds that the convergence speed of this
process degrades very early in the iterations: the first few relaxation sweeps
may substantially reduce the residual
h rh t h A h
r = / L u ,
but, subsequent relaxation sweeps do not. The problem is that relaxation
methods generally do not reduce low-frequency errors very well. In multigrid
13


methods, coarser grids are used to solve for and then eliminate this low-
frequency error that relaxation alone cannot resolve. This is done by way of
the error equation
Lheh = r\ (1.11)
where eh denotes the error in the current approximation:
eh = uh -uh.
Because relaxation does work well for high-frequency error, after a few relax-
ation sweeps one can represent the error accurately on a coarser grid. With
this in mind, we consider a coarse grid version of equation (1.11):
L2hu2h = I2hhrh. (1.12)
Here we think of the operator L2h as coming from the discretization of the
same continuous partial differential equation on a uniform grid, but with
twice the mesh size as the original problem (1.10). The operator I2h is called
the restriction operator and is used to transfer the fine grid function rh to
the coarser grid. One may be able to solve equation (1.12) exactly for u2h,
remembering that this system is on a coarser grid so it is not as large as the
original one. Then one would correct the fine grid approximation by u2h,
which is approximately equal to the error in uh, according to
uh <- uh + I^hu2h. (1.13)
Here, the operator 7^, called the interpolation operator, is used to transfer
the coarse grid function u2h to the fine grid. We now perform a few relaxation
14


sweeps to reduce any high-frequency error that we may have introduced. This
completes one multigrid cycle, the steps of which were:
1) Relax the fine grid problem.
2) Form the coarse grid version of the error equation.
3) Solve the coarse grid version of the error equation.
4) Correct the fine grid approximate solution.
5) Relax the fine grid problem.
In practice, one may not be able to solve the coarse grid version of the error
equation for the same reasons one cannot solve the original problem: it is
still too large, so that relaxation alone does not work and direct methods
are too expensive. One then can use the above multigrid cycle recursively on
this problem, using a still coarser grid with mesh size Ah. This process can
be continued until a coarse enough grid is reached, where a direct method or
relaxation by itself is effective. One multigrid cycle is represented graphically
in Figure 1.5.
The basic components of a multigrid algorithm are the grids themselves,
the discrete operators on each grid (Lh, L2h, L4h,. .), the relaxation process,
and the grid transfer operators (interpolation and restriction). The multigrid
method is an iterative method, but typically one needs only a few cycles to
produce a good approximation. Moreover, the so-called full multigrid algo-
rithms, which start with multigrid cycles applied to coarser levels, can achieve
approximations with accuracy comparable to the exact discrete solution in
the equivalent of just two or three fine grid cycles.
Now to develop a multigrid algorithm designed specifically for solving
15


Figure 1.5: One Multigrid Cycle
16


the discrete set of equations generated by the mixed FVE discretization of
the previous section, consider a family of uniform square grids Qh of mesh
size h that cover the region 0. Figure 1.6 shows a coarse grid with
twice the mesh size of the grid Qh in Figure 1.1. Note that the coarse grid
looks like the fine grid. Simply joining four square fine grid cells produces
a square coarse grid cell (for the general quadrilateral grids in Section 2.2
producing coarse grids from fine grids is not so simple). On each grid Qh, we
can apply the mixed FVE discretization process, and we write the discrete
set of equations that this process generates as
Lhzh = gh} (1.14)
where zh = (uh,vh,hy and gh = with the unknowns uh, vh,
and on
grid Qh. Note that, in the FVE discretization of the previous section,
and were zero. In the multigrid algorithm, on coarse grids these will,
in general, be non-zero when we form the coarse grid version of the error
equation.
We now must define interpolation and restriction operators and a relax-
ation process. For defining interpolation operators, we use the same general
principles as outlined in [14], but we must be specific here. The finite vol-
ume element discretization is based on finite element spaces for the variables
uh,vh, and ment spaces of the different grids to define interpolation. In particular, the
17


Figure 1.6: Uniform Coarse Grid
18


finite element spaces are nested in that one can write a coarse grid basis
function as a linear combination of fine grid basis functions. This nesting
property simplifies the definition of interpolation operators. To define the
interpolation operator for , note that (f)2h is constant on the grid 2h
volume P/,j. Referring to Figure 1.7, we thus have the following characteri-
zation of h = I{(f))^h(f)2h:
bh-
hj
ih
^'+1 ,j
ih
^'j + 1
ih

(1.15)
To define the interpolation operator for it, /(n)^, note that u2h is linear in
x and constant in y on the grid 2h volume P/,j. We thus have the following
characterization of uh = (see Figure 1.7):
ui,j = ui,j+i =
Ui+2,] = Ui+2J + 1 =
(1.16)
i+lj Ui+lj + 1 1/2(mJ+1,J + M/,j)
Dehnition of the interpolation operator for v is similar.
For defining restriction operators, we again use the same general prin-
ciples as outlined in [14]. In the multigrid algorithm, restriction operators
are used to transfer right hand sides and residuals of equations, not the un-
knowns themselves. The definitions of the restriction operators are therefore
naturally based on the relationship between the volumes of the various grids.
The idea is to lump several of the grid h right hand sides to produce the grid
2h right hand sides. To define the restriction operator for the P equation,
J(P)2/l, note that a grid 2h volume P/,j consists of four of the P volumes
19


Figure 1.7: Fine and Coarse Grids, I
20


of grid h. We thus have the following characterization of f2h = I(P)2hhfh,
referring again to Figure 1.7:
r2h rh rh
JI,J Ji,j '
+ fh+1
+ /i+l,j+l
(1.17)
To define the restriction operator for the U equation, I(U)2hh, note that a
grid 2h volume Uyj of Q2h contains all of two of the U volumes of grid
h and half of four others. We thus have the following characterization of
f2h = I(U)2hhfu, referring to Figure 1.8:
1
flkI,J ~ fui,j + fui,j +1 + + + l + fui+l,j + fut + ij + i)- l1-18)
Dehnition of the restriction operator for the V equation is done in a similar
fashion.
For the equations on grid Qh, there are a variety of relaxation techniques
one could use. We will present the two that have proved to be most effi-
cient in our numerical experiments: distributive Gauss-Seidel relaxation and
box Gauss-Seidel relaxation. We consider these two relaxation schemes be-
cause there is no one-to-one correspondence between equations and unknowns
(point Gauss-Seidel makes no sense per se). There are three types of discrete
equations:
U equations:
\ax1 + 6ui,j + Ui+l,j) + ^
V equations:
1 (Gj-i + vi,j + Gj'+i) + h
P equations:
h(ui+i,j - + uG+1 VP) = h2fitj.
kh
h-1 j
0-
M-l
= fh. .
J U % J 5
fh. .
J V % J 5
(1,19)
21


Figure 1.3: Fine and Coarse Grids, II
22


Again, we point out that the labels are used only for convenience; for example,
there is no inherent reason to associate the P equation closely with the
variable and, in fact, it does not even appear in this equation.
The guide [3] contains an excellent description of distributive Gauss-
Seidel relaxation schemes, so that here we will be content with a brief de-
scription of our scheme. We can naturally associate the U equations with
the velocity variable uG of the grid edge that the U volume straddles. Like-
wise, we can naturally associate the V equations with the velocity variable
vG of the grid edge that the V volume straddles. We therefore relax these
equations using standard point Gauss-Seidel relaxation. The basic idea is
to relax the U, V, and P equations together locally in such a way that the
error is smoothed while the U and V equations are not damaged, i.e., their
residuals are essentially unchanged. We can think of the overall relaxation
process (including point Gauss-Seidel for U and V) as a three step scheme.
First, we sweep over all of the uh nodes, change the value of in turn
so that the U equation at (i,j) is satisfied. Second, we perform a similar
Gauss-Seidel relaxation of all of the V equations. Note that these two steps,
U and V relaxation, are independent of each other and could be performed
in parallel. In fact, for full independence, point Jacobi could be used within
each variable without much loss of efficient smoothing because of the associ-
ated tridiagonal mass matrix. Finally, we sweep over each simultaneously changing the value of G and the values of uh and vh that
23


lie on the edge of the volume P8J (namely, , rt^+1 and uG+1) so that
the P equation at (i,j) is satished and so that the residuals of the U equa-
tions at (i,j) and (i + 1, j) and of the V equations at (i,j) and (z, j + 1) are
unchanged. To allow for effective vectorization, the Gauss-Seidel relaxation
process in each of these three steps is done using a red/black ordering.
The second technique we consider is box relaxation. This approach in-
volves only one step, which consists of a sweep over each P volumes in turn
changing the values of the associated five unknowns -,u^-,u^+1 ,and
vf -+i so that the corresponding five equations are satished. We will often
need to use a y-line box Gauss-Seidel relaxation technique, which is the same
as box relaxation except that all of the P volumes sharing the same z index
are relaxed simultaneously as a group. Analogously, x-line box Gauss-Seidel
involves relaxing all P volumes that share the same j index as a group. We
will use the term xy or alternating line relaxation to signify a single relax-
ation step where we perform x-line relaxation followed by y-line relaxation.
As will be evident from the results of numerical experiments the best choice
of relaxation technique, in the sense of the most efficient solver, is problem
dependent. The line relaxation techniques are clearly computationally more
expensive; but for problems with anisotropic coefficients, line relaxation is
often needed to achieve acceptable convergence rates for the multigrid algo-
rithm.
24


1.3 Computational Results
The major motivation in developing the mixed FVE discretization approach
was to have an accurate discretization for irregular grids, where there is no
clear way to use standard cell-centered finite differences. However, a case
can be made for using the mixed FVE discretization even on uniform grids,
as we will demonstrate in this section.
1.3.1 Comparison to Finite Differences for Poisson
Problem
We begin with the simplest of problems of the type to which our discretization
rpplies, namely Poissons equation on the square domain 0 = [1, l]z
f v v(x) = /(x) x e o,
) V(x) r] = g(x) x G <90
(1,20)
One can write this equation as a system of first-order equations in the way
presented in Section 1.1, then apply the mixed FVE discretization. We
will call this the mixed FVE approach. Alternatively, one could discretize
equation (1.20) using standard cell-centered finite differences, then solve the
resulting discrete set of equations for an approximation to . If we were
interested in approximations to velocities or fluxes, we could then apply ap-
propriate standard finite differencing to the approximation of . We will call
this the finite difference approach.
To compare these two approaches, consider the function
(x, y) = cos(kiTrx) + cos(k2Try)}
(1,21)
25


where k\ and k2 are integers. If we set
f(x, y) = 7T2(kf + ^2)(cos(A;i7ra;) + cos(k2Try))}
g(x,y) = 0
in equation (1.20), then the function in (1.21) is a solution to the partial
differential equation. It is not the only solution because we can add an
arbitrary constant to the solution , while preserving the property of satis-
fying the partial differential equation and boundary condition. The function
in (1.21) is the unique solution that satisfies the additional condition that
its integral over 0 is zero. By varying k\ and k2} one can see the effect that
oscillations in the true solution have on the accuracy of the numerical solu-
tions generated by a mixed FVE approach and a finite difference approach.
We measure the error in the solutions in the following way. For the mixed
FVE approach, we calculate the C2 norm of (f)m^ve and call this quantity
^mfve' gecause the mixed FVE discretization relies on finite element spaces,
one has an approximation to (j) everywhere in 0 except at grid interfaces, so
the C2 norm make sense. For the straightforward finite difference approach,
one has approximation to ^ only at the nodes. To create a functional approx-
imation to , we assume that the nodal approximation holds for the entire
cell, then calculate the C2 norm of its difference with the exact solution. We
call this quantity qBecause we are interested in accurate approximation
of velocities, which are here just the components of the gradient of , we
measure the error in them as well, which we do in the following way. For the
26


mixed FVE approach, we calculate the difference between the fluxes,
s
using the exact and the approximate solutions over each grid edge. We
then take the £2 norm of the difference on vertical edges and call this u]£^ve.
Analogously, we take the £2 norm of the difference on horizontal edges and
call this vfve. Note that if the exact velocity is V and the mixed FVE
approximation is v, then ((uenrl^ve)2 + (vem^ve)2)^ is equivalent to a discrete
H(div) norm of the vector velocity error (which incorporates C2 norms of
(v V)^, (v V)^, and div(v V), the last of which is zero by the local
conservation property of the mixed FVE method). For the finite difference
approach, we do the same to calculate v/ed and v£d. To get approximations
for v, we use standard finite differencing of the approximation of .
Tests were run on a grid with 64 cells in each direction and the results
are shown below.
Mixed FVE
h k2 (j/n fve urnfVe Vmfve ue
l 1 4.01E-2 2.52E-6 2.52E-6
l 4 1.17E-1 1.74E-2 4.38E-3
l 8 2.27E-1 7.14E-2 9.16E-3
l 12 3.37E-1 1.44E-1 1.28E-2
l 16 4.41E-1 2.10E-1 1.49E-2
15 13 5.27E-1 2.59E-1 2.92E-1
27


Finite Difference
h k2 <
l 1 4.01E-2 3.89E-12 3.89E-12
l 4 1.17E-1 3.56E-2 8.96E-3
l 8 2.32E-1 1.57E-1 2.02E-2
l 12 3.53E-1 3.62E-1 3.20E-2
l 16 4.84E-1 6.55E-1 4.51E-2
15 13 5.46E-1 8.79E-1 9.92E-1
One can see that both methods produce approximations to of roughly
equal quality. When the solution is smooth, the approximations to the fluxes
are also of equal quality. However, for oscillatory solutions, i.e., increasing
&2 (k2 = 16 is the most oscillatory function that can be represented on this
grid), the mixed FVE obtains somewhat more accurate approximations to
the fluxes. This is also true for the last case where the solution oscillates
in both directions. This is not surprising considering that one can view the
finite difference approach as a midpoint rule approximation of the integral
in the Darcy equation, as mentioned in Section 1.1. This approximation is
accurate for smooth functions, but not so accurate for more oscillatory ones.
To be fair in comparing these two approaches, one should also consider
the work necessary to generate the approximations. For solving the equations
in the finite difference approach, we used black box multigrid [8]. This code
is aimed at solving much more difficult problems than Poissons equation, as
we will say more about later, and is therefore not the most efficient code for
such a simple problem. For solving the mixed FVE equations, we used the
multigrid algorithm of the previous section with distributive relaxation. Both
28


multigrid methods proved to have similar worst-case convergence factors,
approximately .1, so that the amount of work for solving the mixed FVE
equations would be about three times that for the finite difference equations.
Loosely speaking, the factor of three comes from the fact that, in the mixed
FVE approach, one has roughly three times the number of equations as in
the finite difference approach. A careful comparison of relative costs of the
two algorithms should compare relaxation complexity: while the mixed FVE
scheme is more involved, its individual stencils are a little smaller, so that
the factor of three is still approximately correct. If one were to compare
fully optimized codes, we are confident in our belief that the mixed FVE
solver would be no worse than three times slower than the finite difference
solver. For oscillatory solutions we have seen that the error in the mixed FVE
fluxes can be up to three times smaller than the error for the finite difference
approach. So one might argue that the comparison comes out very roughly
even. However, these tests certainly show that the mixed FVE approach has
the potential to generate more accurate approximations for fluxes than the
standard finite difference approach even for simple problems, and for slightly
more complex problems, the difference in accuracy between the two methods
can be much greater as will be shown in the next section.
29


1.3.2 Discontinuous and Anisotropic Diffusion Coef-
ficients
The results for Poissons equation, while encouraging, are not of particular
interest from the reservoir simulation point of view. In reservoir simulation,
the diffusion coefficient A is typically anisotropic and may have large jumps
from grid cell to grid cell. We will now consider these types of problems.
It has been known for some time that to get acceptable multigrid conver-
gence factors for anisotropic problems, one needs a relaxation scheme more
complicated than simple point (or cell) relaxation. The reason is that simple
point Gauss-Seidel relaxation will smooth the error only between strongly
coupled points. For example, consider the diffusion equation
V AV(x) = /(x), x G 0,
where
Standard finite differences and a standard multigrid algorithm with point
Gauss-Seidel relaxation leads to extremely poor convergence factors when
e points in the x direction than in the y direction. Local mode analysis [3]
shows that point Gauss-Seidel relaxation will only smooth the error effec-
tively in the x direction. Therefore, one cannot properly represent the error
on a grid that is coarser in the y direction. One alternative is to accept
30


this smoothing limitation and only coarsen in the x direction, an approach
that has been used in reservoir simulation in [9]. Another alternative is to
relax as a group all of the points that are strongly coupled, which, in this
case, consists of those points that share the same y coordinate. This x-line
relaxation will smooth the error in the weakly coupled direction (the y direc-
tion). Then one can represent the error well on a grid that is coarser in the
y direction. The distributive relaxation scheme as previously discussed is a
single cell relaxation scheme. Only unknowns associated with a particular
cell are updated together, and we cannot hope to get acceptable multigrid
convergence rates with it for anisotropic problems. The alternative we choose
is line relaxation as discussed in Section 1.2. It may often be the case that
the strong coupling is in the x direction in parts of the domain and in the y
direction in other parts. For these problems, we will need to use alternating
line relaxation. An alternative, that we will not consider in this thesis, would
be to use semi-coarsening in one of the directions and line relaxation in the
other as in [9].
Later we will discuss in more detail problems with the multigrid algorithm
when the diffusion coefficient jumps from cell to cell. However, the algorithm
as presented in the previous section yields quite acceptable convergence rates
when the jumps in the diffusion coefficient occur (if at all) at grid interfaces
on the coarsest grid.
We will now present a single numerical example that illustrates the abil-
31


ity of the multigrid solver to deal with both anisotropic and discontinuous
diffusion coefficients. Consider
/ 1 ax u + dc^/dx = 0
ay~1v + d(f>/dy = 0
du/dx + dv/dy = /
with boundary conditions
u(x) = fi'(x)5 x C <9fiwest U <90east
u(x) 9(x)x C ^^south ^ '^nortlr
The domain of the problem, along with values of ax and ay, is shown in Figure
1.9. Note that this problem has regions of strong coupling in the x direction
and regions of strong coupling in the y direction. The diffusion coefficient
also has lines of discontinuity, although they correspond to grid edges on the
coarsest (2 by 2) grid. To obtain asymptotic estimates, the right hand side
/ and the boundary condition g were set to zero, a random initial guess was
used, and ten or more multigrid cycles were performed. We then calculated
the geometrically averaged convergence factor ( factor by which the residuals
are reduced per cycle). The results are shown below.
Grid Size Convergence Factor
8x8 .04
16 x 16 .10
32 x 32 .18
64 x 64 .20
Alternating line relaxation was used with two pre-relaxations and one
post-relaxation. We see acceptable convergence factors in that they are fairly
small and seem to stabilize as the number of unknowns is increased.
32


Figure 1.9: Problem Domain
33


Comparing the mixed FVE method to standard cell-centered finite differ-
ences, the greater accuracy of the mixed FVE method is much more evident
when the diffusion coefficient is discontinuous. Below we present results for
a problem where A is a scalar, but discontinuous. The domain 0 is [1, l]2
and the diffusion coefficient has the following values:
0.05 for x > 0, y > 0
0.01 for x < 0, y > 0
10.0 for x < 0, y < 0
33.33 for x > 0, y < 0.
Boundary conditions specified the normal component of velocity to be 1
on the left and right edges of the domain and 0 (no flow) on the top and
bottom. Because an analytic solution for the problem was not available, the
mixed FVE solution from a 256x256 grid was used as the exact solution in
calculating errors. The velocity errors, measured in the same way as in the
last section, for the mixed FVE solution and the finite difference solution are
shown below. The finite difference discretization used harmonic averaging of
the diffusion coefficient at interfaces as is standard.
Mixed FVE
Grid Size urnfVe Vmfve ue
2x2 3.26E-1 9.53E-1
4x4 2.25E-2 2.86E-2
8x8 1.04E-2 1.94E-2
16x16 5.01E-3 1.03E-2
32x32 2.46E-3 5.19E-3
64x64 1.18E-3 2.49E-3
128x128 4.89E-4 1.00E-3
256x256 0 0
34


Finite Differences
Grid Size 4d *ld
2x2 8.28E-1 9.53E-1
4x4 4.77E-1 5.03E-1
8x8 2.70E-1 3.07E-1
16x16 1.47E-1 1.69E-1
32x32 7.94E-2 9.02E-2
64x64 4.23E-2 4.77E-2
128x128 2.24E-2 2.51E-2
256x256 1.18E-2 1.35E-2
The superiority of the mixed FVE approach is clear; it computes more
accurate velocities on a 16x16 grid than finite differences does on a 256x256
grid. This is not an artifact of using the 256x256 mixed solution as exact,
since the 128x128 and 256x256 finite difference solutions differ from each
other by an order of magnitude more than the corresponding mixed solutions
differ from each other. As in the comparison of the previous section, the two
methods exhibit comparable accuracy for the pressure (not shown). When
comparing solutions on different size grids, the mixed velocities do not show
second-order convergence as they will in all other tests reported in this thesis.
The reason being that the true solution has a singularity at the origin, and
thus we cannot expect second-order convergence.
1.3.3 An Alternative Multilevel Algorithm
A Two-Level Approach
We pointed out in the last section that the multigrid algorithm presented
in Section 1.2 is robust with respect to jumps in the diffusion coefficient as
35


long as those jumps occur (if at all) at grid edges on the coarsest grid. For
a given diffusion coefficient, this limits the number of levels one can use in
the multigrid algorithm: we can coarsen the grid only to the level of the
jumps in the diffusion coefficient. If the diffusion coefficient has only a few
lines of discontinuity, this is acceptable. However, if the diffusion coefficient
has many lines of discontinuity, the coarsest grid we can use may be too fine;
there may be too many unknowns to solve the coarse grid version of the error
equation well enough, and we must clearly seek an alternative.
The problem that the previous multigrid algorithm has with jumps in the
diffusion coefficient that do not occur at interfaces of the coarsest grid can be
traced to interpolation. The piecewise constant interpolation scheme defined
in (1.15) is suitable so long as the diffusion coefficient is continuous within
the coarse grid cell. To illustrate the problem that arises when this is not
the case, consider the extreme situation shown in Figure 1.10.
In the two fine grid cells on the right, there is no flow and <(> = 0; these
are dead cells. In the two fine grid cells on the left, ^ is non-zero; these are
live cells. If we use piecewise constant interpolation from the coarse grid
cell, which will generally be live, then we can expect to interpolate some
nonzero value of ^ to the two dead cells, which is completely incorrect. The
solution to this problem, as discussed in [1], is to use information about the
diffusion coefficient in interpolating . This is the basis of the algorithm for
black box multigrid [8], a robust multigrid solver for the diffusion equation
36


Figure 1.10: Live and Dead Cells
37


with strongly discontinuous coefficients. It may be possible to incorporate
this operator-based interpolation into the previous multigrid algorithm, but
we choose instead a different more convenient approach.
We will explain the approach as a two-level multigrid algorithm, although
this is somewhat of a misnomer since we will have only one grid. Consider the
fine level problem as the mixed FVE discretization of the first-order system
A-1v + Vcf> =0,
V v = /.
We write the mixed FVE equations in matrix form as
M gradh W vh
divh 0 ) \ h
0
h
(1.22)
Here, M is the mass matrix that comes from the integration of Darcys
equation and gradh and divh are the grid h discrete operators corresponding
to the continuous operators grad, and div. Define the residuals according to,
V
0 M gradh \ / vh
fh ) V divh 0 ) \ h
(1.23)
where the variables with hats denote a current approximate solution to equa-
tion (1.22). Define the errors by
eh = vh Yh^
e4> = h~ h-
The error equation is then,
M gradh
divh 0
"V
0h
v
38
(1.24)


Now, rather than using a coarser grid with the mixed FVE discretization
to approximate the error equation, use the same grid with a standard cell-
centered finite difference approximation. This is the coarse level in our
multigrid algorithm. The coarse grid version of the error equation can then
be written in matrix form as
M gradh
divh 0
(1.25)
The only difference between equations (1.24) and (1.25) is the mass matrix.
In (1.25) the mass matrix M is diagonal and computed from the diffusion
coefficient using harmonic averaging at discontinuities. The point is that, in
equation (1.25), we can eliminate the velocity variables. We have
vh = M 1 (gradh(f)h + r(j) . (1.26)
This allows us to rewrite equation (1.25) as
divhM 1gradhh = divhM Xr\. (1-27)
This is just the type of equation for which black box multigrid was developed.
Our approach is then to use black box multigrid to solve this equation for
(A, use equation (1.26) to get v^, and these approximations of the error to
correct our mixed FVE approximation.
This two-level approach is similar to the work in [15], where black box
multigrid was used as a coarse level for a Lagrangian hydrodynamics appli-
cation. In different vein, this two-level approach is similar to preconditioning
39


the mass matrix M by its diagonal. This approach was suggested, in the
context of mixed finite element methods in [18].
Two-Level Computational Results
We present here the results from a single numerical experiment to test the
robustness of the two-level approach with respect to discontinuities in the
diffusion coefficient. The diffusion coefficients ax and ay were separately
and randomly assigned a value between .01 and 100 for each grid cell. The
diffusion coefficient is thus anisotropic and has jumps of several orders of
magnitude between cells. The two-level algorithm described in the previous
section was used, with two alternating line relaxation sweeps on the mixed
FVE equations before using black box multigrid to solve the coarse problem
and one alternating line relaxation sweep after. One point to consider is: how
well do we need black box multigrid to solve the coarse problem? In [15]
only one cycle of black box multigrid was used for this purpose. For this
difficult test problem ( note: black box multigrid convergence factors were
approximately .6 ), we found that performing more than one cycle of black
box multigrid improved the overall convergence of the two-level method. In
the results reported below, five cycles of black box multigrid are used to
approximate the coarse problem solution, although fewer, say 2 or 3, are
often enough to obtain similar two-level convergence factors. The asymptotic
convergence factors for the two-level method are given below.
40


Grid Size Convergence Factor
16 x 16 .43
32 x 32 .44
64 x 64 .46
We see that, although the factors are much less impressive than they are
for Poisson type applications, especially when work is considered, this two-
level approach does exhibit convergence factors that seem to be bounded with
growing problem size. One point about these relatively poor convergence
factors (for comparison, the convergence factors of this two-level algorithm for
constant coefficient problems is .1), these results show that, for this problem,
there is a significant difference between the mixed FVE discretization and
the standard finite-difference approach. This difference was illustrated in the
numerical tests of Section 1.3.2. When these two approaches are close, the
diagonal mass matrix M is an excellent approximation to the mixed FVE
mass matrix M, and the two-level method converges quickly. Alternatively,
when the two-level method converges relatively slowly, the two approaches
are significantly different, and we have seen that the mixed FVE approach
can be much more accurate.
The point of considering this two-level method was to allow treatment
of problems where the coarsest grid for the multigrid algorithm presented in
Section 1.2 is still too fine and has too many unknowns to solve the mixed
FVE equations by itself. In practice, the multigrid algorithm presented in
Section 1.2 could be used down to the coarsest grid that was aligned with
the discontinuities in the diffusion coefficient, where the two-level approach
41


described here could be applied.
1.3.4 Full Tensor Diffusion Coefficients
In this section we present results of two calculations with full tensor diffu-
sion coefficients on uniform grids. The derivation of the mixed FVE Darcy
equation for this case is discussed in Section 3.1.
In the first problem, the diffusion coefficient A is constant throughout
the domain 0 = [1,1]2 and given by
_ / cos 9 sin 9 \ / 1 0 \ / cos 9 sin 9 \
^ sin# cos 9 J ^ 0 0.01 J ^ sin# cos 9 J
where 9 is the angle between the coordinate axes and the principal directions
of permeability. The boundary conditions and source term / were set so that
the exact solution to the problem (with integral over 0 = 0) is
(x, y) = cos [ttx) cos (27ry).
The mixed FVE discretization was solved on uniform grids of various sizes
and for various values of the angle 9. The errors in the resulting approxima-
tions are shown below (measured in the same way as in previous sections).
o O II
Grid Size (j)mfve urnfVe V'mfve
8x8 0.4804 5.942E-3 3.216E-3
16x16 0.2501 2.140E-3 1.091E-3
32x32 0.1263 5.708E-4 2.868E-4
64x64 0.0633 1.449E-4 7.252E-5
42


9 = 15
Grid Size (j/n fve urnfVe V'mfve
8x8 1.836 1.925E-1 8.872E-2
16x16 0.5551 6.959E-2 3.024E-2
32x32 0.1803 1.949E-2 8.268E-3
64x64 0.0712 5.035E-3 2.119E-3
II CO o o
Grid Size (j/n fve urnfVe Vmfve ue
8x8 4.684 1.084E-1 6.844E-2
16x16 1.250 4.353E-2 2.733E-2
32x32 0.3357 1.391E-2 8.679E-3
64x64 0.1006 3.805E-3 2.368E-3
9 = 45
Grid Size (j/n fve urnfVe Vmfve ue
8x8 6.575 4.624E-1 3.577E-1
16x16 1.912 1.987E-1 1.568E-1
32x32 0.5254 6.324E-2 5.076E-2
64x64 0.1453 1.715E-2 1.385E-2
From the results we see, not surprisingly, that the accuracy of the dis-
cretization degrades as the angle 9 grows. However, the results show second
order convergence in the velocities even when the anisotropy is not aligned
with the grid, so that any practical accuracy could likely be obtained by refin-
ing the grid. When the anisotropy is not aligned with the grid, say, 9 = 45,
the convergence factors for the mixed FVE multigrid algorithm degrade (av-
erage convergence factors of about .5 using alternating line relaxation), just
as the discretization accuracy does. Because the strongly coupled points are
at a 45 angle, convergence factors of .1 could likely be recovered by block
relaxation in this direction. However, this specialized treatment (block relax-
43


ation in the dominant direction) would be difficult to implement in practice,
where the dominant direction would vary throughout the domain.
The second problem is similar to the test problem in [2]. Here the grid is
uniform, and A is a full tensor but no longer constant throughout the domain
0 = [0, l]2. The value of the diffusion coefficient is given by
0 < £ < .5 .5 < £ < 1
The boundary conditions and source term / were set so that the exact solu-
tion is
0 < £ < .5 .5 < x < 1
where the value of C was chosen to make the integral of over 0 vanish.
Below are the mixed FVE discretization errors, measured in the same way
as before, for various grid sizes.
Grid Size 0'infve ufve V'mfve
4x4 9.212E-2 4.555E-3 5.748E-3
8x8 4.622E-2 1.110E-3 1.540E-3
16x16 2.313E-2 2.741E-4 3.931E-4
32x32 1.157E-2 6.825E-5 9.871E-5
64x64 5.784E-2 1.755E-5 2.435E-5
In this problem, we clearly see first-order convergence of the pressure and
second-order convergence of the velocities. For all test problems with analytic
solutions, we have observed second-order convergence of the velocities, and,
44


as will be shown in the next chapter, the mixed FVE discretization can be
extended to general quadrilaterals and still retain this desirable property.
45


2 General Quadrilateral Grids
2.1 The Mixed FVE Discretization
The focus now is extending the discretization technique of Section 1.1 to
irregular grids. We will develop the mixed FVE discretization for a logically
rectangular grid of irregular quadrilaterals, an example of which is shown
in Figure 2.1. Again, we consider the following partial differential equation
defined on some domain 0, now not necessarily rectangular, in 1Z2:
V AV(x) =/(x) x, G 0,
VA(x)(x)-77, = 0 x G dft.
(2.1)
The diffusion coefficient A may be a function of position and is not necessarily
diagonal; however, it is assumed to be a pointwise symmetric, positive definite
matrix.
Again we separate out Darcys law,
A_1v(x) + V(x) = 0, (2.2)
from the conservation of mass equation,
V v(x) = /(x).
(2.3)
46


Figure 2.1: General Quadrilateral Grid
47


We will formulate a mixed FVE discretization for this system of equations
that generalizes the mixed FVE discretization presented in Section 1.1.
Important in developing the discretization for general quadrilaterals is the
mapping relating a general to a reference quadrilateral. Consider the general
quadrilateral P with vertices (V00, J/oo), (io, J/io), (oi, J/oi), and (£11,2/11), as
shown in Figure 2.2.
Let the reference quadrilateral P be the unit square. Then there is a
unique bilinear reference mapping of P onto P given by
x{x, y) = xQQ + (xio xQQ)x + (xoi xQQ)y + (xn xw x01 + x00)xy}
y(x, y) = yoo + {yio yoo)x + (y01 y00)y + (j/n yw y01 + y00)xy.
If P is convex, then this mapping has an inverse Since we restrict ourselves
to convex quadrilaterals, for each (x,y) £ P we have an associated point
(x, y) £ P given by the inverse reference map.
Shown in Figure 2.2 are several vectors that will be useful later in de-
scribing the components of our discretization technique: for each (x,y) £ P,
we define the four vectors
X(x, y) is the image of the unit vector (1, 0) in P,
Y(x,y) is the image of the unit vector (0,1) in P,
Vx(x}y) is a unit vector orthogonal to Y(x}y)}
r/y(x,y) is a unit vector orthogonal to X(x,j/).
To extend the discretization of Section 1.1 to general quadrilaterals, we
need use analogs of the control volumes and finite element spaces used there.
For the latter, we use the lowest order Raviart-Thomas elements on the
quadrilateral elements (see [5],[19]), which are defined as follows. The char-
48


Figure 2.2: General Quadrilateral
49


acteristic functions of the quadrilaterals provide a basis for the finite element
space for . The basis functions for v are best described in terms of associ-
ating degrees of freedom with normal components on edges of quadrilaterals.
A typical basis function for the finite element space for v has support on two
adjacent quadrilaterals, a constant normal component on the edge shared
by the quadrilaterals, and a normal component of zero on other edges. The
magnitude of the basis function is such that the flux on the common edge is
one:
/ v yds = 1.
J edge
These conditions alone do not uniquely determine the basis function; the
following additional condition on the finite element space is needed. Within
any quadrilateral P, we assume that
v 77a,|| Y|| varies linearly with x, constant with y,
v 77^ 11X11 varies linearly with y, constant with x.
A typical basis function is represented in Figure 2.3.
We note that, just as in the uniform grid case, the basis functions have
continuous normal components across grid interfaces. So also will our com-
puted solution. The nodes for the unknowns are shown in Figure 2.1. On
each grid cell, which we denote by P4-j, we let denote the value of the
pressure on the cell, and uhi+ij, vhitJ, and u\-j+i the values of the
flux on their respective edges. This is a slight departure from the method in
Section 1.1 where the unknowns uh and vh were defined to be the velocities:
applying the method of this section to a uniform grid results in the same
50


Figure 2.3: Basis Function
equations as in Section 1.1 except for a scaling of the uh and vh unknowns
by h.
We now need to choose the control volumes. By analogy to the uniform
grid case, the quadrilaterals used to describe the grid are the natural choice
for the control volumes for equation (2.3). This has the benefit that it leads
to a scheme with a local conservation property on these quadrilateral grid
cells. To proceed, we integrate equation (2.3), over each grid cell P4-j:
/ V vdxdy = f.
J P i; j P i; j
Applying the divergence theorem, we get
9P.
yds =
/

51


The left-hand side of this equation is just the sum of the fluxes on edges of
Pj-j, so the discretization of the mass conservation equation is
u
*+i ,j
u
hJ
+ V
i,j +1
V
hJ
If we assume that the function / is (approximately) piecewise constant, then
we can replace the integral on the right hand side by:
fax AREA(P8J).
(With more information about /, we could use a more accurate approxima-
tion of the integral.)
The control volumes for Darcys equation generalize the control volumes
we used for the uniform grid case. Consider two adjacent grid cells, P8_i j and
Pj-j. Then U8J consists of the image of (1/2,1) X (0,1) under the mapping
for Pi-ij and the image of (0,1/2) X (0,1) under the mapping for P8J In
Figure 2.4, volume U8J is the shaded region.
We associate this volume with the vertical edge shared by P8_i j and P4-j,
which the control volume straddles. We also have control volumes associated
with horizontal edges. For adjacent grid cells, P;j_i and P4-j, volume V8J
consists of the image of (0,1) X (1/2,1) under the mapping for P;j_i and the
image of (0,1) X (0,1/2) under the mapping for P8J .
We now need to find a method to generalize the integrations of Section 1.1
that led to the discrete form of Darcys law. It is worthwhile to summarize
what we did there, but now interpreting it from a vector point of view. We
52


Figure 2.4: U Volume
53


took the x-component of Darcys law and integrated it over the U volumes,
which is equivalent to taking the dot product of the Darcy law equation
with the vector (1,0) and integrating the result over the U volumes. In
equation (1.7), we were able to integrate out the term, leaving integrals of
along lines interior to the P volumes. Because on P volumes, in equation (1.8) we could evaluate the line integrals in terms
of the nodal values for suitable vector held, with which to dot Darcys equation, that allows us
to integrate out derivatives of , leaving integrals that can be evaluated in
terms of the nodal values of required because no single constant vector will produce the desired result.
Fortunately, simply scaling the vector held X yields the desired property. We
illustrate this process by obtaining the discrete Darcy equation associated
with the U volume shown in Figure 2.4. We proceed hrst by taking the dot
product of equation (2.2) with c;X(x, y) and integrating over the left half of
Uj-j. Similarly, we take the dot product of equation (2.2) with crX(x, y) and
integrate over the right half of U8J. Here, q and cr are constants chosen to
eliminate terms involving values of ^ on the interface between the two halves
where The details of the necessary integrations are in Section 3.2; we will present
here only the form of the equation which this integration gives rise to. Note
that we perform the same kind of integration for the V volumes, only we
54


take the dot product of Darcys law with a scaling of the vector Y. For the
U volume shown in Figure 2.4, we get a discrete Darcy equation relating the
pressure drop between the two cells to the fluxes on cell edges:
clui-l,j + c2 ui,j + c3ui+l,j
-\-C4Vi-iJ + C5Uj_iJ_|_i + CgViJ + CfVi)j +1
^i-lj) = O5
where \E\ is the length of the edge shared by the two adjacent grid cells.
The values of the coefficients ci,. ., C7 depend on the position of the vertices
defining the two grid cells and on the values of the diffusion coefficient within
the two cells. Note that the structure of connections here is the same as that
obtained for the uniform grid case with tensor permeability in Section 3.1.
The cross terms c4,. ., C7 will generally be nonzero even when the diffusion
coefficient A is diagonal.
2.2 Multilevel Solvers
2.2.1 FVE Based Multigrid
In this section, we discuss the changes in the multigrid algorithm presented
in Section 1.2 needed to deal with general quadrilateral grids. Because the
mixed FVE discretization for general quadrilaterals is simply an extension
of the discretization for uniform grids, only slight modifications to this algo-
rithm must be made. Given a logically rectangular grid of irregular quadri-
laterals, one can generate a finer grid of the same type by dividing each
quadrilateral in the given grid into four quadrilaterals. This division is done
55


by connecting midpoints of opposite sides. In this way, a sequence of finer
grids can be generated from a given coarse grid. Figure 2.5 shows a coarse
grid Q1 and its first refinement Q2. On each of the grids Q1,. ., QN, we apply
the mixed FVE discretization technique of the previous section. Then we can
use a multigrid algorithm to generate a solution to the FVE discretization of
the problem on the finest grid QN.
The basic components of the multigrid algorithm remain virtually un-
changed. Now, because the structure of the connections in the discrete equa-
tions remain unchanged, any of the relaxation techniques discussed in Sec-
tion 1.2 can be used. As was the case for uniform grids, the interpolation
operators used in the multigrid algorithm are defined on the basis of the re-
lationship between finite element spaces for the unknowns on different grids.
The finite element spaces are still nested, as we show in Section 3.3. The
relationship between the finite element spaces for is the same as in the
uniform grid case, so the interpolation operator for remains unchanged.
However, the relationship between the coarse and fine finite element spaces
for both u and v is changed somewhat, because these variables now represent
fluxes rather than velocities as before. Since the edge of a fine grid cell is half
the length of the edge of the corresponding coarse grid cell, the interpolation
operators for u and v of Section 1.2 must be scaled by 1/2.
Again, the restriction operators used in the multigrid algorithm are de-
fined on the basis of the relationship between control volumes on different
56


Figure 2.5: O1 and Its First Refinement Q2
57


grids. The control volumes for the continuity equation are related in the
same way as in the uniform grid case (one coarse grid P volume consists of
four fine grid P volumes), and so the restriction operator for this equation
remains unchanged. The relationship between coarse and fine control vol-
umes for Darcys equation can be different than it was in the uniform grid
case. A coarse grid U volume contains all of two of the fine grid U volumes
and fractions (no longer necessarily half) of four others. Figure 2.6 depicts
the relationship between coarse and fine control volumes. Strict adherence
to the FVE principles requires the calculation of the fraction of a fine grid
control volume included in the coarse grid control volume; however, we have
chosen to use the same restriction operator as in the uniform grid case. It is
certainly possible to calculate the fractions and use them in the definition of
the restriction operator, but, for the present, we have used this more easily
implemented definition (approximating the fractions by |). From our expe-
rience with numerical testing of the algorithm, this approximation does not
appear to adversely affect the convergence of the multigrid algorithm. In
fact, because the convergence factors are generally quite small, it is not clear
that there is much to be gained by stricter adherence to the FVE principles
in defining the restriction operator.
Note that in the multigrid algorithm presented here, the process of devel-
oping its components begins with the coarsest grid. Input to the algorithm
is a coarse grid that models the large-scale irregularities of the reservoir ge-
58


Figure 2.6: Coarse and Fine Darcy Volumes
59


ometry, and a specified level of refinement necessary to resolve the flow. The
algorithm then generates a solution to the mixed FVE discretization on the
finest grid. If a direct method is used to solve the equations on the coarsest
grid, then this places a practical limitation on the number of quadrilaterals
used to model the reservoir geometry.
2.2.2 A Two-Level Approach
There are two practical problems with the multigrid algorithm discussed in
the previous section. The first is that, as in the uniform grid case, the jumps
in the diffusion coefficient must occur (if at all) at grid edges on the coarse
grid. The second is that the grid irregularity is constructed based on a given
coarse grid, which is then refined by bilinear coordinates to generate finer
grids: one cannot apply the multigrid algorithm to the equations on the
coarsest grid, so they must be solved some other way. In practice, both of
these problems represent limitations on the coarsest grid: it must be fine
enough to capture the reservoir geology and the jumps in the diffusion co-
efficient. This limitation may result in the set of discrete equations on the
coarsest grid being too large to solve directly.
The two-level approach described in Section 1.3.3 provides one possible
remedy to these two problems. The only component of the two-level algo-
rithm that needs additional discussion for the general quadrilateral case is
M, which is the mass matrix that multiplies the discrete fluxes in Darcys
equation. For uniform grids, M is defined by equation (1.9), and since its
60


diagonal, we invert it to eliminate the flux variables. In essence, we can use
a standard cell-centered finite difference discretization as the coarse level
for the fine level mixed FVE discretization. We would like do something
similar for the case of a general quadrilateral grids. However, one of our mo-
tivations for looking at the mixed FVE discretization is that it can be applied
in a clear and direct way to general quadrilateral grids, where standard cell-
centered finite differences cannot. Unfortunately, an effective cell-centered
finite difference scheme for discretization on general quadrilateral grids does
not currently exist. However, our demands on this scheme are not very strict
because the objective of the discretization on the coarse level is simply to
correct the solution from the fine level mixed FVE discretization, which
will serve as our final approximation. In other words, we would like to use
the coarse level discretization only to accelerate the relaxation process on
the fine level. For this purpose, our experience shows that equation (1.9)
serves as an effective way to define M even in the general quadrilateral case.
There are perhaps more sophisticated ways of defining M that we may ex-
plore later, but we have found that this simple definition works well for many
applications. It is clear that, for very distorted grids, M defined this way
will provide a poor approximation to Af; however, we will see in the next
section that, for moderately distorted grids, the two-level method works as
well as in the uniform grid case.
61


2.3 Computational Results
2.3.1 Multigrid Performance
This section contains results of numerical experiments designed to test the
robustness of the multigrid algorithm of Section 2.2.1. For all problems
considered here we used two pre-relaxations and one post-relaxation in the
multigrid algorithm. For the first test problem, we used the mixed FVE based
multigrid algorithm with a coarsest grid of eight cells in each direction. We
began with a uniform 8x8 grid on 0 = [1,1] X [1,1] and then moved each
interior vertex in both the x and y direction separately by a random number
between .2h and .2h, where the original mesh size is h = .25. The resulting
grid is shown in Figure 2.7.
We then solved the mixed FVE discretization of Poissons equation, where
the diffusion coefficient is the identity, using the multigrid algorithm of Sec-
tion 2.2 with varying levels of refinement. Ten or more cycles were performed
on the homogeneous problem and the average convergence factor was calcu-
lated. The results are shown below.
Grid Size Convergence Factor
16 x 16 .12
32 x 32 .17
64 x 64 .22
Here we used distributive Gauss-Seidel relaxation and we see convergence
factors that are quite acceptable for this distorted grid in that they are rel-
atively small, do not grow significantly with problem size, and are not much
62


Figure 2.7: Grid with up to 20% Distortion
63


different from the convergence factors for the uniform grid case. The algo-
rithm allows us to solve this problem with about as much work as is needed
to solve the uniform grid problem.
For the second numerical experiment, we solved the same problem, only
now the 8x8 grid was distorted by moving each interior vertex in both the x
and y direction by a random number between .5h and .5h. The resulting
grid is shown in Figure 2.7.
This grid is likely to be much more distorted than one would encounter in
an application, and, in fact, some of the quadrilaterals are no longer convex.
The results for the mixed FVE based multigrid algorithm are shown below.
Grid Size Convergence Factor
16 x 16 .03
32 x 32 .05
64 x 64 .08
Because of the significant distortion of the coarsest grid, we used alternat-
ing line relaxation for this test example: analogous to the anisotropic uniform
grid case, the structure of the resulting discrete mixed FVE equations con-
sists of cells that are strongly coupled to some neighbors and weakly coupled
to others. We need alternating line relaxation to get the error smoothing
that is needed for acceptable multigrid convergence factors. With this more
expensive algorithm, compared to the cell by cell distributive Gauss-Seidel
relaxation, we obtain quite good multigrid convergence factors for this very
distorted grid.
The third numerical test combined moderately distorted grids with a
64


Figure 2.8: Grid with up to 50% Distortion
65


diagonal but anisotropic and discontinuous diffusion coefficient. The coarsest
grid here is the same as that used in the first experiment, an 8x8 grid with up
to 20% distortion. On each cell in this coarsest grid the diffusion coefficients
ax and ay were separately set to random values between .01 and 100. The
results for the mixed FVE based multigrid algorithm are shown below.
Grid Size Convergence Factor
16 x 16 .04
32 x 32 .06
64 x 64 .08
Again we used alternating line relaxation, now because of the anisotropic
diffusion coefficient. Again we see quite good multigrid convergence factors.
This numerical test is close to the types of problems that might actually
occur in an application and the performance of the multigrid algorithm is
well within an acceptable range.
The two restrictions on the coarsest grid for our mixed FVE based multi-
grid method that were discussed in Section 2.2.2, namely, capturing the ir-
regularity and discontinuities of the diffusion coefficient, have not been vio-
lated in the numerical experiments thus far. However, even without violating
these limitations, one can pose problems where the mixed FVE based multi-
grid method fails to perform as well as in the above tests. For example,
consider the coarsest grid used in second numerical experiment: 8x8 with
up to 50% distortion. On each cell in this grid, the diffusion coefficients ax
and ay were separately set to random values between 1/a and a. The results
for the mixed FVE based multigrid algorithm using four grids (finest grid is
66


64x64) are shown below.
a Convergence Factor
1 .08
10 .13
100 .15
1000 .20
These last results may be of little practical interest as the grid is much
more distorted than one would likely use in an application, but they do
demonstrate the robustness of the multigrid algorithm. For all numeri-
cal tests encountered thus far, the mixed FVE based multigrid algorithm
with alternating line relaxation has yielded acceptable (rarely worse than .3)
convergence factors if the lines of discontinuity in the diffusion coefficient
are captured on the coarsest grid. With full tensor coefficients, when the
anisotropy is not aligned with the coordinates, the convergence factors can
be significantly larger (at times, as large as .6), and we have no rigorous
theoretical results, so we can make no general claims. However, it appears
that the method has potential as a fast solver in many important practical
applications.
2.3.2 Two-Level Performance
This section contains results of numerical tests designed to test the robustness
of the two-level approach of Section 2.2.2. This is an approach that one could
use alone or as the coarsest grid solver in the mixed FVE based multigrid. For
the results of this section, we use alternating line relaxation, with two sweeps
67


performed before solving the coarse problem and one sweep performed after.
Five cycles of black box multigrid were used to solve the coarse problem,
although often two or three are sufficient. In the first numerical experiment,
the mixed FVE discretization of Poissons equation on 0 = [1,1] X [1,1]
was solved on the moderately distorted grid obtained from a uniform grid
by the up to 20% distortion described in the previous section. Note that
the distortion here takes place on the finest grid, so the mixed FVE based
multigrid algorithm cannot be used. Average convergence factors for the
two-level approach on different size grids are shown below.
Grid Size Convergence Factor
16 x 16 .07
32 x 32 .08
64 x 64 .08
The convergence factors are remarkably good, given the quality of the
approximation used on the coarse level. As discussed in Section 2.2.2, we
basically assume the grid is uniform in constructing the mass matrix M for
the discrete Darcy equations on the coarse level. This appears to work well
for the moderately distorted grids like those in this numerical experiment
and, quite probably, the grids one would use in practical applications. When
the grid is very distorted, like the up to 50% distorted grid of the previous
section, the two-level algorithm can fail to converge and may even diverge.
The reason is that the very poor approximation of the mass matrix results
in a correction from the coarse level that has little, if anything, to do with
the fine level error. It is possible that this could be remedied by a more
68


sophisticated choice for M, but this has yet to be investigated.
In the next numerical experiment, we use the same grids as in the previous
test and solve the mixed FVE discretization to the diffusion equation with
diagonal diffusion coefficient, where on each cell in this grid the diffusion
coefficients ax and ay were separately set to random values between .01 and
100. The results, in terms of average convergence factors, are shown below.
Grid Size Convergence Factor
16 x 16 .43
32 x 32 .40
64 x 64 .38
While the convergence factors are not as small as for Poisson-like prob-
lems, they are likely acceptable in most practical applications, especially if
one is using the two-level approach only on the coarsest grid of the mixed
FVE based multigrid algorithm. There the amount of work on finer grids
in the mixed FVE based multigrid algorithm will typically be much larger
than the work of the two-level algorithm on the coarsest grid, even if sev-
eral cycles of the two-level algorithm are required. We should note that the
convergence factors for this problem are roughly equal to those obtained for
the similar uniform grid problem: our belief is that the poor convergence
factors are again due to the difference between the mixed FVE methods and
the finite difference methods treatment of discontinuities in A, not related
to the moderate irregularity of the grid.
69


2.3.3 Accuracy of the Mixed FVE Discretization
In this section, we present results of numerical experiments to investigate the
accuracy of the mixed FVE discretization on general quadrilateral grids.
First, we will present the result of a simple numerical test designed to
illustrate the effect that the shape of the grid cells has on discretization errors
and multigrid convergence rates. We consider Poissons equation posed on
Q = [-l,l]2:
v + V = 0,
v V = /,
with boundary condition v rj = g. For this test, / and g were determined
so that the exact solution is,
(x, y) = cos(7rx) cos(27ry) + constant.
The solution from the mixed FVE discretization was normalized to corre-
spond the exact solution with constant= 0. The coarsest grid we use in our
tests is shown in Figure 2.9.
To investigate the effect of the angle f3 on the accuracy of the mixed FVE
discretization and on multigrid convergence factors, we use varying /3. We
measure the error in the same way as in Section 1.3, using the C2 norm of the
error and C2 norm of the u and v errors. In the tables below, we compile
our observations on these errors and the average convergence factors for the
mixed FVE based multigrid algorithm. For f3 = 80 and 70, distributive
Gauss-Seidel relaxation proved effective in terms of convergence factors. For
70


Figure 2.9: Coarsest Grid
71


f3 = 60 and 50, the stronger alternating line relaxation scheme was needed
to get acceptable convergence factors. This is again due to the distortion of
the grids that yield couplings of different strengths to neighboring cells.
f3 = 80 distributive Gauss-Seidel
Fine Grid ue Ve Convergence Factor
8x8 4.841E-1 1.575E-1 7.595E-2 .054
16x16 2.514E-1 5.090E-2 2.538E-2 .073
32x32 1.267E-1 1.347E-2 6.752E-3 .083
64x64 6.354E-2 3.414E-3 1.713E-3 .086
f3 = 70 distributive Gauss-Seidel
Fine Grid ue Ve Convergence Factor
8x8 4.881E-1 2.327E-1 1.106E-1 .067
16x16 2.541E-1 6.541E-2 3.330E-2 .090
32x32 1.285E-1 1.689E-2 8.809E-3 .107
64x64 6.516E-2 4.256E-3 2.235E-3 .143
f3 = 60 alternating line relaxation
Fine Grid Ue Ve Convergence Factor
8x8 5.009E-1 3.506E-1 2.140E-1 .019
16x16 2.619E-1 9.111E-2 5.723E-2 .031
32x32 1.335E-1 2.305E-2 1.483-2 .035
64x64 6.967E-2 5.779E-3 3.745E-3 .054
f3 = 50 alternating line relaxation
Fine Grid Ue Ve Convergence Factor
8x8 5.318E-1 6.060E-1 4.169E-1 .022
16x16 2.760E-1 1.737E-1 1.199E-1 .038
32x32 1.406E-1 4.493E-2 3.189E-2 .057
64x64 7.382E-2 1.135E-2 8.142E-3 .097
These results, show good convergence in the accuracy of the discretization
as the grids are refined. Each refinement reduces the norm of the error in
by a factor of 2 and the error in the fluxes by a factor of 4. Although
72


we do not have rigorous proofs of the order of convergence of the mixed
FVE discretization, these results, and all other results from experiments with
analytic solutions, appear to indicate that, at least in the norms we used,
the convergence is first order in the pressure and second order in the fluxes.
We also see that the accuracy of the discretization decreases as f3 decreases,
which is not surprising since most finite element approximations degrade
as the grid becomes less uniform. However, one can still likely achieve a
solution of any practical accuracy by refining the grid yet further. Note also
that the multigrid convergence factors degrade slightly as the grid becomes
less uniform, but the factors are still certainly acceptable.
The next problem is similar to the previous one in that the same type
of grids are used, and the boundary condition and right hand side were
determined so that it has the same exact solution; however, this problem has
full tensor permeability. The diffusion coefficient A is given by
_ / cos 9 sin 9 \ / 1 0 \ / cos 9 sin 9 \
^ sin# cos 9 J ^ 0 0.01 J ^ sin# cos 9 J
where 0, the angle between the coordinate axes and the principal directions
of permeability, is 45. The errors in the mixed FVE approximations for
various size grids and varying angle of distortion f3 are given below:
'CO II 00 o o
Fine Grid Size ue ve
16x16 2.261 2.722E-1 1.998E-1
32x32 0.6356 9.497E-2 7.194E-2
64x64 0.1743 2.714E-2 2.078E-2
128x128 0.0520 7.042E-3 5.442E-3
73


'CO II cs o o
Fine Grid Size ue Ve
16x16 3.092 4.527E-1 3.021E-1
32x32 0.9468 1.963E-1 1.398E-1
64x64 0.2698 6.375E-2 4.568E-2
128x128 0.0788 1.747E-2 1.267E-2
As in the uniform grid case, the accuracy degrades due to the anisotropy
not being aligned with the coordinates, but the results still indicate second-
order convergence in the velocities. Note that for f3 = 60 and 6 = 45, a
case with severe distortion and non-grid aligned anisotropy, the asymptotic
regime is not reached until the grid is quite fine.
The final test problem involves a non-uniform grid and full tensor, discon-
tinuous permeability. The domain for the problem, which consists of three
regions, is depicted in Figure 2.10.
The coefficient A is defined in each region as follows:
In region I In region II In region III
/ 2 1
l -
4 4
U 4
/2 |
1 i
2 2
VI 1
The boundary conditions and right hand side were defined so that the exact
solution is given by,
In region II In region III
In region I
6 = x2 + C
= T-y+c
= b2 + c
Again the value of C is chosen to correspond to the solution whose integral
over the domain vanishes. Note that, and the normal component of AV
are continuous across the interfaces between regions. The coarsest grid (a
74


Figure 2.10: Problem Domain
75


2x2 grid), was formed by splitting, along the x-axis, region I into two cells,
and all finer grids were generated by bilinear refinement of this coarsest grid.
The errors in the mixed FVE approximation, calculated in the same way as
previously, are presented in the table below.
Grid Size ue Ve
4x4 0.2524 7.065E-3 8.269E-3
8x8 0.1302 1.747E-3 2.197E-3
16x16 0.0658 4.567E-4 5.859E-4
32x32 0.0333 1.213E-4 1.582E-4
64x64 0.0174 3.234E-5 4.294E-5
Here we clearly see first-order convergence in the pressure, and again we
see second-order convergence in the velocities, which has been the case for
all experiments, ranging from Poissons equation on a uniform grid to this
significantly more complicated case, with analytic solutions.
76


3 Technical Derivations
This chapter contains all the derivations needed to assemble the mixed
FVE equations. Many of the results from this section have been briefly
presented, or at least alluded to, in previous sections.
3.1 Uniform Grids and Tensor Permeability
In Section 1.1, we presented the FVE discretization assuming that the coef-
ficient A in equation (1.1) is a diagonal tensor. Here we consider the case
where A is a full tensor, but we assume that it is a pointwise symmetric,
positive definite matrix. In this section, we assume the grid is uniform (see
the following section for the general quadrilateral case). We write
A
a b
b c
where a, 6, and c may now be functions of position. Darcys law, equa-
tion (1.2), takes the form
v = -AV(>,
77


or
A 1v + V = 0.
Now, adding the conservation of mass equation and noting that v = (it, n)4,
we have the following first-order system:
(detA)-1 (cit bv) -\-d(f)/dx =0 (u equation)
< (detA)-1 ( bu + av) -\-d du/dx -\-dv/dy = f (p equation).
The main components of the mixed FVE discretization process remain the
same: the choice of volumes and finite element spaces are the same as pre-
sented in Section 1.1. The change in the form of the u equation, comparing
equation (3.1) to equation (1.4), is reflected when we integrate the u equa-
tion over the U volumes. In this section, we will assume that A is constant
over the integration domain, otherwise the technique of Section 3.3 can be
applied. As an example, let U8J be the volume in hi that is centered at the
interior uh node (i,j) as shown in Figure 3.1.
We then have
/ ((detA) 1(cu bv) + d(f>/dx^ dxdy = 0. (3-2)
Integration of the terms involving u and d(p/dx is the same as in Section 1.1.
The additional term involving v is handled by splitting the integral into two
parts, one for each half of U8J. Then, because our finite element approxima-
tion vh is linear in y and constant in x on each half, by direct calculation of
78


Figure 3.1: U Volume
79


the integral over the left half, we have
T

Combining this result with the analogous one for the other half and adding
the results from Section 1.1 for the other terms in equation (3.2), we get the
discrete equation
V(c(det A)"1) (tw +
-T (fc(det A)-1) ("tw+i + ?-W + <-+. + C)
+ h (jt, ~ ti-u) =
We get a similar equation when we integrate over V volumes. This change
in the form of the discrete Darcy law is the only change in the discretization
process from Section 1.1.
3.2 Derivation of FVE Stencils for General
Quadrilaterals
3.2.1 Scalar Diffusion Coefficient
Here we develop the details of the integrals encountered in deriving the dis-
crete Darcy equation for general quadrilateral grids. We first consider the
case where the diffusion coefficient is a scalar. The case where the diffusion
coefficient is a tensor is covered in Section 3.2.2. Shown in Figure 3.2 by
dotted lines is a typical U volume, which we will use to discretize Darcys
equation.
80


Figure 3.2: U Volume
81


P; and Pr will denote the two P volumes who share the edge with which
the U volume is associated. The continuous Darcy equation is
a 1v(z,i/) + V0(x,y) = O.
(3.3)
As mentioned in Section 2.1, we take the dot product of equation (3.3) with
QX(x,y) and integrate over the left half of the U volume. Similarly, we
take the dot product of equation (3.3) with crX(x,y) and integrate over the
right half of the U volume. We then add the two integrals to get the discrete
Darcy equation corresponding to this U volume. Letting U; and Ur denote
the two halves of the U volume we have the following:
fu, (a Ma V) + y)) ciX(x} y)dxdy+
fur (a_1v(x, y) + V(x, y)) crX(x, y)dxdy =
(3.4)
Q/\ja 1v(x} y) X(z, y)dxdy + q/U; V^(z, y) X(x, y)dxdy+
Crhj^vix, y) X(x, y)dxdy + crfUrV(f>(x, y) X(x, y)dxdy.
Thus, we have two types of integrals to evaluate, ones involving v and ones
involving Xcf>.
First consider the integral involving v. For the moment, we ignore the
scaling parameters q and cr. We evaluate this integral for the finite element
basis functions for v, which are the formulas needed to assemble the stencil
for equation (3.4). Consider the evaluation of
where w is the basis function with flux equal to 1 on the edge straddled by
the U volume and 0 on all other grid edges. It is important to recall the
82


properties used to define our finite element space: within each grid cell,
w(x, y) yx(x, y) || Y(x, y) || varies linearly with x, constant with j/,
w(x, y) r]y(x, y) ||X(x, y) || varies linearly with y, constant with x.
(3.5)
The first condition implies, that within the volume P/,
w(x,y) rjx(x,y)\\Y(x,y)\\ = x,
where x refers to the pre-image of x in the reference space. Similarly, y will
denote the pre-image of y in the reference space. The second condition, along
with the fact that fluxes are zero on the top and bottom edges of P/, implies
that, within the volume P/, we have
w(x,y) T]y(x} 2/)||X(m, y)\\ = 0.
Thus, w is perpendicular to rjy and, hence, parallel to X. We thus have,
w(x,y) X(x,y) = ||w(x, y)|| ||X(x, y)||
and
w(x,y) rjx{x,y) = ||w(x,y)|| sin(0(x,y)),
where 9 is the angle between X and Y as shown in Figure 3.2. Solving this
last equation for the norm, we get
w(z, y)|| = (w(x, y) yx(x, y j) / sin(0(a;, y j).
(3.6)
83


Putting these results together, we have,
fu ,a_1w(a;, y) X(x, y)dxdy
= /u,a_1M>2/)lll|X(z,2/)|| (3.7)
= /uffi 1 ((w(a;, y) dx{x, y)) / sin(0(x, y))) ||X(x, y)\\dxdy
= /o/i_1||X(x, y)\\2dxdy.
The last step is the result of the change of variables from the actual space to
the reference space. The Jacobian for the map is
J = 11X111|Y|| sin 0.
This final integral can now be evaluated analytically in the reference space.
We had defined the vector held X(x,j/) by taking images of the unit vector
(1,0) under the bilinear mapping
x(x, y) = xQQ + (xio xQQ)x + (xoi xQQ)y + (xn xw x01 + x00)xy}
y(x, y) = Voo + {yio yoo)x + (y01 y00)y + (yn yw y01 + y00)xy.
We can express the vector held X as follows:
X = (dx/dx7dy/dx)
= ((x10 Xoo) + (xn x10 Xoi + oo)y, (3.8)
(yio yoo) + (yu yio yoi + 2/00)2/)-
Thus, the last integrand in equation (3.7) is a polynomial in x and y whose
coefficients depend on the coordinates of the vertices of P; and this integral
can easily be evaluated analytically.
Next consider the evaluation of
/ a~1w(x, y) X(x, y)dxdy,
84


where w is the basis function with flux equal to 1 on the top edge of the left
grid cell and 0 on all other grid edges. The second condition in (3.5) implies
that, within the volume P/, we have
w(x,y) T]y(x} 2/)||X(m, y)\\ = y.
The first condition in equation (3.5), along with the fact that the fluxes are
zero on the left and right edges of P/, implies that, within the volume P/, we
have
w(x,y) f]x(x,y)\\Y(x,y)\\ = 0.
Thus, w is perpendicular to rjx and, hence, parallel to Y. We thus have,
w(ab2/)/llw(ab2/)ll = Y{x,y)/\\Y{x,y)\\
and
w(x,y) rjy(x,y) = ||w(x,y)|| sin(0(x,y)).
Solving this last equation for the norm, we get
\\w(X, j/)|| = (w(x, y) f]y(x, y)) / sin(0(a;, y)). (3.9)
Putting these results together, we have,
/U;alw(x, y) X(x, y)dxdy
= llw(^7 2/)II (w(x,j/)/||w(x,j/)||) X(x,y)dxdy
= IU;a_1 ((w(5 V) Vy{x, y)) / sin(6'(x, y))) (3.10)
(Y(x,y)/\\Y(x,y)\\)-X(x,y)dxdy
= /o/l_1yX Ydxdy
85


Again, the last step is a result of the change of variables. This integral can
now be evaluated analytically in the reference space. The vector held Y can
be written in terms of reference space coordinates as follows:
Y = {dxjdy, dy/dy)
= ((xoi Xoo) + (xn x10 Xoi + x00)x, (3-11)
(s/oi yoo) + (yn 2/10 2/01 + yoo)x).
Thus, the last integrand in equation (3.10) is a polynomial in x and y whose
coefficients depend on the coordinates of the vertices of P/, and this integral
can easily be evaluated analytically. Integrals for other basis functions are
analogous to one of the two that we presented because they can be evaluated
in the same way with a simple change of variables. The basis functions that
will have nonzero contribution to the integral
j a 1w(x, y) X(x, y)dxdy
will be those that are associated with one of the edges of either P; or Pr.
There are seven such edges, so the discrete version of Darcys law on the
volume U will involve seven nodal values for fluxes.
Now consider the terms involving V, for example,
/ V(:r, y) X(x, y)dxdy.
J U,
86


We again express X in of terms reference space variables to get
fu, y) X(a;, y)dxdy
= ful (d(f>/dx, d

= + l^%dxdV (3.12)
= f xjydcf)/dxdxdy
= /o/iff J dxdy.
Here, J is again the Jacobian of the mapping, which can be expressed as
j dx dy dx dy
dx dy dy dx
Note that, expressed in terms of reference space variables, we have
J(x}y) = a + bx + cy7
where a,b, and c are constants on the quadrilateral P; and depend on the
location of its vertices. With this and the last result in equation (3.12), we
have
fu, V(£, y) X(x, y)dxdy
= dxdy
= foil HdS:dS:
= flJ(i,y)(i,y)dy SlJ{\,y){\,y)dy bflf\dxdy.
We have now accomplished the goal of integrating out the derivatives of .
The second integral is a line integral interior to the volume P; and can be
87


written in terms of nodal values of ; also, the third integral is a volume
integral within the volume P; and can be written in terms of nodal values
of . Because the approximation of is constant on Pi (constant= i,j), for
the third integral we have
However, the first integral is a line integral along the edge shared by P; and
Pr, and because we approximate with a piecewise constant function, we
seek to remove the first integral. In terms of the Jacobian of the map for the
right grid cell in Figure 3.2, similar integrations over Ur produces
JurV(K, y) X(a;, y)dxdy =
foJ(b $)(%, y)dy foJ(0, y)(0, y)dy bflf$ (f)dxdy.
In this expression, the first and third integral can be expressed in terms of
nodal values of ; the second integral involves integration along the edge
shared by P; and Pr. It is at this point that we see the utility of the scaling
parameters q and cr. We choose them so that when we add the results from
the two halves, the terms involving on the interface between the two grid
cells cancel. In fact, choosing q = \E\/J(1, |) and cr = \E\/J(0, |) gives us
the result
ci V(f)(x,y)-X(x,y)dxdy+cr 'V(x,y)-X(x,y)dxdy = \E\{cj)l+1^-cj)hl).
J TJ i J TJ r
Here, \E\ is the length of the edge shared by P; and Pr.
88


3.2.2 Tensor Diffusion Coefficient
For a full tensor coefficient, much of the results from the previous section
are unchanged, including, for example, all of the results for the part of the
integral involving . We must only reconsider the evaluation of
/ A-1 w(A, y) X(x, y)dxdy
J u
when w is one of the two basis functions discussed in the previous section.
Again we split the volume into a left and right half and integrate over each
half separately. We will use the same notation as in the previous section, so
that Figure 3.2 still represents the U volume that we integrate over.
We first consider the case where w is the basis function that has flux
equal to 1 on the edge straddled by the U volume and 0 on all other grid
edges. We proceed as follows:
fu, Alw(x, y) X(x, y)dxdy
= fu,w(x>y) ATX(x, y)dxdy
= /uJlw(A2/)ll (w(m, 2/)/||w(m, y) ||) A ~TX(x,y)dxdy
= fu, ((W(A v) Vx(x, y)) / sin(0(a;, y))) (3.13)
(X(a;, y)/\\X(x, y) ||) ATX(x, y)dxdy
= foflxX(x,y) A~TX(x,y)dxdy
= foflxA^Xix, y) X(x, y)dxdy.
The hrst step follows from the adjoint property; the third step uses equa-
89


tion (3.6) and the fact that w and X are parallel; the fourth step is the
result of the change of variables to the reference space; and the final step
again follows from the adjoint property. We can then evaluate the final inte-
gral in equation (3.13) in the reference space using equation (3.8) to express
X in terms of the reference variables. The integrand, as in the previous sec-
tion, is just a polynomial in x and y, so this integral can easily be evaluated
analytically.
We now consider the case where w is the basis function that has flux
equal to 1 on the top edge of the left grid cell and 0 on all other grid edges.
We proceed as follows:
fu, Alw(x, y) X(x, y)dxdy
= fu,w(xiy) ATX(x, y)dxdy
= /uJIw(a?/)|| (w(ab2/)/llw(ab2/)ll) A~TX(x,y)dxdy
= fu, ((w(> V) ' l/)) / sin(0(a;, y))) (3.14)
(Y(x, j/)/||Y(x, j/)||) A-TX(x,y)dxdy
= Jo/lj/Y(x,y)- A~TX(x, y)dxdy
= JlJ'ixA-1Y(x, y) X(x, y)dxdy.
The steps are the same as in the previous case except for step three, which
uses equation (3.9) and the fact that w and Y are parallel. Again, the
final integrand is a polynomial in the reference variables and can easily be
evaluated analytically.
90


Note that for given A, one could consider choosing the quadrilaterals such
that A-1 Y X, making this integral vanish. This is unlikely to be useful in
complex problems, but if A is constant for example, it could be of interest.
Integrals for other basis functions can be evaluated in a manner similar
to one of the two we presented here by a simple change of variables.
3.3 Using Fine Grid Integrals to Generate
Coarse Grid Integrals
In this section, we will show how to generate coarse grid integrals from fine
grid integrals, which is important when the permeability coefficient A is
piecewise constant on the fine grid cells but not on the coarse grid cells.
We begin by examining the relationship between fine and coarse grid
finite element spaces. As mentioned in Section 2.2, the finite element spaces
are nested; we will now make this statement more precise. In Figure 3.3, we
depict a coarse grid volume Pc and the four fine grid volumes Py i = 1, 2, 3, 4,
that are generated by bilinear refinement.
Letting Wc denote the coarse grid velocity basis function with nonzero
flux on the north edge of Pc and similarly letting w4- denote the fine grid
velocity basis function with nonzero flux on the north edge on Py we will
prove the following result: within Pc, we have
Wc = ^ (wj + w2) + ^ (w3 + w4). (3.15)
We prove this result by verifying that this linear combination of fine grid
91


Figure 3.3: Coarse and Fine Volumes
92


basis functions has all of the properties that the coarse grid basis function is
to have, namely:
1) The fluxes are zero on the east, south, and west edges of Pc.
2) The flux is one on the north edge of Pc-
3) The basis function is parallel to the vector held Y£/),
i.e., the basis function is aligned with the bilinear coordinates of Pc
that correspond to vertical lines in the reference space.
4) The quantity W t^HXH varies linearly with y, and is constant
with x, where all quantities are coarse grid quantities.
(3.16)
The first property follows from the fact that each fine grid basis function,
wi = 1, 2, 3, 4, has zero flux on the east, south, and west edges of Pc- Thus
the linear combination also has this property.
We now verify the second property in (3.16):
fpn (§ (wi + w2) + |(w3 + w4)) yds
= 2 Jp" (wi y) ds T ^fpn (w2 y) ds
= |(1) + |(1) = 1-
The first step follows from the fact that both w3 and w4 have zero flux on
the north edge of Pc, and the second step follows from the fact that both
w4 and w2 have flux = 1 on the north edges of P4 and P2, respectively.
The third property in (3.16) follows from the agreement of coarse and
fine grid bilinear coordinates.
We are left with the fourth property in (3.16), to verify this, we will need
the following result:
If (x, y) G Pi, then Xc(z, y) = 2X8(x, y). (3-17)
93


Full Text

PAGE 1

#2)314)5%6))6))4)+%3 1%0#//0#)/%*#%4%1163 5)6%/)74*%0%7)3# 8)8 749.:, 749.:. # 1;7 "/3" 3* # 9..< %01%"22"13"40"%01"$%"0224

PAGE 2

3* 8)8 6 71/! # === 5> 8?0@' 10 8 %01%"22"13"40"%01"$%"0224

PAGE 3

88)&*3#' #15)#/ 15* = *71/! # 70#/ #! ""6( A ""/! B"0" "$" 7"" """ $ "*CA %01%"22"13"40"%01"$%"0224

PAGE 4

A """ C 7 71/! %01%"22"13"40"%01"$%"0224

PAGE 5

2 )8%24444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444) 9915)3$9 9D#15) #9E 9D/0D< 9D9/13**D< 9DD3#3/FFFFEG 9DE##""#E< 9DH13/HD /890%8%244444444444444444444444444444444444444444444444444444444444444444'* D915)3$H, DD"7"<< DD915) << DDD#(6"#,G DE/0,D DE9*,D %01%"22"13"40"%01"$%"0224

PAGE 6

DED(6"*,I DEE#15)3$IG 1"1:6244444444444444444444444444444444444444444444444444444444444444444444444444444-E9;*II ED3"15)7;@:G ED973/:G EDD3/:. EE1;//;FFF.9 '0%00;<444444444444444444444444444444444444444444444444444444444444444445$3"4444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444)!! %01%"22"13"40"%01"$%"0224

PAGE 7

#/J4%?6)3;))47 "( 1!*7"/!" !B! */!" !*0 !?" "!!*8 0*7"7?"" ((!B! !!3*0)0 746"A 1! (31!( (! " %01%"22"13"40"%01"$%"0224

PAGE 8

)8%2 )4)"=%>0:2? 17 "A $K" +" $" ? ( (5(#&'5>&'LM&' 2N &'"O&'P .&' &' "A ( 9 %01%"22"13"40"%01"$%"0224

PAGE 9

$" (A """ "LQ#5O&9D' &99' 5P"LM&9E' "&9D'3C &9E'"" (R99SR9DSR9:S "$#A + 7E9 $#6 L" (&99' / %01%"22"13"40"%01"$%"0224

PAGE 10

( (9T LG&' !" ( LG&"'&9H' ( &' &'LUV&'W?< &'L<&' #$#%& "&9H' $"&15)' "RISR9ESR9HSB15) $""A 1"15)$A & $ & "? "&9H' )"199? "$&9H' 75 "5* ""4"*A !#" "#"5$ (& E %01%"22"13"40"%01"$%"0224

PAGE 11

0@AB 6@ACB 0@CAB 6@AB D60 60 ,E>60 30)4)F8% %01%"22"13"40"%01"$%"0224

PAGE 12

415)$1 0"( R9,SR9.S" & *N & *N* )& *N 4 & & $ !199 &99'" & & 4 O# & 19D 4 $* $? "" 7D9( "(" R9IS < %01%"22"13"40"%01"$%"0224

PAGE 13

19DK 1 ?$7! &9H'"N "3C#(B" & & ) 19E ?" *)+ &9<' /1 &9<'(B M&(9' &>(9O' C @)4*B $ X %01%"22"13"40"%01"$%"0224

PAGE 14

19EK(B5 I %01%"22"13"40"%01"$%"0224

PAGE 15

&(BRB"" 4 & &&Y9' L((F9&G&E8(TPFZ['P #"B 4" "(9&M' " 7EE 4 &9<'" @ &M )&)/+)) &9I' & G L & +)) G & @H++ @ 9 4 : B ( %01%"22"13"40"%01"$%"0224

PAGE 16

*"&9<' IJKLM@NOB@ CJC&&BCJ4PQ@ $-& CJR+ ) &BBS (9 C 0& CJOOC&BCKJ $/$) GT+ 4 (3CK $ T 1$/$/) LG(&9(.' 15) ""A "%"&9.' &9,'7A 9E 15)"A 3C ( "&9H'"5" """ & & %01%"22"13"40"%01"$%"0224

PAGE 17

$A 9H"" & ? "*[B719H "1 *)2 5 / 3$ "K M*5P 43/)# 56 /457 T 458 ( 45 >9 9 :9 B" & & 4 2457 ( 4 *\ Q 45 ;<9 > &=*$/$9*$$*$!'$)/ "&9H'B ($*$!$ C)+G>U+ +"M "M & ) # !K" [M "1 )! %01%"22"13"40"%01"$%"0224

PAGE 18

C 19HK*5 99 &TB' ] %01%"22"13"40"%01"$%"0224

PAGE 19

!M "; $"3C ; "5 ; 5 "A $"" ; 1*""&' "&!'" ""15)$ & '" A "" 1" $A $# 3A &99' UW&'LW&'; R9GS15)$*C 3 )/ %01%"22"13"40"%01"$%"0224

PAGE 20

)4/>2%03%3" ?""" ; ^CR,S"^CRES# R9HS&""' B15)$ /A .&& L & &99G' ?!&99G'$A $ & &RHS'!" ;(7 & &99G'%" "K >&&/.&& (" 9E %01%"22"13"40"%01"$%"0224

PAGE 21

"( .&6&>& &999' K 6& L &2& !(A ? "&999'K .&& &&>& &99D' +! .& $ $&99G'IM_ >& %"&99D' & & & `( &* ?&& &99E' +IO & ? 9H %01%"22"13"40"%01"$%"0224

PAGE 22

(" K 9'0 D'1" E'7"" H'/ <'0 "" "K "%"" $H "% 19< .&.&. H' &' "(A "" "B 4"" 9< %01%"22"13"40"%01"$%"0224

PAGE 23

> > ; ;D a V /; 06L 3"6L 1" 19
PAGE 24

15)$ & $ & "7I19, @& $ & 1994 ^!C!7B &7DD '% & 15)$ .&A& L 4& &99H' A& & &&& 4& L &B & 9 CB 415)$" $ ($" ?A 1 R9HS"A $" & & A 9I %01%"22"13"40"%01"$%"0224

PAGE 25

5&8T9' &8' &T8' = *" 474" 5" 19,K/; 9: %01%"22"13"40"%01"$%"0224

PAGE 26

) b & D & "* $ 019I"A $ D& L ?$ G ?*$?*$?* C 9*C M&'O & D & "* $ ?" $ &) && &19I'K 9 L $9*C $* !$* *! ,C)&UN@)K)*B $*C9! W>&ACGG>/@;,C)&UC;,UAB4 3 1A R9HS A !" "" & D & 3)& b3c & "*MB*" 9. %01%"22"13"40"%01"$%"0224

PAGE 27

; &TBT' &TB' &B' ;D &8' &T8' &8' 19IK1/; /! &TDBT' &TDB' %01%"22"13"40"%01"$%"0224

PAGE 28

& ?"$ & L 3)&& 19IK MM8L TMBTTMMTBT(&Z(9I' )E D & " & ?"$ & F 19:K &$ L9T *$/ T9TdT\eTT9'(&9(9:' 35 1" ?""A K";(7 ;(7?A ((! &;(7!' K K >T ($ TT & (L / 5K 2 9 &/ 9 C ($ CCBC &&// 9'L $ &9 M)5B *K T T9( L /) %01%"22"13"40"%01"$%"0224

PAGE 29

; ;D BT9 19:K1/; // %01%"22"13"40"%01"$%"0224

PAGE 30

#"[ "" ^CRES";( 7A ? ""WE("6!A 5"" "WE(5"? ;(7 5* 5f?!" &;(75' 1" & "(B $) 7 ;(754 5 8 "A 1" & "" & & DE %01%"22"13"40"%01"$%"0224

PAGE 31

"*(B&BOT9' $) A & ) $ (9 ) 5 $) & $ T9' ""$;(7 M! A """*" ""! $ ZT $ (BT9"? (;(7 *" $ #(;(7 ""*" ? A (( #" "[ ""A DH %01%"22"13"40"%01"$%"0224

PAGE 32

)4 0202 B""15)$ "$ (+" 15)$" )4 4)2:12D22 D$ ?$ *CLRQ9SDK 79915)$? 15)#"$ &9DG'(" "A ? 6+6H@=BG,@=B=& >H@=BO G =8 @ ) 4 /! B %( G12@-FBC12@F/-,B @)4/)B D< %01%"22"13"40"%01"$%"0224

PAGE 33

B ) L B T'&&I['T 12@1/-A,BB& 4) LG &9DG'&9D9' b"A &9D9' "$ B D "A 15) ?1 15) @ 2)F6 OM` 15)$ D @ !1 A D @ ? "B ) 1 /* %01%"22"13"40"%01"$%"0224

PAGE 34

15) "P >)# "? E #! E $ G6 4"515) "&& "'DT& "'D'W" ($) ""& &"Q5'&"Q5'"&"Q5'$ "15)'1 H ,H 15) K) B 6 F6 IM 6 6 99 HG9)(D D
PAGE 35

13 dD 9 9HG9)(DE:.)(9DE:.)(9D 9 H99I)(9E<,)(D :.,)(E 9 : DED)(9 9
PAGE 36

""(" 9!"15) 6! 15) #" K15) """ $15) "1"15) 7" "+"15) D. %01%"22"13"40"%01"$%"0224

PAGE 37

)4 4/:21002%21:02? 12 *C """ #"B ? !"A &' ;(7 1 Q5P#5O&'LM&'9I 7 ;(7" 9 6RES ;(7A %" EG %01%"22"13"40"%01"$%"0224

PAGE 38

"R.S#" ( !&A "" %! "" 79D 1 #" ( R.S 6 B+" "" B&' ?( E9 %01%"22"13"40"%01"$%"0224

PAGE 39

" / F9T ) LG F9(8( LG T LM @=BGX-@=B=8HH2H2 0@=BG =* #$ X5"4 1 9.4 &DD' M 4 $ ? ""& ;7$ /"1 :: GH 9,9, 9G EDED 9: ,H,H DG #( (?" $! ED %01%"22"13"40"%01"$%"0224

PAGE 40

19.K*3 EE %01%"22"13"40"%01"$%"0224

PAGE 41

/15)(A 15)" # RQ9SD "K GG< bGbG GG9 KUGbG 9GG $UGUG EEEE KbGUG "9 G&' 15)D<,D<,f" 15) $" 15) ;7$ F6 /F6 DD ED,)(9 .
PAGE 42

13 ;7$ ( DD:D:)(9 .
PAGE 43

B&'1 "" K" B +" [ !"" !" "B &99<' 199G LG[ ^C ) ($[ ^"C ^"C $" ^C R9S !R:S" E, %01%"22"13"40"%01"$%"0224

PAGE 44

g7 LL/ LLG 199GK^6"C^3C/ EI %01%"22"13"40"%01"$%"0224

PAGE 45

(" ?(" "/ "15)$( #(9"T5OLG >O6G,4 ?15) + I 3C 4>& $& & 4> $ 3 & 9 '(&'(&^'&h'h Ueib "A &9DD'3 G>KZ f D&/ i >'& < ('(&' E: %01%"22"13"40"%01"$%"0224

PAGE 46

415)$ ( ^C" ^C" &9DH'&9D<' &9D<' I &9D<'""?" &I!< & 4>&& TW'&9D,' &9D<' ( $&I!4>& & Lj $&I!> W&9(DI' B!" %!" &9D,'" 15) ("!R9
PAGE 47

I R9:S +Y60202 ? (" "G99GG B" ("" 15)!"^C %K !"^CVR9
PAGE 48

;7$ /"1 9,9, HE EDEDHH ,H,H H, ?" *!( "" $%"f-" &"(" 9' 15)$ ( 79ED?fI 15) I (""!#" (""" "15) (" 79D!" 15) 79D (" H9 %01%"22"13"40"%01"$%"0224

PAGE 49

)4 4'02:0212 A "15)3 7E9 # ZLRQ9SD" Og=O9G= J%#0 QG= OQ + + OGGG9OG + C + M &"LG' D$' 'L 7) &DV ) 15)$""$ "" + A &"' %%& ;7$ F F6'F6 (=(!4'(!' .45'/+ 4/)*+ )*=)* !4/.!) /4)'!+ )4!5)+ /= / !4)/* .4-!(+'/4(*(+' *'=*' !4!* )4''5+' -4/./+. '/ %01%"22"13"40"%01"$%"0224

PAGE 50

0 L9
PAGE 51

' RDS+ # LRGSD"" GUKU</A@/, A GUKU<
PAGE 52

15)$ H< %01%"22"13"40"%01"$%"0224

PAGE 53

/890%8%2 /4)"=%>:217 $799 ?"15)$ 1D9# I M +, =&& \ 5#&'O&'P N LG O # ["" #3C #F9"&'T5UW&'LG&DD' 5P"&'LM&'&DE' H, %01%"22"13"40"%01"$%"0224

PAGE 54

&TB' T&B' 1D9K;@; HI %01%"22"13"40"%01"$%"0224

PAGE 55

?15)$ $15)$799 "$ / *"&>b'&>'&>'&>' 1DD 6* 3 K@KN K!! C@KZKBKC@KZK C@KZKZKCKBK ,A @K] $$ C $ BKC@ZBC@ZZCBK+ 3 ""7" "&>'N 3 &>'N*"" 71DD""A $K&'N* 2&>'"&9G'* k&>'"&G9'* "k&>' "2&' $7 )4) 10"( &R
PAGE 56

& 1 / / K;@ '5 %01%"22"13"40"%01"$%"0224

PAGE 57

" "A #" B $ K [ ? 3 5 P8VK__k__" "P -B__2__" #1DE ?B" 7A !1D9% *8 &$ =T ) B &$& VUC) " 7 )4) & & "K
PAGE 58

1DEK 1 7 )4) & & & ?" "&DE' &DE'"* #" <9 %01%"22"13"40"%01"$%"0224

PAGE 59

(B *B$" &' K &?MA /B*FB *B8(&9MD ) '&G9' *(B&G ) M / '&G9'*B 1DH"B ?"^"C*(B*B "?"" ^$C1B*B(*BB"5[B & ) '& ) M / ) 3$/$ &G9'&G9MD'*(B ?$799 3C$ ""? $"OP"3$) "3C$"
PAGE 60

1DHK5
PAGE 61

!(3C"" "!3 "&9G'"" &9I'" *" & *"&9:'" & 1 "3C "" "" & ;"" 1"2? 3 "1DH?! &DD'M2& ) "^C B7!&DD'2& ) "^C[B+ J """ & ? 7ED[ "4 !5"
PAGE 62

!3C"k1 "1DH3 K /9T J$ TE`TB *@D'$/$ C @#'$/$*$ C @0'$ C @>'$*C CZG P B "/ " 4 7E9 ^CH/ $" # /4/0662 /4/4)>2%03% 79D 15)$ $A ;"A "" << %01%"22"13"40"%01"$%"0224

PAGE 63

"1D< ) D%9 Q 15)$" 15)$ Q "A 4A 7A 9D# A 7EE +" "7 79D9MD #A <, %01%"22"13"40"%01"$%"0224

PAGE 64

3) 3M 1D
PAGE 65

" &*" *"' "A 3C #"" &'1D, "7 15) ""["" &_'1A "" 15) 4"A ("( <: %01%"22"13"40"%01"$%"0224

PAGE 66

23%:1>0 3%:1>0 1D,K/135 <. %01%"22"13"40"%01"$%"0224

PAGE 67

"" 15)$ " /4/4/+Y61" "B &' K K "BA ("79EE" ("A I 3C 1 I &9.'C ,G %01%"22"13"40"%01"$%"0224

PAGE 68

"" ($^C" ^C"15)$?! +"A "!15)$ ( "( $ +"" B"$^C" ^C"15)$ "! ^C"$ ^C"1&9.' "" I I A "! I I ("! ,9 %01%"22"13"40"%01"$%"0224

PAGE 69

/4 0202 /4 4)03%D1 7DD91 (( 115) ? ( ( LRQ ) ) SR( ) ) S" QD & D & $ & D< 1DI ?"15)$*C 7A DD"" ""A ;7$ /"1 9,9, 4)/ EDED 9I ,H,H 4// +";(7" A "$ ,D %01%"22"13"40"%01"$%"0224

PAGE 70

( !4/ (GH ( !4* ( !4( 1DIK;DGl3 ,E %01%"22"13"40"%01"$%"0224

PAGE 71

"A "! 1" ( ( "" Q< & R& 1DI 15) ;7$ /"1 9,9, GE EDED G< ,H,H G: A K 15)A ? "? "";(7 "" ,H %01%"22"13"40"%01"$%"0224

PAGE 72

(9(G:(G,(GH(GDGGDGHG,G:9 1D:K;
PAGE 73

( ( DGl% "G99GGO 15) ;7$/"1 9,9,GH EDED G, ,H,H G: # #" 15)A 7DDDA ""A +""" 15)A "1 K ( (
PAGE 74

,H,H' /"1 ) G: )! 9E )!! 9< )!!! 4/! 1A 15) & ? & '" !+" /4 4/+Y6D1 ("7DDD "15)1 ,I %01%"22"13"40"%01"$%"0224

PAGE 75

"^C 1"!"^C 15)$*CLRQ ) 9SRQ99S DGl"4 !15) #"" ("$ ;7$ /"1 9,9, GI EDED G: ,H,H G: "!" ^C"#7DDD I 3^C"! ? "!
PAGE 76

I " "15)$ 4!) 9GG"" ;7$ /"1 9,9, HE EDED HG ,H,H E: ?"*(!A (" 15)! 15) !("""A ("? Kf-" 15)C C# ,. %01%"22"13"40"%01"$%"0224

PAGE 77

/4 4 1101"=%>:217 15)$ 1 $ "?*C G GI+ 9 & 9 S/F 6C>X^G!& 5("LM "P > L 4 1M 4 ) L&I'&DI'T 15)$$A LG 1D. "MV15) $""MV? 79E @ E """ 15)1MVL:GcIGc" ;(7"""1 IG %01%"22"13"40"%01"$%"0224

PAGE 78

@+ 9 & 9 B@ 9 & 9 B @+)&!B & ! &9G' @+ 9 &+ 9 B @ 9 &+ 9 B 1D.K/; I9 %01%"22"13"40"%01"$%"0224

PAGE 79

MVL,Gc
PAGE 80

"" 15)$ ?$ 9 +"!" 4 " [" #" S O =M ) =O 02 0 = Qg 0 OGGG9 g + C L H
PAGE 81

G< % %1;7$ 6 6 6 9,9,EG.D H
PAGE 82

&D<9' &99' "2h @+ 9 &+ 9 B &(<(9' @ 9 &+ 9 B 1D9GK*3 I< %01%"22"13"40"%01"$%"0224

PAGE 83

DD'( 15) ;7$ 6 ? '6 '=' !4/./' -4!*.+ (4/*5+ (=(!4) !/ )4-'-+ /4)5-+ )*=)* !4!*.( '4.*-+' .4(.5+' /= / !4! )4/) +')4.(/+' *'=*' !4!)-' 4/ '+.'4/5'+. +(" ("" *C I, %01%"22"13"40"%01"$%"0224

PAGE 84

1"1:62 15)" 4)8%2%2D$ 79915)$A #&99'+ # "& '? U "!CUJ' U 3CA &9D'! "LQ#5UW II %01%"22"13"40"%01"$%"0224

PAGE 85

#(9"T G 4""L ) "(K &#'(9 J2U)* LG&' `&#'(9 2U T )* LG&"'&E9' LM&' 15)$ K"A 799 &E9'&9H'A ""# "7EE #" & $) 1E9 ?" M&&#'(9&Q U) T J$) G&ED' "" 799 "" O(A & I: %01%"22"13"40"%01"$%"0224

PAGE 86

1E9K5 I. %01%"22"13"40"%01"$%"0224

PAGE 87

"" / 799&ED' I&&#'(9' $&*0 T -* ) !(6")/C) C*C*'<$ TT ) *&! gF\P' ?"5" 3$ 799 4/:6>128 90%2 4/4)1:021 +""A 3? "7EDD71ED "$3C :G %01%"22"13"40"%01"$%"0224

PAGE 88

1EDK5 :9 %01%"22"13"40"%01"$%"0224

PAGE 89

*[**" "3 j9"&'T5U&'LG&EE' #7D9!&EE' J$,) "^C"7 !&EE' @>,) ^C"? 3"6M ^"C""K M&(9"& ) T5O& )) P J,)* > &(9"& )*')) P @>,) L &EH' J$ "&K ) P2&[ )*J$'' ) P2& )* M`h"&[b )V &[b ) T @>>') e2&K ) """"" "" 1"""1 M"?" &EH'/" $')/,) $ 9 "G :D %01%"22"13"40"%01"$%"0224

PAGE 90

K & )W '__k&K'__" ) PIMB&KM'__2&K'__" &E<' D,& &K'PI&>i'__k&>'__ L ( 7 ( $ D,& D,& & ) P $N '__2&K'__LG ) 2?" )/,) L__&'____2&'__ )/X) __&'__&G&'' + 2k1ED7" __&['__L&&['P N)) M&G&[''&E,' :E %01%"22"13"40"%01"$%"0224

PAGE 91

*" MF ) &BM'P ,) LMGY9h&Ki'&Zii'^>i &EI' LMc&&& ) P N )) M&G& ))) __2& ) L YY!,) 8 UG__`____a__W4 ?" ,) !" &9G' ) 2 F T&m ( $'mT&m(m )? T&m(m ( mTm ) ) LDMGGT&DM9G(DMGG'mT&(DMGG'DMT&DM99(DM9G(DMG9T'm ?"2K 2L ) &&m9GQm'T&mQm9GQmG9Tm'& ( ( &DM9GQ / M'T&DM99j / MQDMG9TDMGG'DM' &EI' "*M 4" M !
PAGE 92

9 G&E<' DA& &'Pl&m'__2&B'__L &E<' $ Db& D,& ) P ->&> / M'__k&K'__LG > k?" &$'M__&['__Lk&VM'M__k&KBM'__ &>'P ) L__&'__&G&'' 7" &>i'L&> )< l&$b ''M&G& )) &E.' *" !7)/,) M&F9__&'__&&VM'M__&K'__'P ,) LM @c) &&&m ) PBB&B )) M&O&K & )! P&k&m ) '__'P ,) Y/C,/ :< %01%"22"13"40"%01"$%"0224

PAGE 93

#" ""k K kL ) &&>(> !! 'T @K)) ( K)! ( K!) T> !! '>&E99' & /,!) Q / M'T& /,)) Q / MQ /,!) T / M !! '>'P &E9G' "*M " "$ +) &$ ) P2& ) *M* ""3C """"" 4""5UO ')/,) M (* %01%"22"13"40"%01"$%"0224

PAGE 94

?2" 85> )W 2&$ ) LM& ) P& ) L YI$ T M &E9D' L H L +88 < 4"" )*U*J *M "?&E9D' ) P ,) YY; 8MB /0 1 L YIZ$/Y&FU)/U ?""" "*[ :I %01%"22"13"40"%01"$%"0224

PAGE 95

"UW[" "*M" ) *&L +"*M !"8 1ED" > ?&>i ) P&>b ) #2331 "" *[* M J ? """ ) ML P) L P[) \ ) )<,)*J> ')/,) L P$*&/$$) +_DD_*M* (( %01%"22"13"40"%01"$%"0224

PAGE 96

4/4/2:021 1" ""UO?" # +) & )W 2& ) #"" ?" 1ED"" ? 9"G ?K "!C7 'P2&K ) / M'P "!,) LM&>i'&&>i'M&>i''P "!,) @@;d ON,K@Ke BB M&G& & ) P&2&'M__2& __'P#F2& ) L ,) P# !,) M8M#( ) 2&'P ,) B[( :. %01%"22"13"40"%01"$%"0224

PAGE 97

&E,'2[ "[ B?"A &E9E'&E:' 2""A B ? ) ?K M&# +) &K ) P2& ) L 7)<"!,) GU0& &>'&&>ii'M&>b''>#F2&K ) GU0&@@@Ke W' i''M&G&'''&E9H' P&k&'M__k&>'__'P "!,) L ) P "!,) "/<) P ,) &E.'k# " .G %01%"22"13"40"%01"$%"0224

PAGE 98

4"# # +) k(62!"! # " 4 238%328 28%32 # ? #7DD [!1EE "* J "*( $ L ) DEH 6? J "$ J U "$*( "K*" ?LO&BTD'TO& TH'&E9<' ?"" .9 %01%"22"13"40"%01"$%"0224

PAGE 99

1EEK/15 .D %01%"22"13"40"%01"$%"0224

PAGE 100

"K ) '$* J/ D' 3J E'" J) H' ; PVB__2__" &E9,' L9DEH$* J/ ?"&E9,'K .4 @/ &TD'T Ð''P 5# _Ue@O 6# C, & / O--B # L9&9'Tn&Z'L9( "$ ? / "L ) *B*D" &E9,' ?&E9,'" K )P *\2&m 'LD2&K ) &E9I' .E %01%"22"13"40"%01"$%"0224

PAGE 101

+2&9G' *2( *(0&E9I' ?"&E9,'1 1EE?"K ?&K'P $X '__2&mM'__ L&W &K'T#O''PO&mb'&D__2&>'__' &E9:' L& &$b / M'P 5 i'__2&m'__T& $) PBM>'__2&>'__ LJ 9 "* L 7 +& ) *D"&E9I' & / $* Pb__ "" 1EE" K ?& ) P '__2&'__ L&? &m'O&K''D__2 ) &m'__&E9.' L &>i'(O_2&i'__L ?9J .H %01%"22"13"40"%01"$%"0224

PAGE 102

"&[ ) * 0&E9:'&E9.'" &E9,' &E9<' ? #3 "* 1EE K M# +) ?&m ) P2& ) *D*H ? ) *K / K ?L ? T_ / `1G/`/ K L _' 2L / 2 .< %01%"22"13"40"%01"$%"0224

PAGE 103

K ?O? 'P2& ) 5 #( ) ZJ) P ,J) *5$"!CZJ)/,J) Lo / #j ) &'e2 / &' *5"!C7) P ,) *T5D"/C7D)],D) #"A 3 3 1# !# ., %01%"22"13"40"%01"$%"0224

PAGE 104

'0%00;< 15)$" "%A 15)A (! (15) A "(" """ "" "15)#15) ""! 5%01%"22"13"40"%01"$%"0224

PAGE 105

" "$B RES?$ $ ?" "" 15)$ ?15) "! *! 15)A (""#A 15) 1" 15)!A "#15)" .: %01%"22"13"40"%01"$%"0224

PAGE 106

"A ( .. %01%"22"13"40"%01"$%"0224

PAGE 107

$3" R9S0)## 8)38?*A A Y"IYJ$Y@%F5 DKHEG(H%J66$K4#%&6 C9&Y3PYF5%#$F%KO6#6>%$>Y$F$%K <9(<.9..< RES# I$4>$6J&K$^6#CL_D4$67$&55$J$%K#%$ KF$J# ;3(74:<;! 3"7#9.:H RHS# 71/!8?0# &#;'38)" Y5>#$ K#"55$J$%K# /"*/9.:H R$`$K$6P6F6KI6&%# "49< Y5>$K46>Y6>$6#$K@%F5$%KI&6F$J# 7(59..9 R S?6 "I$4>$%>$ 7# !*9.:I RISp/871/!" Y"IQ/ F6>"K D:KE.D(HGD9..9 R ( S8)3 @%F53&# H:KE,,(E:,9.:D R.S8)371/!8?01077 (" 3>%J66$K4#%&6 C+& Y3PYF5%#$F%KO6#6>%$>Y$F$%K 9.(D<9.:. )!! %01%"22"13"40"%01"$%"0224

PAGE 108

R9GS/388#*"( (" 43#71/! Y$& @%556>I%K$K@%K6>6KJ6%KI$4>$I6&%# "/*EDDH 9<.(9IE+5#9..E4#7# R99S0))01+ A 3>%J66$K4#%&6&Y3PYF5%#$F%KO6#6>%$>Y$F$%K EH9(EHI9.:E R9DS0))101?7A A 3>%J66$K4#%&6&Y3PYF5%#$F%KO6#6>%$>Y$F$%K I9(:99.:E R9ES/671/!"&15)' 3>%J66$K4#%&6CC&K6>K$%K@%K6>a 6KJ6%K@`b ?5#9.:: R9HS71/! I$66"5$6I6&%#%>3>$b$6>6K$ P^$%K# ", `>%K$6>#$K"55$6I&6F$J# 7# *9.:. R9%$>Y$F$%K E,(9G,7# 9.:E R9.S8 Y>VK#6KF6>$^66#F6&%6#V66F6K#$K$# &U>$6#6F$6# *3C"*/ 9.II 9G9 %01%"22"13"40"%01"$%"0224