Citation

## Material Information

Title:
Discontinuous enrichment methods for computational fluid dynamics
Creator:
Oliveira, Saulo Pomponet
Publication Date:
Language:
English
Physical Description:
146 leaves : ; 28 cm

## Thesis/Dissertation Information

Degree:
Doctorate ( Doctor of Philosophy)
Degree Grantor:
Degree Divisions:
Department of Mathematical and Statistical Sciences, CU Denver
Degree Disciplines:
Applied Mathematics
Committee Chair:
Franca, Leopoldo P.
Committee Co-Chair:
Bennethum, Lynn S.
Committee Members:
Briggs, William L.
Knyazev, Andrew
Trapp, John A.

## Subjects

Subjects / Keywords:
Discontinuous functions ( lcsh )
Fluid dynamics -- Mathematics ( lcsh )
Stokes equations ( lcsh )
Diffusion -- Mathematical models ( lcsh )
Heat -- Convection -- Mathematical models ( lcsh )
Diffusion -- Mathematical models ( fast )
Discontinuous functions ( fast )
Fluid dynamics -- Mathematics ( fast )
Heat -- Convection -- Mathematical models ( fast )
Stokes equations ( fast )
Genre:
bibliography ( marcgt )
theses ( marcgt )
non-fiction ( marcgt )

## Notes

Bibliography:
Includes bibliographical references (leaves 132-146).
General Note:
Department of Mathematical and Statistical Sciences
Statement of Responsibility:
by Saulo Pomponet Oliveira.

## Record Information

Source Institution:
Holding Location:
|Auraria Library
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
53993875 ( OCLC )
ocm53993875
Classification:
LD1190.L622 2002d O54 ( lcc )

Full Text
DISCONTINUOUS ENRICHMENT METHODS FOR COMPUTATIONAL
FLUID DYNAMICS
by
Saulo Pomponet Oliveira
B.S., State University of Campinas, 1997
M.Sc., State University of Campinas, 1998
A thesis submitted to the
in partial fulfillment
of the requirements for the degree of
Doctor of Philosophy
Applied Mathematics
2002

This thesis for the Doctor of Philosophy
degree by
Saulo Pomponet Oliveira
has been approved
by
William L. Briggs ^
Andrew Knyazev
Dec. /*?. apco-
Dat^

Oliveira, Saulo Pomponet (Ph.D., Applied Mathematics)
Discontinuous Enrichment Methods for Computational Fluid Dynamics
Thesis directed by Professor Leopoldo P. Franca
ABSTRACT
This dissertation develops and analyzes discontinuous enrichment methods
Stokes problem. These methods fit into the class of hybrid finite elements, where
continuity of the unknowns are enforced with Lagrange multipliers. The flexibil-
ity of employing globally discontinuous weighting functions allows us to choose
fundamental solutions of the differential equation as shape functions (the enrich-
ment functions), in analogy with the method of residual-free bubbles (RFB). The
enrichment functions derived herein were crucial in capturing boundary-layers on
advection-diffusion problems and stabilizing a low-order velocity-pressure pair
for the Stokes problem. Numerical experiments contrast DEM with stabilized
continuous and discontinuous Galerkin methods.
in

This abstract accurately represents the content of the candidates thesis. I
recommend its publication.
IV

ACKNOWLEDGMENT
I would like to thank my advisor, Leo Franca, for his support and guidance
during all stages of my doctoral journey.
I am grateful to Lynn Bennethum, Bill Briggs, Andrew Knyazev, Tom Rus-
sell and John Trapp for kindly joining my committee, Isaac Harari, Charbel
and Cristina Cunha, Che Sun and J. Tinsley Oden for providing valuable refer-
ences.
A significant part of this dissertation was developed in the National Labo-
ratory for Scientific Computing (LNCC), Brazil.
I would also like to acknowledge and thank the financial support of CNPq,
the Organization of American States and the Department of Mathematics at
the University of Colorado at Denver.

CONTENTS
Figures .............................................................. ix
Chapter
1. Introduction........................................................ 1
2. Hybrid Methods for the Poisson Problem........................... 6
2.1 Functional Settings............................................... 6
2.2 The Primal Hybrid Formulation.................................... 10
2.3 Mixed and Hybrid Methods......................................... 13
2.4 The Reference Element............................................ 17
2.4.1 Partial Derivatives in the Reference Element................... 20
2.4.2 The Laplacian ................................................. 22
2.5 Primal Hybrid Method............................................. 24
2.6 Other Discontinuous Formulations................................. 27
2.7 Mortar Element Method............................................ 28
2.8 Discontinuous Enrichment Method.................................. 31
3.1 Finite Element Methods........................................... 36
3.1.1 Discontinuous Galerkin......................................... 38
3.1.2 Stabilized DG Method (SDG)..................................... 39
3.1.3 The Mortar Element Method...................................... 40
3.2 Discontinuous Enrichment Method.................................. 41
vx

3.2.1 The Subspace .................................................. 42
3.2.2 The Subspace VÂ£................................................ 43
3.2.3 Computation of JV|,............................................ 45
3.2.4 Domain Splitting / Computing VIV|, ............................ 47
3.2.5 Existence and Uniqueness of Solutions.......................... 49
3.3 Numerical Experiments.......................................... 55
3.3.1 Advection Skew to the Mesh..................................... 55
3.3.2 Unit Force..................................................... 57
3.3.3 Elements with Three Inflow Edges .............................. 61
4.1 Finite Element Methods......................................... 63
4.1.1 Stabilized Methods............................................. 64
4.1.2 Stabilized Baumann-Oden Method................................. 65
4.2 Discontinuous Enrichment Method................................ 67
4.3 Matched Asymptotic Expansion in the Unit Square................ 68
4.4 Boundary-layer Functions....................................... 72
4.5 A Modified Gaussian Quadrature................................. 73
4.6 Numerical Experiments on Rectangular Domains................... 74
4.6.1 Advection Skew to the Mesh..................................... 75
4.6.2 Unit Force..................................................... 77
4.6.3 Advection in a Rotating Flow Field............................. 80
4.6.4 Thermal Boundary Layer ........................................ 81
4.7 Asymptotic Expansion in Quadrilaterals......................... 81
vii

4.8 Generalized Boundary-layer Functions............................... 84
4.9 Numerical Experiments on Quadrilateral Meshes................... 85
4.10 Limitations of DEM ............................................... 86
5. Stokes Problem....................................................... 94
5.1 Finite Element Methods............................................. 95
5.1.1 Stabilized Finite Elements....................................... 96
5.2 Stabilized Baumann-Oden Method..................................... 97
5.3 Discontinuous Enrichment Method.................................... 99
5.3.1 Free-Space Solutions of the Stokes Equations ................. 100
5.3.2 Spaces Qf and V* ............................................... 103
5.3.3 Stability of the Pair X Qh...................................... 103
5.4 Numerical Experiments............................................. 118
5.4.1 Smooth Velocity-pressure Field.................................. 119
5.4.2 Driven Cavity Flow.............................................. 121
5.4.3 The Stick-slip Problem.......................................... 123
6. Conclusions......................................................... 130
References............................................................. 132
vm

FIGURES
Figure
2.1 A representative element in a partition of fl........................ 17
2.2 The reference element................................................ 18
2.3 Mortars and nonmortars of a subdomain interface...................... 29
3.1 Inflow, outflow and tangential boundaries of K....................... 35
3.2 Traces of Vh from the inflow (v^) and outflow (v^) .................. 37
3.3 Mortar and non-mortar sides of a 2x2 uniform mesh ................... 40
3.4 Approximate solutions of Prob. 3.3.1 (Sec. 3.3) by DG (left) and
mortar (right) methods ............................................. 44
3.5 Location of xp....................................................... 46
3.6 Edges that do not contain xp. From the left to the right, we have
t > 1 in the first case, s > 0 in the second, and A = 0 in the third . 47
3.7 One inflow vertex.................................................... 48
3.8 Two inflow vertices.................................................. 48
3.9 The points x+ and xc used to evaluate N^ x........................... 49
3.10 A partition satisfying assumptions 3.1-3.3.......................... 50
3.11 A sketch of the function vfc ....................................... 53
3.12 Meshes Ml (left) and M2 (right), n = 20............................. 56
3.13 Average CPU time of DEM, SDG, and DG versus n ................. 56
3.14 Statement and reference solution of Prob. 3.3.1, 0 = 30............ 57
3.15 Error plots of Prob. 3.3.1 (30), meshes Ml (left) and M2 (right) . 57
IX

58
3.16 Error plots of Prob. 3.3.1 (45), meshes Ml (left) and M2 (right) .
3.17 Error plots of Prob. 3.3.1 (60), meshes Ml (left) and M2 (right) . 58
3.18 DG solution for Prob. 3.3.1: uJL (left) and (right)........................ 58
3.19 SDG solution for Prob. 3.3.1: % (left) and uÂ£ (right) ........................ 59
3.20 DEM solution for Prob. 3.3.1: u(left) and (right).......................... 59
3.21 Statement and reference solution of Prob. 3.3.2, 6 = 30................... 59
3.22 Error plots of Prob. 3.3.2 (30), meshes Ml (left) and M2 (right) . 60
3.23 Error plots of Prob. 3.3.2 (45), meshes Ml (left) and M2 (right) . 60
3.24 Error plots of Prob. 3.3.2 (60), meshes Ml (left) and M2 (right) . 60
3.25 Error plots of Probs. 3.3.1 (left) / 3.3.2 (right), mesh M2 (15) . . 61
4.1 DEM solutions for Prob. 4.6.1 (30) with Neumann (left) and Dirich-
let (right) outflow boundary conditions.................................... 68
4.2 A domain with (^-boundary ................................................. 69
4.3 and gxp for 1 < 0l/ < 103............................................. 74
4.4 Average CPU time of DEM, SBO, and GLS versus n............................. 75
4.5 DEM elevation plots for Prob. 4.6.1, n = 40, 9 = 30, k = 10-m,
m = 0,3,6 (left to right).................................................. 76
4.6 SBO elevation plots for Prob. 4.6.1, n = 40, 9 = 30, k = 10~m,
m = 0,3,6 (left to right).................................................. 76
4.7 GLS elevation plots for Prob. 4.6.1, n = 40, 6 = 30, k = 10m,
m = 0,3,6 (left to right).................................................. 76
4.8 Error plots of Prob. 4.6.1, 6 = 30, 45.................................... 78
4.9 Error plots of Prob. 4.6.1, 6 = 60 ........................................ 78
4.10 Error plots of Prob. 4.6.2, 9 = 30, 45.................................... 78
x

4.11 DEM elevation plots for Prob. 4.6.2, n = 40, 9 30, k = 10 m,
m = 0,3,6 (left to right)............................................ 79
4.12 SBO elevation plots for Prob. 4.6.2, n = 40, 9 = 30, k = 10-m,
m = 0,3,6 (left to right)............................................ 79
4.13 GLS elevation plots for Prob. 4.6.2, n = 40, 9 = 30, k 10-m,
m = 0,3,6 (left to right)............................................ 79
4.14 Statement of Problem 4.6.3.......................................... 80
4.15 DEM, SBO and GLS elevation plots for Prob. 4.6.3 (left to right) . 80
4.16 Error plot of Prob. 4.6.3.............................................. 81
4.17 Statement and mesh pattern for Prob. 4.6.4.......................... 81
4.18 DEM, SBO and GLS elevation plots for Prob. 4.6.4 (left to right) . 82
4.19 Error plot of Prob. 4.6.4........................................... 82
4.20 Domains D2 and 03 with respective meshes when n = 20.......... 86
4.21 DEM elevation plots for Prob. 4.6.2 (D3), 8 = 30, 45, 60 (left to
right) ................................................................. 87
4.22 SBO elevation plots for Prob. 4.6.2 (f23), 8 = 30, 45, 60 (left to
right) ................................................................ 87
4.23 GLS elevation plots for Prob. 4.6.2 (D3), 9 = 30, 45, 60 (left to
right) ................................................................. 87
4.24 DEM contour plots for Prob. 4.6.2 (f23), 6 = 30, 45, 60 (left to
right) ................................................................. 88
4.25 SBO contour plots for Prob. 4.6.2 (fi3), 9 = 30, 45, 60 (left to
right) ................................................................ 88
4.26 GLS contour plots for Prob. 4.6.2 (fi3), 9 = 30, 45, 60 (left to
right) ................................................................. 88
4.27 Error plots of Prob. 4.6.2 (30), meshes M2 (left) and M3 (right) . 89
xi

4.28 Error plots of Prob. 4.6.2 (45), meshes M2 (left) and M3 (right) . 89
4.29 Error plots of Prob. 4.6.2 (60), meshes M2 (left) and M3 (right) . 89
4.30 Statement of Problem 4.10.............................................. 90
4.31 Velocity field a used in Prob. 4.10, which was found by solving the
incompressible, steady-state Navier-Stokes equations with Re = 60 90
4.32 Boundary conditions of the velocity field a............................ 90
4.33 DEM elevation plot for Prob. 4.10, k = 10............................. 91
4.34 SBO elevation plot for Prob. 4.10, k = 10........................... 91
4.35 GLS elevation plot for Prob. 4.10, k = 10........................... 91
4.36 DEM elevation plot for Prob. 4.10, k 10-3......................... 92
4.37 SBO elevation plot for Prob. 4.10, k = 10-3.......................... 92
4.38 GLS elevation plot for Prob. 4.10, k = 10-3.......................... 92
4.39 DEM elevation plot for Prob. 4.10, k = 10-6......................... 93
4.40 SBO elevation plot for Prob. 4.10, k = 10-6.......................... 93
4.41 GLS elevation plot for Prob. 4.10, k = 10-6.......................... 93
5.1 Two elements of %, satisfying d(Q\ QJ) = 4........................... 105
5.2 Functions |^i,..., Ui \ka (left ro right) in a patch of elements . 113
5.3 Nodal values of jV,-|/q, VV3, ^23,f and ^2?^ (left ro right) ......... 113
5.4 Meshes Ml (left) and M4 (right) for Q = (1, l)2...................... 118
5.5 Average CPU time of DEM, SBO and GLS versus n......................... 119
5.6 Pressure errors of Prob. 5.4.1, meshes Ml (left) and M4 (right) . 120
5.7 Velocity errors (L2) of Prob. 5.4.1, meshes Ml (left) and M4 (right) 120
5.8 Velocity errors of Prob. 5.4.1, meshes Ml (left) and M4 (right) . 120
5.9 Statement of Prob. 5.4.2.............................................. 121
Xll

5.10 Velocity vectors (left), pressure contours (center) and pressure ele-
vations (right) of the reference solution for Prob. 5.4.2......... 121
5.11 Pressure errors of Prob. 5.4.2, meshes Ml (left) and M4 (right) . 122
5.12 Velocity errors (L2) of Prob. 5.4.2, meshes Ml (left) and M4 (right) 122
5.13 Velocity errors of Prob. 5.4.2, meshes Ml (left) and M4 (right) . 122
5.14 DEM, SBO and GLS pressure elevations for Prob. 5.4.2 (left to
right), Mesh M4, n = 20........................................... 123
5.15 Flow leaving a pipe (left) and stick-slip approximation (right) ... 123
5.16 Statement of Prob. 5.4.2............................................ 124
5.17 Statement of Prob. 5.4.2 accounting for symmetry ................. 125
5.18 Velocity vectors of the reference solution for Prob. 5.4.3........ 126
5.19 Pressure contours of the reference solution for Prob. 5.4.3....... 126
5.20 Pressure elevations of the reference solution for Prob. 5.4.3..... 126
5.21 Horizontal component of the velocity along y = 1. A zoom over
x = 0 is shown on the right....................................... 127
5.22 DEM pressure elevations for Prob. 5.4.3, Mesh M4, n = 10.......... 127
5.23 SBO pressure elevations for Prob. 5.4.3, Mesh M4, n = 10......... 128
5.24 GLS pressure elevations for Prob. 5.4.3, Mesh M4, n = 10......... 128
5.25 Pressure errors of Prob. 5.4.3, meshes Ml (left) and M4 (right) . 129
5.26 Velocity errors (L2) of Prob. 5.4.3, meshes Ml (left) and M4 (right) 129
5.27 Velocity errors of Prob. 5.4.3, meshes Ml (left) and M4 (right)
Xlll
129

1. Introduction
The Navier-Stokes system of equations [75, 91] is an established continuous
model for Newtonian fluids, in particular incompressible ones (see [16] for a
discussion on the validity of this model). Its dimensionless version states that
the velocity u of an incompressible fluid and the pressure p satisfy
du 1 _ , . ,
+ u Vw + Vp V e(u) = /
at He
V u = 0 ,
in addition to initial and boundary conditions. The function / represents nor-
malized external forces and the constant Re, known as the Reynolds number,
provides a measure of the ratio between inertia and viscous forces. Moreover,
e(w) = (Vit + (Vu)f)/2 represents the deformation rate tensor [34].
Unfortunately, exact solutions of Navier-Stokes are known only for a few
limited cases [17], so one cannot count yet on a general formula for the velocity
field. Experimental simulations are generally expensive, while the exact solution
of a simplified model may result in inaccurate predictions. An alternative is
(CFD), the simulation of flow fields by numerically
solving the differential equation that models the fluid. CFD also contributes to
theoretical and experimental fluid dynamics by providing clues on the properties
of real fluids or on the qualitative behavior of exact solutions of Navier-Stokes
and other fluid equations.
One of the most popular techniques in CFD is the finite element method
(FEM), mainly because of its flexibility with the domain geometry combined
computational fluid dynamics
1
!
i
i

with a solid mathematical background [57, 64]. FEM generally partitions the
domain of a problem into a finite number of regions (elements) and constructs
approximate solutions from piecewise polynomial spaces defined on these ele-
ments.
Many CFD problems pose several difficulties to the finite element method,
as well as other numerical techniques. Some of these difficulties are:
The incompressibility condition, which constrains the choice of the finite
element spaces used for the velocity and the pressure [18].
Propagation of discontinuities and shocks, typical for first-order hyper-
bolic problems. They may produce oscillations on continuous approximate
solutions [64].
The development of boundary layers, sharp variations of the solution near
the domain boundary, which may cause oscillations in numerical solutions
as well [86].
In general, a standard (Galerkin) finite element solution cannot capture
aspects of the exact solution which are restricted to regions that are small com-
pared with element sizes (the fine scales [58]). Some finite element methods
have been developed to capture fine scales. For instance, we have:
Residual-free bubble (RFB) methods. The finite element space is enriched
within each element by bubble functions that minimize the residual at the
element level [22, 20, 49, 87].
2

Two-level FEM. An approximate RFB method where the residual-free bub-
bles correspond to a finite element solution evaluated on a submesh at
each element [40, 48, 41].
Other methods attempt to capture global characteristics of the solution (the
coarse scales [58]), filtering the pollution of under-resolved fine scales, such as
Stabilized methods. The Galerkin variational formulation is modified with
mesh-dependent residual terms to guarantee stability in the sense that the
augmented bilinear forms satisfy a coercivity condition [19, 23, 33, 60].
Discontinuous Galerkin (DG) methods. This class of methods relaxes the
continuity of the approximate solution across element boundaries, impos-
ing a weak continuity of numerical fluxes [14, 29, 69, 77].
This thesis aims to extend the discontinuous enrichment methods to prob-
lems in CFD. The discontinuous enrichment method (DEM) is a hybrid finite
element method that has been applied to problems in acoustics [35]. The finite
element spaces are enriched within each element with fundamental solutions of
the differential equation. These enrichment functions are able to represent the
fine scales of the problem, improving the solution of standard polynomial fields.
In particular, we address the following problems:
Stokes equations.
3

The thesis makes new contributions in all these models and carefully con-
structs meaningful enrichment functions for each situation.
Research in discontinuous finite element methods has increased after the
development of new DG methods for second-order problems [29, 77]. One of the
open questions in this area is how to derive numerical fluxes [5] in a systematic
manner. Unlike DG or stabilized methods, DEM seeks stability and accuracy
using knowledge about the differential operator and does not depend on ad-hoc
perturbation terms defined on the interior or at the boundary of each element.
The outline of this thesis is as follows:
The next chapter reviews hybrid finite element methods and uses the
Dirichlet problem for the Laplace operator to introduce the DEM. Current
discontinuous finite elements [77, 14] are reviewed as well.
Linear, first-order problems are considered in Chapter 3. This chapter
presents a link between the standard discontinuous Galerkin method [69]
and the mortar element method [14], which together with the RFB method
[20] provide the basic ingredients for the development of a discontinuous
enrichment method.
The theoretical and experimental results for the hyperbolic problem mo-
tivate an extension of DEM to the advection-diffusion in Chapter 4. The
method is enhanced with exponential boundary-layer functions formerly
derived by the method of asymptotic expansions [62].
The objective of Chapter 5 is to stabilize the bilinear-velocity-pressure
4

pair for the Stokes problem. Homogeneous, harmonic polynomials provide
such stabilization without resorting to perturbation terms.
The chapters above also review the literature and include numerical error
estimates. A summary of the main results is presented in Chapter 6.
5

2. Hybrid Methods for the Poisson Problem
This chapter has three goals:
Review finite element notation and concepts employed throughout this
thesis;
Introduce the foundations of the discontinuous enrichment method; and
Formulate a discontinuous enrichment method for the Poisson problem,
one of the most common elliptic problems in CFD.
The material presented herein is based on established textbooks in finite
elements [18, 27, 57, 64] and the articles of Farhat et al [35] on the discontinuous
enrichment method, and Raviart and Thomas [81] on the primal hybrid method.
Some aspects of the mathematical theory of finite element methods are
needed. They are briefly described on sections 2.1-2.4. In particular, the primal
hybrid formulation (Sec. 2.2) is further explored. The discrete approximations
of the hybrid formulation, that is, the primal hybrid methods (Sec. 2.5) admit
several variants, such as the mortar element method (Sec. 2.7) and eventually
the discontinuous enrichment method (Sec. 2.8). The hybrid formulation can
be further related to non-hybrid discontinuous methods (Sec. 2.6).
2.1 Functional Settings
Let us introduce some standard spaces of functions and respective inner
products and norms. These functions are defined in an open, bounded set e
6

Mn (or an open subset of Cl) with Lipschitz-continuous boundary [27]. We
represent the boundary of Cl by T. This thesis is particularly concerned with
two-dimensional problems (n = 2) and real-valued functions.
We denote the space of measurable, square integrable functions in Q by
L2(Ci). The space L2(Cl) is a Hilbert space with respect to the inner product
Mo ,a = /
J n
uv dCl
(2.1)
and a Banach space with respect to the norm Iloilo n = (u, ufJn- Let us also
consider the following subspaces of L2{Cl), known as Sobolev spaces:
Hm(Cl) = (ue L2(Cl) ; M e L2(Cl) 0 < i +j < m } ,
dxidyi
equipped with

(M V^m,n ? \ d&dyi d&dyi / 0>n
h3
Another subspace of interest is the space of functions that are continuous
over fl, denoted by C(fl). We have that (70(fl) is a Banach space with respect
to the norm
IMIoo,n = sup|u()| .
If a function v is differentiable in Cl and its partial derivatives dv/dx and
dv/dy are continuous in Cl, we say that v is of class C1. The space of all functions
of class C1 is denoted by In general, the space of functions whose partial
derivatives of order < k are continuous (i.e., functions of class Ck) is denoted
by
7

Remark 2.1 : The subscripts are dropped unless the domain is a subset of
or m > 0. For example, IHI^- = ||n||0 K and Wu^ = little n.
All definitions above may be extended to vector-valued functions. For in-
stance, the inner product of (L2(Q))2 is defined as follows:
(u,v) = I u v dQ,.
Jn
Since Vu G (T2(f2))2 for any u G H1^), (u,v)1 can be written as
(u, v)1 = (it, v) + (Vit, Vv) .
One can also introduce a semi-norm in Hx(fi):
We frequently employ the following subspace of (T2(fi))2:
H(div; ft) = {q G L2(fi)2 ; V q G L2(ft)} ,
with the norm
I9IU =(IIII2 + I|V'I|2)1/2
Remark 2.2 : Scalars are written in italics, while arrays are written in
bold italics, unless they are represented by their components. For example, a
point x G ]R2 and a matrix A can be respectively written as (x, y) and A^.
Remark 2.3 : |a| represents the Euclidean norm of the vector o; that is,
\a\ =
+ a
2
2
8

The space of measurable, square integrable functions on T is denoted by
L2(r) and associated with the following inner product:
(u>v) = f
uv dT .
Let tr be the trace operator from H1^) to Z>2(r), denoted by tr v = v |r-
We define
H',2(r) = fjel2(r); |r=9 De/f'ffi)} ,
llslli/2 = iff IMIi I V1 = {veff,(fl) ; t)|r=s} >
' VeVg
H](Q.) = {v e H^fl) ; v |r= 0} .
Remark 2.4 : The inner product (-, ) may be used to represent any line
integrals over T. A subscript is employed when the integration domain is other
than T.
We denote the dual space of Hl/2(Y) by H 1//2(F) and define its norm as
lk'll-i/2 = SUP
{o, 9)
sup
(
(2.2)
|| || || ||
geHi/2(r) II0II1/2 veH^n) Fill
Remark 2.5 : The supremum or the infimum of a fraction excludes the cases
where the denominator is zero.
Let n be the unit vector that is outward normal to the boundary of (a subset
of) Q. Greens formula [18, Lemma III.1.1] relates H{div;Q) with R-1/2(r) :
Lemma 2.1 : If q e J?(div; Q) then q n |re r) and
(V q, v) + (q, Vv) = (q n, v) V v G H1^) . m
9

2.2 The Primal Hybrid Formulation
Given a function / : Yl > 1R, let u : > M be the solution of
A u = f in ,
(2.3)
u 0 on r .
Problem (2.3) is known as the Dirichlet problem for the Laplace operator, or
the Poisson problem. The existence, uniqueness and regularity of the solution
u have been studied in several contexts [71]. In particular, it is well known
[27] that if / G L2{fl) and Y is Lipschitz continuous, then there exists a unique
generalized solution u E V = Hq(Q) of problem (2.3) in the sense that it is also
a solution to the variational problem
u G V ; a(u,v) = F(v) V v G V , (2.4)
a(u,v) = (Vu,Vv) F(v) = (f,v).
Since a(-, ) is a symmetric form, we have that (2.4) is equivalent to the
following minimization problem [27]:
J(u) = min(J(v) ; v V} J(v) = ^a(u,u) F(v) . (2.5)
Z
Consider the abstract minimization problem
f(x*) = min{/(a;) ; x G A fl B} . (2.6)
We can also write (2.6) as a constrained minimization problem:
f(x*) = min{/(a:) ; x A}
subject to x B .
10

One can use the same argument to replace V in (2.5) by a larger space X of
admissible functions. In particular, X may be related to the domain partition.
We define a partition Th of O, to be a subdivision of into nel elements ile such
that
nel
1. n = |J
e=l
2. If i ^ j, IT n is either an edge, a vertex or the empty set.
A function v e L2(fi) belongs to A if v |oeÂ£ H1^) for 1 < e < nel, i.e.,
nel
x=n^1^6) (2-7)
e=l
Given a fixed order of the elements Qe, we define the jump {v} of v G X as
[I(a5) = v()|n-t)(aj)|n# ccGfTnfV i>j. (2.8)
We also define
nel nel
= , f-|Jre , = .
e=l e=l
Notice that functions v G X may be discontinuous across interelement
boundaries. Thus, we must consider extensions of the operators a(-, ) and /():
nel ^
a{u,v) = J2(Vu,Xv)ae = (Xu,Vv)fl J{v) = -a(v, v) F(v) .
e=l
The solution of (2.5) also satisfies the constrained minimization problem [35]
J(u) = min{J(u) ; v G X}
subject to
H = o on f\r ,
v = 0 on r .
(2.9)
11

Problem (2.9) is then transformed into an unconstrained minimization problem
with the aid of Lagrange multipliers. Let us introduce
f nel m = JJfr1/2(re) 1 z>i ; Bq G H (div; fi) ; q n = g on f > . , (2.10)
V ei / nel
&(>m) = - ^2 (v, m)r = (v> f e=l (2.11)
We have that (2.9) is equivalent to the following:
C(u, A) = minmaxÂ£(u,//) C(v, g) = J(v) + b(v, g) . (2.12)
vX
The saddle-point equations for Â£(v,g) yield the hybrid formulation [27]
(, A) e X x M ;
a(u, v) + b(v, A) = F(v) V v E X ,
b(u, g) =0 V g M .
If non-homogeneous boundary conditions are imposed on u, i.e.,
-Au = / in Q ,
u g on T ,
then the corresponding primal hybrid formulation is
(, A) G X x M ;
a(u, v) + b(v, A) = F(v) V v E X
b(u, ii)
= -( (2.13)
(2.14)
(2.15)
Formulations (2.13) and (2.15) fit into the framework of mixed methods,
which we discuss in the following section.
12

2.3 Mixed and Hybrid Methods
We review the abstract theory of mixed and hybrid methods presented in
[18, Chap. 2]. Let V, Q be Hilbert spaces and / G V\ g G Q'. Consider the
continuous bilinear forms a : V x V > M and b : V x Q M with associated
linear, bounded operators A : V > V' and B :V Q' such that
a(u,v) = (Au,v)V'xV ,
b(u,q) = (Bu,q)Q'xQ = (u,Btq)VxV>
We consider the following variational problem:
(2.16)
o(,) + p) = (/, v)v>xv
Ku,q) = (9,q)q'xq
VveV ,
V q G Q .
(2.17)
Let IMlQYKerjB* = inf{||g + 9II, 9 Â£ KerB1}. The following theorem [18,
Th II.l.l] describes the existence and uniqueness of solutions of (2.17).
Theorem 2.1 : If there exist positive constants a, (3 such that
( r a(u,v)
mi sup Tj rj > a ,
uEKerB j,gKer.B ll^llv ll^llv
< (2.18)
. , a(u,v)
mi sup > a ,
vEKer^ B ||^||y ll^lly
SUP ~ & H9H(3\KerB , (2.19)
v<=V |P||y
then there exists a unique solution (u,p) eVx Q \ KevB1 to problem (2.17) for
any / G V' and g G ImB.
Let 14 and Qh be finite-dimensional spaces of V and Q, respectively. A
Galerkin approximation of (2.17) is given by the solution (u^Ph) of the discrete
13

problem
/
a(uh,vh) + b{vh,ph) = {f,vh)V'xV V vh G Vh ,
<
b(uh, Qh) = (9, qh)Q'xQ VqhtQh-
V
We define Ah : Vh ^ Vf[ and Bh : Vh - as in (2.16).
following [18, Th II.2.1]:
(2.20)
We have the
Theorem 2.2 : If there exists Vh 14 such that b(vh, qh) = (g, qh)q'xQ for
any G Qh and there exist positive constants o;o and k0 independent of h such
that
a(uh,Vh)
Uglily H^/illy
inf
> a0
(2.21)
inf sup
b(vh, qh)
> k0
qhtQh vhÂ£Vh ll^/illy 11 11 Q\KerjS^
then (2.20) has a unique solution (Uh,Ph) Â£ 14 x Qh \ KerBlh and
ll MW + \\P ~ Pfc||g\Kerfl < C(/|f Vh\\V + IIP ~ ftllg )
(2.22)
Condition (2.22) is known as inf-sup or Babuska-Brezzi condition. One may
consider inequality (2.22) with ||-||g instead of the quotient norm |Mlg\KerS* when
spurious modes qh G KerBl [18] are difiicult to eliminate. In this situation, the
constant k0 may depend on h and the following result [18, Prop. II.2.1] applies:
Theorem 2.3 : If there exists t)/, 6 14 such that b(vh,qh) = ( any G Qh and there exist positive constants ah and kh such that
a(uh, vh)
inf sup
U/jGKerfl/j VflÂ£KerBh 11 ^h 11 y 11 w/i 11 y
> h
(2.23)
14

(2.24)
inf sup
b{vh,qh)
> kh
Qh.^QhVhevh H^/illv II^aIIq
then (2.20) has a unique solution (Uh,Ph) Â£ Vh x Qh\ KerB^.
Remark 2.6 : Condition (2.18), which corresponds to the invertibility of a(*, )
in KerB, is fulfilled if a(-, ) is coercive in KerR; that is,
a(u, u) > a\\u\\y VuG KerR ,
and a similar argument holds for (2.23).
Let us analyze (2.13) according to Theorem 2.1. Set V = X, Q = M and
nel
Mlv = ElMli;
e=l
Because
sup
,n >
b(v, fi)
nel
<2 = E
e=l
Q >
-i/2,re
ueV |P||y
Ker_Bs = {0} and (2.19) follows. The following lemma [81] shows that Ker.6 =
nan).
Lemma 2.2 : A continuous linear functional L : X > JR vanishes on Hq(Q)
if and only if there exists a unique n 6 M such that
L(v) = b(v, fi) V v G X .
Therefore the existence and uniqueness of the solution of (2.13) follow from
Theorem 2.1. Raviart and Thomas [81] proved this result and presented an
15

Theorem 2.4 : Problem (2.13) has a unique solution (u, A) G X x M. More-
over, u G Hq(Q,), u is a solution of (2.3) and
, du
A = 7 = vu n on r e = 1,..., nel.
on
Remark 2.7 : The standard norm for the space M [18, 27, 81] is given by
IHU = inf{llg|ldiv ;getf(div;0) q-n = p} ,
which is equivalent to ||-||q [81, Remark 4].
We now consider finite element approximations of (2.13), which are known as
primal hybrid methods. Given finite-dimensional spaces Xh C X and Mh C M,
we want to find (u^ Xh) G Xh x Mh such that
{d{uh, vh) + b(vh, Xh) = F(vh) V vh G Xh , .
(2.25)
6(ttfc,A*fc) =0 VftGl),.
Notice that
KerRft = {vh ; &(wA, p,h) = 0 V/i/,G Afft} ,
KerS^ = {/*fc G Mh ; 6(vft, p,h) = 0 V G Xfc} .
One can eliminate the Lagrange multiplier by restricting Xh to KerR^. In
this case, (2.25) reduces to finding Uh G KerBh such that
a(uh, vh) = F(vh) V^G KerBh . (2.26)
The spaces Xh and Mh must be designed such that conditions (2.23)-(2.24)
hold. Similar conditions are given as follows [81, Th. 2]:
16

Theorem 2.5 : If a(-, ) defines a norm in KerÂ£>/t; that is,
{vh G KerBh ; a(vh,vh) = 0} = {0} , (2.27)
then problem (2.26) has a unique solution Uh G Ker Bh- Furthermore, if Ker Blh =
{0}, then (2.25) has a unique solution (Uh, A/J G Xh x M^. m
Before introducing the classical primal hybrid approximation, we describe
the domain partition in the following section.
2.4 The Reference Element
We focus on convex polygonal domains Cl and partitions into quadrilateral
elements. The boundary of an element Cle is denoted by Te, while the i-th side
and the i-th vertex of Te, i = 1,..., 4 are respectively denoted by 7 and xf, as
shown in (Fig. 2.4). An arbitrary side of T is denoted by 7.
Figure 2.1: A representative element in a partition of Q
Bilinear interpolation is associated with quadrilateral elements. We define
the space of bilinear polynomials over an element Oe as follows:
Q1 (Cle) = {p : Cle -> JR ; p(x, y) = o0 + axx + a2y + a3xy a0, , a3 G M} .
17

Most of the finite element computations are performed on the reference
element, which can be formally defined [27] as {K, Qi(K),^}, where
K = [-1,1] x [-1,1] ,
X = ,
*1 = (-1,-1) 2 = (1,-1) *3 = (1,1) *4 = (-1,1)
y x%

i
*1 x\
Figure 2.2: The reference element
The basis for Qi(K) is formed by the polynomials Nu i = 1,..., 4, satisfying
Ni(xj) = Sij] that is,
ft) = (1 + ?)4(1 + n) , m, n) = (1 ;)4(1 + ,,)
Let Cle G %, be an arbitrary element. We denote by F the isoparametric
transformation which maps each point (Â£, rf) G K to (x, y) G Qe, where
x(^,T]) = xe1Nl(f,Ti) +...+xtNi(^,r]) ,
y(Â£,v) = yfNi{,y) + --- + (Â£,r?).
(2.28)
18

The Jacobian of F is given by
J = J(f,v) =
xAtV) x,v(,ri)
y,d^v) Vxi&v)
(2.29)
We assume that f)e is such that J is invertible, unless otherwise stated.
Remark 2.8 : We frequently employ the indicial notation for derivatives, as
seen in (2.29).
Example 2.1 : Let Â£le be a square element of length h; that is, fle = [x0, Xq +
h\ x [yo,yo + h\ for some (x0,y0) G JR2. Substituting x\, yf, i 1,..., 4 into
(2.28) yields
Although the map F is generally not affine, it is affine at the element bound-
Further, J = Ih/2 and \ J\ = h2/4.
Given u : Â£le JR, let
&(Â£> V) = u{F(Â£, rj)) = (a;($, /?), y{Â£, rj)) - We can also write u = u o F 1, as seen in the integration formula ary f = OK. We define Sm to be the space of piecewise polynomials of degree m over T. 19 Proposition 2.1 : The restriction F\f. f > Te is a piecewise affine mapping. Furthermore, given v G L2(Fe), there exists fro G S0 such that [ vdT= l Jre Jr vfro dV . Proof : Consider the restriction of F to the side \r\ = 1. From (2.28), *(0=*(Â£,-l) = ^^ + Â£^-^ (2.31) ,r\ (r i\ 2/2 + 2/1 rV2-yi 2/(0 = 2/(0l) = 2 + ^2 which shows that F |-y, is affine. Consider the parametrization of 7 given by (2.31). We have ds y/dx2 + dy2 ^y/{x2- xfr)2 + (y2 yiY df, =$(%.
The proof follows by proceeding as above with 7*,
> 2 .
Remark 2.9 : Consider the space
Qi(fie) = {v G H\Fle) v = voF~1 v G Qi(K)} .
Although Qi(Qe) = Qi(Qe) only if F is an affine mapping, we also denote
Qi(f2e) by <3i(^e)-
2.4.1 Partial Derivatives in the Reference Element
We wish to express partial derivatives of u G C2(Fle) with respect to partial
derivatives of u. The independent variables (x, y) and (Â£, rj) are also referred to
20

as (:ri,Â£2) and (Â£i,Â£2), respectively. The indices i and j (k and /) are associated
with the variables (x,y) ( (Â£, rf) ). The summation symbols are omitted.
From the chain rule, the first-order derivatives of u can be written as
du d^k du
(2.32)
dxj dxj dÂ£k
Let A be defined as Akj = d^k/dx3. We can write (2.32) in vector notation
Vw = A^u = J Vti.
(2.33)
(2.34)
Second-order derivatives are computed with the chain and the product rules:
d2u d (dÂ£k du
dxidxj dxi l dxj dÂ£k
_ du
dxi ydxj dÂ£k
= & (dÂ£k\du d& dÂ£k d2u
dxi d& ^ dxj J dÂ£k dxi dx3 dÂ£kd&
= 9 fdÂ£k\du d& d f dÂ£k\du d& dÂ£k d2u
dxi ydxj J dÂ£k dxi d& J dÂ£k dxt dx3 dÂ£kd& '
Identity (2.34) is written in indicial notation as
'a,ij A\i^kjiu^k + A2i^kj2^,k T utki.AhAkj .
Let us denote the Hessian of a function / by H(f); that is,
H(/) = f,xx f,xy
f,xy f,yy
Let B and C be respectively defined as Bkj = Â£kji and Ckj = ^kj2] that is,
dJ-1 dJ-1
B , C x .
dr]
21

Therefore,
^,ij WJtlAnAkj + AuBkjU^k + A^iCkj'Qjk (2.35)
Let J-1^ be the ft-th row of J L We write (2.35) in vector notation as
H(u) = J-'-iJ(w)-J~l+J-1(1)
+J~1(2)( ^-Vti ) (2.36)
2.4.2 The Laplacian
The Laplacian of a function u can be written as the trace of its Hessian:
'U'j.j u,,ki-AiiAki T AiiBikutk T , or
'dJ-* ^
(2.37)
Au = (J"1 J"4) : H{u) + Vu)* J4(1) + -Vm)'- J~4(2) .

The inverse of J can be written as
J l = A7
drj
V,n x,n
-VA xA
Thus,
(J 1J 4) : H(u) = 2 [(^2 + y,5), + (zj + J/J),-
m
2(x,vxa + y
Since z and y are bilinear functions of Â£ and rj,
dJ
-t
d( \jy
-t
dJ
dr\
\jf
VÂ£v 1 Jl y, v l*d>C i/,d Jl.f
~{xAn \J\ ~ x,n l*d>Â£) -z,dJk
~y,v I'd? ~{yAn l*d y, ? l*d?)
X,V I'd? xAn I'd XA l*d?
(2.38)
(2.39)
(2.40)
(2.41)
22

The dot product of (2.40)-(2.41) and Vu yields
/
dJ~* 1
aT'v = T
(vavlJ\ v,v lJU)* + vaIJUfl,?
\
(2.42)
dJ~l 1
----V =
9rj \J\
\ (XAn 1*^1 l^1>Â£ )^A j
~V,v I*^Ij7 1*^1 y,(, \J\in )^,fi ^
(2.43)
^ x,n \J\rri + (2:,^ |7| x^\J\,v)u^ y
Let us multiply (2.42)-(2.43) by the columns of J l\
dJ
-t

Vfi
UA
y,v y,tv v,n
\Jk
+ X,T) I #,Â£?? X;7)
\J\*
u
,v
\J
|J|,Â£ |J|*
V^~\J\ + X*X*~\J\
x,nxAv + V,ny,tn ~ T7T (xl + yl)
7jk(xAx,v + y,ty,v)Tff
(2.44)
dJ
-t
dr)
Vu\ J~t{2)
UA
\J\2
^.n
*^577 1 *^57?
VMfi-rFi + x>(x*
\J\
+
\J\,v \ ( \JU
y,z 1 vtn y,s-rf\ + x& ~x^
T^2(xAx,r, + yAy,v)-
\iv
+
U.
xAxAn + VAVAn ~ ^ (^| + 1/j)
(2.45)
23

Substituting (2.39) and (2.44)-(2.45) into (2.37) yields
A u
(*, S + - 2(x,vx,t + V*V,n)uAri +
x,vx,tv + yytv + -jj^(x,$x,v + y*y,v) Tjf(x* + y) K? + x*xÂ£v + VtVtTt + Tjf(x,^x,v + y&n) - + y$ ) ^
. (2.46)
A .
The results above are applicable when u = Nf = N o F i = 1,..., 4.
From (2.33) and (2.38) we have
/
VNf
\
(2.47)
y,r]NiÂ£ y,i^i,V
^ ~x,nNiÂ£ + x^NitV ^
Since each JV* is bilinear, we have = Nim = 0 and (2.46) reduces to
1
A Nf =
~ 2(v,? + y^y,n)Nify +
\J\
(2.48)
x,vx,tv + y,r,y&i + T/f + y*y*) njfW + v*)) +
xÂ£x*v + v&tn + T7f(xÂ£x,n + y^n) ~ -^n t"2 1 "2
W + vJ) ^
2.5 Primal Hybrid Method
We now present the finite element spaces employed in the classical primal
hybrid method in the particular case of bilinear elements and piecewise constant
multipliers. For triangular and/or higher order elements we refer to [81]. Let S0
as in Sec. 2.4. We define
24

SI = , Â£e&} ,
(2.49)
nel
Mh = jyiifc e fj S; fih |r +/jth |rj= 0 where PnF^ 0 j , (2.50) nel Xh = l[Qe Qe = {v eH\ne) ; v^voF-1 v e Q} , (2.51) e=l where Q is a subspace of Z/1 (KJ to be defined later. We can also write as [35] Mh = j/* 6 nH~mT) ; 39a BDFAf, ; BDFMi = |qh E H(div;ft) ; qh |n#= (ci + c2a;,C3 + c4y)f , ci,... ,04 E JR and 1 < e < nel j . Lemma 2.3 : a(-, ) defines a norm in KerS^ Proof : a{vh,Vh) = 0 implies that Vh |n= is constant for any e. Given an interior edge 7 = PnP, choose fiu e Mh such that 1, e = i , fJ>h |ne= 1, e = j , 0, otherwise . s. If vh G KerBh, in particular b(vh,fj,k) = vh |ni -vh |m= 0. Since 7 G T \ T was arbitrary, Vh = no for some constant vq. Choosing 7 6 T and repeating the argument above yields v(} = 0. The condition KerB^ = {0} can be replaced by a stronger one which is easier to check ([81, Lemma 2]): 25 (2.52) Lemma 2.4 : KerBlh = {0} if Ker Bl = j /} e So ; J vftdT = 0 V v G Q ^ = {0} . Proof : Let ///t G KerBlh. In particular, [ (vfxh) oF~1dr = 0 VveQ. J re From Prop. 2.1, there exists /*o S0 such that J vjihfio dT = 0 V v G Q Since /i/j/io G S0, fthfro = 0, which implies fih = 0. Thus, /i* |re= 0. Since CP was arbitrary, /j,h = 0. To select a space Q G H1(K), we must first notice that Qi(K) x Mh may not satisfy (2.24). Indeed, if Q = Qi(K), then Â£*(*) = 1, x G 7i ri73 , -1, JC G 72 n 74 . would satisfy 0 ^ ft* G Ker iK In order to assure (2.52), we define n2 Â£2 Q = Qi(K) U span{ur} vrt{Â£,v) =--------; (2.53) (2.54) Lemma 2.5 : The pair Xh x Mh x Q given by (2.54) satisfies condition (2.24). Proof: Let p,h G Ker BL Since <2i(K) C Q, we have in particular J v/Xfi dt = 0 V v G Q . (2.55) 26 Condition (2.55) implies /ih cÂ£l* for some c M [81, Lemma 7]. Then, by choosing v = vrt yields a , which implies c 0; that is, Ker B = 0. The proof follows from Lemma 2.4. From Lemmas (2.3)-(2.4) and Theorem 2.5, The hybrid approximate solu- tion defined by Xh x Mh exists and is unique. The inf-sup condition was proven [81, Lemma 10] for mesh-dependent norms under the assumptions that the mesh is regular and the mappings from K to Â£le are affine for any e. 2.6 Other Discontinuous Formulations Hybrid methods are often criticized because of the additional degrees of freedom associated with the Lagrange multiplier. This motivated the design of formulations where the multiplier unknowns are eliminated. The classical elimination procedure is to restrict Xh to KerE^, which yields (2.26) Since in general KerBh KerH = H\(Q) (see Sec. 3.2.1), we have that (2.26) is a non-conforming method. In fact, several non-conforming methods can be written as conforming, primal hybrid methods [81]. Another approach is to approximate the saddle-point problem (2.12): C(u) = minÂ£(u) = minÂ£(u, ft(v)) Â£(w, A) . vÂ£X v&X 27 If fj, depends linearly on v, the optimality condition C(u + ev) C(u) lim 1------------ = 0 V v e X . e0 yields the following variational formulation: Find u 6 X such that nel a(u, v) ^2 (u> Hv))r + (v-> Ku))r^ ~ F(v) V v G X . (2.56) e=l Theorem 2.4 suggests choosing p(v) = Vu n |pc, e = 1,... ,nel. However, p,(v) g M in general. Instead, we consider the average flux {Vu} n, 7nri ) , 7 = Fl n Tj i j , , 7 = P n r. {Vu}|7nr<= < K' Vv 7nri +Vu (2.57) Vo |-(nr' V Method (2.56) with pt(v) given by (2.57) corresponds to the global element method (GEM) of Delves and Hall [31, Th. 4], This relationship was mentioned by Oden et al [77], who proposed the following modification of GEM: nel a(u, v) + ^ (it, {Vu} n)re (v, {V} n)re = F(v) Vuel. (2.58) e=l Arnold et al [5] proposed a common framework which includes (2.58), the local discontinuous Galerkin [29, 90] and interior penalty methods [76, 6]. 2.7 Mortar Element Method The mortar element method [14] is a non-conforming domain decomposi- tion method that weakly imposes continuity between subdomain interfaces by introducing a finite-dimensional space of functions (the mortar space) defined on these interfaces. We consider a primal hybrid version [12] where the mortar space is the subspace of Lagrange multipliers. 28 Consider a partition of f1 into non-overlapping subdomains Dk, 1 < k < Np\ that is, ~Nd n = (J Dk. k=1 Let Xh be a finite-dimensional subspace of nd X = ]jH\Dk) . k=1 Two adjacent subdomains Dk and Dl may induce different partitions of their interface ^{k,i} = dDk fl dDl (Fig. 2.3). Let us denote the partition of T induced by Dk and Dl as and T/*, respectively. The mortar space presumes a unique classification of these partitions as mortars or nonmortars. The finer partition (r^ in Fig. 2.3) is usually selected as the mortar [14, Remark 3.2], r'ik l Dk 7 Dl 7 r ki Figure 2.3: Mortars and nonmortars of a subdomain interface Let Mf be a one-dimensional piecewise-polynomial continuous space defined over the nonmortar partition. The continuity of a function Vh Xh on T{^,1} is 29 imposed by the following orthogonality constraint: (vh\Dk -vh \DhVh)r{k l} =0 V nh G Mj? . In general, assume that the union of all interfaces is decomposed as nd s=UU r - k=1 l where the sets JA{k) C {/ G {1,,..., No} ; / 0} satisfy either l G Al(fc) or k G M(l) whenever T^i} 7^ 0- We define Mh = K (2-59) rfcJes Consider the following bilinear forms: nd aD(u,v) = ]T(Vu,V?V , (2.60) k=1 bD(v,fj) = (M>/*)rw (2-61) r mÂ£S The mortar element method for (2.3) becomes finding (u^, u/t) G 14 x Mh such that aD(uh,vh) = (f,vh) Vu G 14 = {v G Xh ; 6c(u,^) = 0 V [i G M/J . (2.62) Let us denote the endpoints of a partition Fkl as xh and xbkl and an arbitrary element (segment) of as 7. Note that each element 7 corresponds to an edge of some element fle G Dk (Fig 2.3). Ben Belgacem [12] considered the mortar element method (2.62) on trian- gular meshes, proving stability of the pair Xh x Mh given by (2.59) and 30 nd x* = n*t > k=1 xkh = {n G <7(Â£>fc) ; t; |n-e P,(f2e) V Qe G D*, n U*nr= 0} , Mh = {v> G Wfc' ; ^l7e A-i(t) if e 7 or G 7> , W? = 0* e C(f{H}); /x|7G Pg(7) V7 g r} , where Pq{A) represents the space of polynomials p : A > JR of degree < q. 2.8 Discontinuous Enrichment Method The primal hybrid method enriches the space Qiifi6) within each element with the function vRT o F~l to guarantee that Xh and Qh satisfy (2.22). One could actually enrich Qi(Qe) in order to reduce the residual Res(v) = / (An) within each element, according to the RFB philosophy. For instance, consider the space Qeres = span{Res(n), n G Qi(Oe)}; we could enrich Qi(^e) with Xres = span{n G H1^6) \ Qi(fie) ; -An = w w G Qeres} . In particular, we have the following subspaces of X^es: X0e = {n G H^{Qe) ; -An = f + Aw, w E Qi(tte)} , XeH = {n G \ Q^Q?) ; -An = 0} . The former relates primal hybrid methods to RFB methods (although it may lead to a pair Xh x Mh such that KerBLh ^ {0}), The latter allows us to capture the homogeneous component of the solution. A discontinuous enrichment method for (2.3) can be defined by considering the primal hybrid formulation (2.13) with the spaces Mh and Xh defined as in 31 (2.50)-(2.51), with Qe defined as Qe = Q1(Qe) UV| , where VÂ§ is a finite-dimensional subspace of (2.63) Example 2.2 : Assume that Th is composed of square elements of length h. From (2.30) and (2.46), ve = vrt o F~l satisfies Ave 4}i~2Avrt = 0 , that is, ve e Xfj C X^es. Therefore we can set VJ = span{ue}, which recovers the classical primal hybrid method. Let us consider &(,) as hi (2.11). Given a general differential operator C associated with a bilinear form a(-, ), we define the DEM for the problem 1Cu = f in , u = g on T . as finding (Â£, uf, Ah) Vf x V x Wh = 14 x 144 such that a(uh, vh) + vf) + Kv*i= (/<) V vÂ£ e Vf , < > vh) + ) + Kvhi xh) = (/, ) V uf e Vf , k &(& > /%) + KuhMfc) = (^fc, The space 144 is a finite-dimensional subspace of M, while nel Vh = IlVr Vp = {ve H\ne) ; v = vo F~l veVP} e=l nel vf = m e=l (2.64) (2.65) (2.66) (2.67) 32 where Vp is some polynomial interpolation space, V% is a finite-dimensional space of Vrees, Ve r res spanj/u G i?1(f2e) \ Vfi ; Cu v . . VtQles) oe res = span{u G H1 (Q,e) ; v = f Cw , weVfi} . and C, / are local approximations of C and /, respectively. For instance, C could be an asymptotic approximation of C, while / could be the element average of /. Sometimes C = C and/or f = f as well. Remark 2.10 : Formulation (2.65) may be adapted to allow Neumann and/or Robin boundary conditions [35, Sec. 3.3]. The approximation space of the primal variable is split as 14 = Ifif U V)f. We refer to spaces V* and V as the polynomial field and the enrichment field, respectively. Farhat et al [35] also associated V/p and V)f with the coarse and fine scales, according to the variational multiscale terminology [58]. The original DEM formulation [35] is distinct from (2.65) in two aspects: the polynomial field is piecewise continuous; the enrichment field is a subspace of Vfj = span{w E Vfies ; Cu = 0} . The functions v E Vfi are referred to as free-space solutions of the operator C. They satisfy the homogeneous differential equation Cu = 0 but are not tied to boundary conditions or to the geometry of the element [35]. 33 The proposed extension of the enrichment field allows DEM not only to capture the fine scales associated with the underlying homogeneous problem, but also to decrease the residual of the polynomial field at the element level. Since the enrichments can be statically condensed as described in [35], the amount of computer memory does not change from the original approach. Relaxing the continuity of the polynomial field increases the computational effort, but seems more appropriate for (nearly) hyperbolic problems and allows the use of similar arguments from the theoretical analysis of the primal hybrid method (for instance, Lemma 2.4). 34 3. Linear Advection Equation Let a = (ai(x,y),a2(x,y)Y be a vector field in Q. The following boundary value problem is a linear model for the non-diffusive transport of a contaminant with concentration u = u(x, y) in a fluid which moves with velocity a. The boundary of a subset K C O is classified according to a as follows: We refer to 3AT, dK+ and dK0 respectively as the inflow, outflow and tangential boundaries of fi. They represent the portions of dK where the fluid enters (3AT), leaves (dK+), or flows along the boundary (dK0). The concentration is modelled by the solution of a first-order, linear partial differential equation subject to a force / and an inflow boundary condition u g\ 3AT {x G dK ; a n < 0} , dK+ = {a? e dK ; a n > 0} , dKo = {x 6 dK ; a n = 0} . dK0 dK+ dK_ Figure 3.1: Inflow, outflow and tangential boundaries of K (3.1) 35 The right-hand side / represents a source/sink of the contaminant, whose concentration is prescribed as ^ on T_. The equation a Vu = / is known as the (linear) advection equation. If a, f are of class C1 and T_ can be parameterized by a C1 curve, then problem (3.1) has a unique solution in the neighborhood of T_ [96]. The solution may exhibit discontinuities across the characteristic curves, which are given by ai(x,y)x + a2(x,y)y = 0 . (3.2) Let 7 = (7i(t),72(t)) be a solution of (3.2) and let u(t) = u(7i(t),72(i)). The first equation of (3.1) can be reduced to an ordinary differential equation: ^ = /(7i(*),72(*)) (3.3) Problem (3.1) constitutes a valuable tool for validating finite element for- mulations, since explicit solutions may be available and the same numerical difficulties arising from discontinuities are present in more complex transport problems. 3.1 Finite Element Methods The finite element methods described below allow the test functions vh to be discontinuous across the element boundaries. Given an interior edge 7 G f, the trace v^ |7 may be undefined. However, one can define approximate traces based on the well-defined traces of restrictions of Vh to each element of % These approximate traces are known as numerical fluxes in the discontinuous-Galerkin literature. 36 For instance, let Qe be an arbitrary element and 7 Te be an interior edge. Let fle> be the element such that re< PI Te = 7. The trace VfL |7 may be defined either from the interior (vmt) or from the exterior (vext) of the element Q.e [15]: (3.4) Vint(x) = v(x) |ne , vext(x) = v(x) |ne/ x G 7 . We may use the velocity field to define globally the traces at the edges that are not aligned with a (Fig 3.2): vh(x) = tonvh(x + ta) , t^0 Vh (*) = ~ ta) (3.5) Figure 3.2: Traces of Vh from the inflow (vh) and outflow (vÂ£) Note that and can be locally expressed as v/>()|r* = xere+ , vÂ£(x) , x e n. If 7 is a boundary edge, we set vext(x) = 0. Analogously, v (x) = 0 and v+(x) = 0 if 7 e T_ and 7 e T+, respectively. Remark 3.1 : We omit the superscript of the traces evaluated from inside a fixed element (that is, the ones corresponding to the superscript int as in (3.4)) from here on. 37 3.1.1 Discontinuous Galerkin The discontinuous Galerkin (DG) method was proposed in [82] for the linear advection-reaction equation and analyzed in [67, 69]. The DG method for (3.1) can be found, for instance, in [64]. Let Qe be an interior element and u be the solution of (3.1). Multiply a Vu by a test function v that vanishes in D \ fte and integrate by parts: ua Vv dQ+ uva ndF + / uva ndF = fvdFl. (3.6) Jae Jy%. J rt We approximate u in De by Uh Â£ Qi(Â£le). The trace of is approximated as in (3.5); that is, Uha Vvh dVL + UhVhd n dF + uhVhd ndF = dfl . Another integration by parts yields (a Vufc, vh)a, (uh uh, vha n)re_ = (/, vh)ne . (3.7) If Tf = Ti fl r_ 7^ 0, we have (a Vuh, vh)ne (uh %, vha n)r, = (/, vh)^e (g, vha n)r. . (3.8) Let a(u, v) = (a Vu, v)q and nel ui==n,3'(ne)- <3-9) e=l The sum of equations (3.7) over all elements Fte yields the DG method: the goal is to find uh Â£ V* such that nel a(fc, Vh) ~^2(uh- , VfcO n)ri = (/, vh) {g, vha n)r_ V vh Â£ VÂ£ . (3.10) 38 Remark 3.2 : One may compute the DG solution for (3.1) element-by- element if the elements are ordered from inflow to outflow [69, 84]. Remark 3.3 : Note that the line integrals in (3.10) vanish at the edges aligned with a, so u^ is always well defined. 3.1.2 Stabilized DG Method (SDG) Some authors have added a streamline diffusion stabilization term to (3.10): nel nel /'ll l Dt. n. n. . (pe a(uh, vh) + ra (a Vu, a Vu)ne ^2(uh-uh, vha n)l e=l nel e=l (3.11) = (/>)+ 53 T (/>a' V)i>. i'J, l>a >r_ V ti* Vi P h e=l Houston et al [56] noticed that methods (3.10) and (3.11) have similar per- formance on quadrilateral elements and that the stabilization term rendered slightly more accurate results when high-order polynomials were employed. The stabilization parameter was originally proposed to be ra r = he, where he is the diameter of the element fle [65]. It was recognized later that ra should be proportional to he [33, 93], or in general to he/p, where p is the degree of the polynomial space [15, 56]. In particular, let a be the velocity evaluated the element centers. We consider the stabilization parameter proposed for the advection-diffusion equation [23, 39] in the advective limit for bilinear elements: he Tn 2 Idl where a is the velocity field evaluated at the element center. 39 3.1.3 The Mortar Element Method Let us consider a mortar element method where each element is a subdo- main, so that each edge is an interface. The velocity field provides a natural criterion for assigning the edges as mor- tars or non-mortars. Consider an interior edge 7 = r'+nr,_. The (non)mortar sides are associated with Cl1 (Â£T). Thus, the boundary of an element fle can be split into a mortar side and a nonmortar side Te_ (Fig. 3.3). mortars non-mortars Figure 3.3: Mortar and non-mortar sides of a 2x2 uniform mesh Let Xh = V* as in (3.9), a(u,v) = (a Vu, v)^, and nel S = U Fi . e=l We define a mortar element method for (3.1) by finding G 14 such that a(uh,vh) = (f,vh) M vheVh , (3.12) Vh {'U'h Â£ Xh j b(v,h, fJ>h) (9t Abi)r_ ^ f^h ^ Mh\ nel nel b{uh, fih) = Y] / [uh]nh dT = V / (fc uh)/Jth dr . e=i M e=i Jt- The choice of Mh is motivated by the DG method. Notice that the approx- imate solution of (3.12) would satisfy 40 if nel Wh) -^2(uh~ ul,ha n)ri = (/, wh) (g, wha n)T_ V wheVh e=l (3.13) / (ttA ul)vha ndÂ£ = 0 V vh e Qi(tte) e = 1,,... ,nel , (3.14) J n which is achieved by letting the space of multipliers Mh be nel Mh = T1 Mh > Mh = spanf^a n |ri, vh G Qi(fte)} . (3.15) e=4 Remark 3.4 : Method (3.12) with the multiplier space (3.15) resembles the method in [32], which incorporates a n into &(, ) (into the variational formu- lation rather than into the multiplier space). 3.2 Discontinuous Enrichment Method DEM was intended for solving second-order problems [35]. Since we are interested in advection-dominated problems, it is desirable that the method be able to solve first-order problems as well. In particular, problem (3.1) is the hyperbolic limit of the advection-diffusion equation considered in the next chapter. We propose a DEM for (3.1) motivated by the mortar element method (3.12) and by the RFB method for the advection-dominated advection-diffusion equation [20, 43]. The goal is to find (ttÂ£,itf, A/t) G 14 x Wh such that / > vh) + a(uh> vh) + KvL xh) = (/X) V vh e vh , a(uhWh) + a(uhWh)+b(vh^Xh) = (/X) V Kuh,nh) + Kufl,iih) V = (/bn 9) r_ V/^G Wh. 41 with a(u, v) = (a Vu, v)^ Q nel ^ (v, M)re\rj. ~ (w>/i)f\r+ e=l We let Wh = Mh as in (3.15) and 14 as in (2.66)-(2.67), with Vp = <2i(^e) and Vp = span{lV|}, where 1V| satisfies a V = 1 in Vte , u = 0 on n . Notice that &(, ) and Kercan be respectively written as nel 6( /*) = - > ^)ri > nel (3.17) e=l KerBj = < uh e 14 ; (tp, u,)ri = 0 V e 17 (3.18) e=l 3.2.1 The Subspace Wfe Choosing W4 as in (3.15) naturally raises the question of why we do not use piecewise constant multipliers as in the Poisson problem. Unfortunately, the space of constant multipliers does not form a stable pair with V^. Indeed, consider a mesh of a single element (the reference element K) and Vp = Qi(K) Vi = span{JV^} , Wi, = span{xi+Xl+J *,(?,>) = 1 > e((,/) = , 0 9(t,v) 50. In particular, v = Ni + + aNs + W4 (3Np satisfies Kv,Xi+v) = j *(*,-!) = f tdt-0 , &(t5,Xi+e) = f v(-l,fl)dT} = J r]dr] = 0 . 42 Prom the bilinearity of &(, ), v e KerF^. Further, we have that a Vv = -a 4 1 + a Srj 1 V + a ( \ y \ vJ /? 1 + a 3Â£ Choosing a; = 3 and fi = i + ei2 yields a Vv = 0. Therefore, a(u, v) is not invertible in KerBh. The same conclusion applies to a general domain Q, with quadrilateral ele- ments. Indeed, let Q.e C% such that Fj_ C T+ [69]. Then, v(x,y) v(F 1(a;, j/)) (x,y) e Qe , 0 (ar, y) & Qe satisfies a(v,v) 0, although v e Ker B \ {0}. On the other hand, renders a stable pair with V/t, as shown in Theorems 3.1-3.2 below. 3.2.2 The Subspace Let us resume the comparison between DG and the mortar element method (3.12). Unlike (3.10), additional stabilization is needed for method (3.12), as reported in [32] for the advection-reaction equation. Note that DG does not strictly satisfy restriction (3.14), but rather tries to balance it with the fulfillment of the differential equation, represented by the first integral of (3.7). Therefore, DG should perform better in the interior of the elements, while (3.12) should have an advantage on the interior edges and the inflow boundary, as shown in Fig. 3.4. We would like to enrich the polynomial space of the mortar element method so that a Vuh is closer to zero. We follow an approach based on [43], seeking 43 a function u% such that a W(uh + uf) = 0 in Qe , ul 0 on n (3.19) Assume for simplicity that the velocity field is constant and Ti = 7 U 74. In particular, if is the reference element K, we can write Uh\o?=ueh as U1N1 4- U2N2 + U3N3 + U4N4. . We have that a-Vueh + (Ui U2 + U3 114)0, VA3 (3.20) Notice that a Vueh = a V(C'i^i + 62^2), where a = bi in K , 4>i = 0 on f _ , bi = 1 62 = A3 Figure 3.4: Approximate solutions of Prob. 3.3.1 (Sec. 3.3) by DG (left) and mortar (right) methods 44 Thus, ul = C\4>i C22- Note that2 = -Ns is already spanned by the
polynomial field, so that only 1 needs to be added to the enrichment space.
For general elements, we add iVf as in (3.17). The following sections describe
how to calculate iVJ. and its partial derivatives.
Remark 3.5 : If Fe is a single side, then 4> 1 is also polynomial and no enrich-
ments need to be added at the element Q,e.
Remark 3.6 : The function 02 is not bilinear if Te is composed of three edges,
since a nonzero bilinear function vanishes at most on two edges.
3.2.3 Computation of IVJ,
Let us determine N% from (3.17). Consider Xp G Ye_ and the characteristic
curve 7 (t) defined by
x = Xp + a\t
y = yP + a2t.
(3.21)
Notice that w(t) = is the solution of
w'(t) = 1 , t > 0 ,
u;(0) = 0 ;
that is, w(t) = t. Solving (3.21) for t yields
(3.22)
NeE(x,y) =
x-xp y-yp
ai
0>2
(3.23)
for any {x, y) over the characteristic 7. We can also write (3.23) in vector form
[20] as
N% (x,y)
(3.24)
45

In general we need to determine xp from x. Consider two adjacent vertices
xA and xB of an element and the lines li and l2 as shown in Fig. 3.5.
Parameterize l\ and l2 as follows:
x' = x + ais x' = xA + (xB ~ xAt
h : , k : { (3.25)
y' = y + a2s [ y' = yB + (yB yA)t.
Figure 3.5: Location of xp
The intersection point (x\y') of li and l2 is given by the pair (s,/;) that
satisfies the linear system
(xB xA)t + (Qi)s x xA
{yB -yA)t + (-a2)s = y-yA
(3.26)
By Cramers rule, if
XB XA di x_xA cii xB xA X-XA
At 1 As
Vb VA -Cl2 y-yA d2 yB va y~yA
and A / 0, then s As/A and t = At/A. We have xp = (x',y') if the
intersection is between xA and xB (i.e., 0 < t < 1) and s < 0. Figure 3.6 shows
the other choices of xA and xB, which fail the conditions above.
46

3.2.4 Domain Splitting / Computing ViVJ,
Let Qe be an element such that Li is composed of two or three edges. The
intersection of two adjacent edges in Ye_ is called an inflow vertex. Since Ye_ is
non-smooth on the inflow vertices, we expect ViVJ to be discontinuous along the
characteristics starting at these points. We must partition Qe into subelements
Qf such that ViVJ. is continuous at each flf.
Suppose that there is a single inflow vertex, say x\. Let x5 be the intersec-
tion of the characteristic starting at xi and the element outflow boundary (Fig.
3.7). The line passing through x\ and *5, given by the equation
Qi(x,y) = (xs-x1)(y-yi)-(y5-yi)(x-xi)=0 , (3.27)
splits the vertices of fle in two subsets:
V1 = {xi = {xi,yi) ; qi(xi,yi) < 0} , (3.28)
V2 = {xi = (xi,yi) ; qi(xi,yi) > 0} . (3.29)
For the element shown in Fig. 3.7, we have V1 = {xi,X2,x5} and V2 =
{i, *5, x3, *4}. Notice that x5 could be relabelled as *2.5 to preserve the num-
Figure 3.6: Edges that do not contain xp. From the left to the right, we have
t > 1 in the first case, s > 0 in the second, and A = 0 in the third
47

Figure 3.7: One inflow vertex
bering order. The subelements are now defined by the vertices in V12. The
triangular subelements may be treated as degenerate quadrilaterals [57].
An edge e of a subelement Qf is defined as an external edge if e Otherwise, it is an internal edge. Internal edges, such as the one from X\ to 375
in Fig. 3.7, must not be taken into account for element boundary integrals.
When we have two inflow vertices, two characteristic curves split the vertices
into three sets. For the element shown in Fig. 3.8, V1 = {aq, *2, *5}, V2 =
, 375, Xqj 374} and V {374, 376, 373}
Figure 3.8: Two inflow vertices
Let us focus on N% x. Since N%>x{x,y) is constant within each O, it can be
determined by N%(x,y) and N%(xc,yc), where (xc,yc) is also an interior point
of Of.
48

Let x+ be the intersection of a horizontal line passing through the point x
with r*. (see Fig. 3.9). In particular, xc may be the midpoint of x and x+. It
follows that
neAx^)
N%{x,y) -N%(xc,yc)
X X r
(3.30)
Figure 3.9: The points x+ and xc used to evaluate N% x
3.2.5 Existence and Uniqueness of Solutions
We show in this section that a(-, ) is invertible in Kerand that the pair
Vh x Wh defined on Sec. 3.2 satisfies the conditions of Theorem 2.3 under the
following assumptions for any e = 1,..., nel (Fig. 3.10):
Assumption 3.1 : a is piecewise constant in Â£1 and a 0.
Assumption 3.2 : There exists an ordering of the elements Qe such that
1. rf c r+
2. Given an interior edge 7 = r+ n Tj_,
(a) i < j
(b) If x\ = 7 fl 7' with t'cT^, then x\ G Tl implies l > i
49

Assumption 3.3 : There are no opposite vertices of Qe aligned with a |ne.
Assumption 3.4 : There exist at least two outflow edges per element
Figure 3.10: A partition satisfying assumptions 3.1-3.3
Lesaint and Raviart [69] proved that there exists an ordering satisfying as-
sumptions 3.2.1 and 3.2.2a if the velocity field is constant. Assumption 3.3 has
been used in the analysis of the RFB method [20, 43] and implies that N% |pe
cannot be written with respect to traces of functions in Q\ (fle).
(3.31)
Lemma 3.1 : If Uh e KerB/t, then
a(uh, uh) = i J (ttfc ulf \a-n\dT + j u\a n dr
Proof : Let Uh Â£ KerBh be written as
UhW= ueP + aeN% uePÂ£Q1(Qe) e = 1,..., nel.
Let Vh G V* such that Vh |n== 2ueP, e = 1,..., nel. Since N% |n= 0,
Vh. |ri= 2Uh |n= (u/i uh + uh + Uh) |n
Therefore,
/ {uh uf)vha n dY = / {uh uf)[(uh uf) + (uh + u^)]a n dT
J n J ri
= f [{uh-ulf+ {u\-(ul)2)]a-ndT .
J re
50

Summing over 1 < e < nel and taking (3.18) into account yields
nel nel
(Uh~ uh?a n dT = ^2 / (uh ~ KT) |a n| dT . (3.32)
e=l e=l >T-
On the other hand,
nel
n dT
uei a
a(uh,uh) = Y] / a-Vuhuh>
e=l Jne
nel
= / *4
e=l
(nel - N
XJ yre (ft )2)a 'ndr + Jr ula'n dr
Substituting (3.32) into (3.33) yields (3.31).
Theorem 3.1 : Let uh e KerF^. We have = 0 if Uh satisfies
a(uh, vh) = 0 V vh E KerBh .
(3.33)
(3.34)
Proof : Let Uh G KerF^ satisfy (3.34). In particular, a(uh, Uh) = 0. From
Lemma 3.1,
[uft] = 0 on ri U \ T e = 1,..., nel ,
Uh ~ 0 on r_ U r+.
Let e = nel and
(3.35)
uh \ae= a\Nt + ... + ot%Nl + ae5NeE .
From Assumption 3.2.1, Uh |r.= 0, which implies a% = 0 together with As-
sumption 3.3. From Assumption 3.4, Uh |n vanishes on at least two edges.
51

Without loss of generality, assume that Uh |oe vanishes on 7, i = 2,3. Then,
Uh(x1) = 0, i = 2,3,4, so that Uh |n= a\N{. Let vj; G be defined as
/
fc() = <
Ni(x)
0
Since Te_ c 7 U 7I and JV| |7eU7e = 0,
cc G ,
x Qe .
[ K - n dr = [ N$vha n dT = 0 V vh G Vl i=i 'ri dri that is, G KerB^. Furthermore, 0 = o(fc, O = a\ (a VAT, jV3%e = ^ |a x (x2 4)| . From Assumption 3.3, |a x ( 4)| / 0. Thus, Uh |n= 0. Suppose that Uh\ne= 0 e > k + 1 . (3.36) From Assumption 3.2.2 and (3.35), we have Uh |r* = 0. Assuming 72 U 73 C r* and proceeding as above, we have 1^= a^N^. Consider the piecewise-linear function (Fig. 3.11) S() = N}{x) 0 x g {fi* e % ; aj| = xlj g T*} , otherwise . Since a;* G r*, vfc |n= 0 for any e < k from Assumption 3.2. Moreover, vfc |r_= 0, and since v\ is continuous, G KerBh. From (3.36), 0 = (*, vl) = Â£ (a Vu, vl) = a* (a VJVf, JV*)nl , e=l which implies uÂ£ |nt= 0. The proof follows from the induction principle. 52 Figure 3.11: A sketch of the function vfL Lemma 3.2 : Given s G So, let rs and Ws be defined by fs = {7j C f ; s < 0 at 7*} , Wa = {Ae5i; A |f\fs= 0} . If A G Wo satisfies [ vXsdf = 0 VfieQi(K) , Jr* (3.37) then A = 0. Proof : Let us introduce the following shape functions over an interval [1,1]: Mt) = pp Mt) = pp Let A G Ws satisfy (3.37). We first consider the case where fs is a single side % Then A G Ws can be written as Ai(t) = ai^i(t) + a24>2(t) , j = i , 0 , j Choose v in (3.37) such that v |^= p We have 0 = f (t) + a2(j>2(t)) dt (3.38) A1^. s l7i \ 3al + 3^2 J 53 Since s 0, (3.38) implies a2 = 2ai. Analogously, setting v |^.= (p2 yields cq 2a2. Thus, A = 0. Now assume that r,s is the union of two adjacent sides. Without loss of generality, let rs = 7i U 74. There exist constants a*, i = 1,2,4, such that A = (a\Ni + a2N2 + a4N4) |^; that is, / A |^. Ai (t) a\(j)\(t) + ct22{t) ^(t) + 0 3 = 1 > 3= 4 , 3* 1,4 Let s|^= Â§i, i = 1,4. Condition (3.37) can be written as follows: si J v(t-l)(iM0+a2fa(0)d+h J v(-l,7])(a1^1(rj)+a4h(v))dv = 0 . Choosing v = Ni, i 1,2,4, yields the equations 2(Si + S4)Q!i -\-Â§i(x2 +540:4 0 , < Si&i +2sia2 0 , +2s4Q:4 = 0 , which imply a* = 0, i = 1,2,4, and the result follows. Theorem 3.2 : If fj,h G Wh satisfies b(vh, nh) =0 V Vh G 14, then = 0. Proof: Let Â£ KerB^, e = nel and veh G Qi(Qe). Choose G Vh as 0l Vh i = e > vh In* < 0 , % 7^ e . Let G Qi(fie) such that Hh |ri= w^a n |r* and let A^ G L2(re) be 54 / K\rt-Wh\rt > < A^|r=\n= 0 . V From Prop. 2.1, there exists /io So such that 0 = / v^w^a ndT = / vAafao dT , J n where v = v^o F_1, A = A^o F-1 and d = (a n) o F_1. From Assumption 3.1, a n is piecewise constant on Te. Thus s = apL0 Â£ *% and A G Ws. Since was arbitrary, we have A = 0 from Lemma 3.2. Thus, p,/,. |re= 0. Proceeding as above with e = nel 1,..., 1 yields fih = 0. m 3.3 Numerical Experiments We compare the approximate solutions of (3.1) using methods DG, SDG, and DEM in two test problems. Let a = (cos(0), sin^))* 9 = 30, 45, 60 . We consider a domain D = [0, l]2 and two meshes of n x n elements, n = 20, 40, 80. One mesh (Ml) is composed of squares, while the other one (M2) is made up of non-rectangular quadrilaterals (see Fig. 3.12 and [7]). This setting aims to test the sensitivity of the methods to the discretization of the domain. Figure 3.13 shows the average CPU time spent on DEM and the SDG and DG methods. The computational cost of DEM is the highest because of extra degrees of freedom associated with Lagrange multipliers. 3.3.1 Advection Skew to the Mesh This test is based on the classical problem of propagation of a discontinuity from the inflow boundary without a source term [23]. A problem statement and 55 Figure 3.12: Meshes Ml (left) and M2 (right), n = 20 Figure 3.13: Average CPU time of DEM, SDG, and DG versus n the reference solution for 6 = 30 are shown in Fig. 3.14. The approximation errors in the L? norm are shown in Fig. 3.15-3.17. The estimate error curves suggest that DEM and DG/SDG have similar convergence rates. The errors for DG and SDG are also consistent with the experiments of Houston et al [56]. Figs. 3.18-3.20 present the graphs of the approximate solutions of DG, SDG and DEM along the cut x = 0.5 for the mesh Ml, with 9 30 and n = 40. All methods performed similarly, as seen in the estimate error curves of Fig. 3.15. 56 u = 0 v 0 x Figure 3.14: Statement and reference solution of Prob. 3.3.1, 6 = 30 Figure 3.15: Error plots of Prob. 3.3.1 (30), meshes Ml (left) and M2 (right) 3.3.2 Unit Force In the following problem we set g = 0 and / = 1. Its solution is (Fig. 3.21) u(x,y) = < y esc 6 if y < x tan 6 , x sec 6 otherwise . As in problem 3.3.1, DEM and (S)DG seem to have similar convergence rates (Fig. 3.22-3.24). When 6 = 45 DEM reproduces the reference solution at mesh Ml. Notice that the convergence rates for this problem are higher than the rates shown in Sec. 3.3.1. The higher regularity of the solution of Prob. 3.3.2 led to this improvement in the approximate solutions. 57 h h Figure 3.16: Error plots of Prob. 3.3.1 (45), meshes Ml (left) and M2 (right) h h Figure 3.17: Error plots of Prob. 3.3.1 (60), meshes Ml (left) and M2 (right) : : ; : 1 ! h ; Figure 3.18: DG solution for Prob. 3.3.1: uh (left) and (right) 58 r-N / 1 j 1 m Figure 3.19: SDG solution for Prob. 3.3.1: uh (left) and (right) Figure 3.20: DEM solution for Prob. 3.3.1: uh (left) and (right) 59 Figure 3.22: Error plots of Prob. 3.3.2 (30), meshes Ml (left) and M2 (right) Figure 3.23: Error plots of Prob. 3.3.2 (45), meshes Ml (left) and M2 (right) Figure 3.24: Error plots of Prob. 3.3.2 (60), meshes Ml (left) and M2 (right) 60 3.3.3 Elements with Three Inflow Edges When we consider problems 3.3.1-3.3.2 on the mesh M2 with 9 = 15, there are elements Qc such that Fe has three inflow edges, so that another enrichment would be needed (see Remark 3.6). The DEM approximate solutions with this extra enrichment (2ENR) are compared in Fig. 3.25 with solutions with the single enrichment (1ENR) given by (3.17). Notice that the additional enrichment did not significantly improve the accuracy the method. Figure 3.25: Error plots of Probs. 3.3.1 (left) / 3.3.2 (right), mesh M2 (15) 61 4. Advection-diffusion Equation We now consider the extension of (3.1) to the second-order problem of advection-diffusion: kAu + a Vn = f in Q , (4-1) u = g on T , K. where k, the diffusivity coefficient, is a positive constant. Problem 4.1 is an extension of the linear model of transport of a contaminant that admits an isotropic diffusion of the contaminant within the flow and represents a particular case of the linear, second-order partial differential equation -V (k(x, y)Vu) + a(x, y)Vu + a(x, y)u = f(x, y) . Due to diffusion, u{x) depends not only on its value along the streamline (characteristic curve), but on all points on a neighborhood of ax In particular, the concentration in the interior of Q, now depends on the concentration at the tangential and outflow boundaries. Therefore, one needs to prescribe u along the entire boundary T so that (4.1) is a well-posed problem. When k |o|, the solution is more affected by diffusion phenomena; that is, problem (4.1) is diffusion dominated. Similarly, if |a| >> ac, we have an advection-dominated problem. Given a characteristic length L of the domain, e.g., L = diam(Q), we define the (global) Peclet number as K 62 The Peclet number provides a measure of diffusive/advective dominance; the cases Pe 0 and Pe > 0 characterize a diffusive and an advective regime, respectively. Notice that problem (4.1) under the limits Pe > 0 and Pe - oo corresponds respectively to problems (2.14) and (3.1). Further, the order of the advection- diffusion equation does not change at the limit Pe Â¥ 0, but drops one unit if Pe oo. One should expect that the latter limit is more challenging. Indeed, problem (4.1) with 0 < \a\ < M and k > 0 is an example of a singularly perturbed problem with respect to k. A boundary-value problem is singularly perturbed with respect to a parameter e if its solution u(x, e) does not converge uniformly to u(x, 0). An asymptotic approach is considered in sections 4.3 and 4.7. In the following sections we consider some finite element approximations of problem (4.1). We focus on the advective regime. 4.1 Finite Element Methods There has been significant progress on continuous finite elements for the advection-diffusion equation. Since the approximate solution by the standard continuous Galerkin method may present spurious oscillations if advection dom- inates, some modification of this method must take place to ensure stability. Some of the proposed modifications are the introduction of mesh-dependent residual terms to the variational formulation (which we refer to as stabilized methods, see Sec. 4.1.1), enrichment of the polynomial field with residual-free bubbles [22, 20] or L-splines [55, 78], and the use of parameter-dependent meshes 63 that guarantee uniform convergence [70, 85, 92]. The extension of the DG method (3.10) for problem (4.1) is usually based on the weak continuity of the normal flux kVu n. One of the earliest works in this direction was done by Richter [84]. Baumann and Oden [10] extended the method (2.58) to the advection-diffusion equation in a similar manner. Stabilizing terms have been proposed to non-conforming [63, 72, 88], domain- decomposition [1, 2, 4, 94], and the Baumann-Oden method ([33] and Sec. 4.1.2). 4.1.1 Stabilized Methods A bilinear form g(-, ) defined on Vx V is coercive with respect to a norm ||-||v if a(v, v) > C \\v\\y for some C > 0. The coercivity, or stability of the Galerkin method in the Hl seminorm depends directly on the diffusivity coefficient. Let Vh be a finite-dimensional subspace of The Galerkin method for (4.1) with g = 0 is: find Uh G Vh such that a(uh, vh) = k (Vuh, Vvh)n + (o Vufc, vh)^ = (/, vh) V vh e Vh . (4.2) From Lemma 2.1, Stabilized methods modify the bilinear form a(-, ) to reduce the dependency of the coercivity on k, or in general, on the element Peclet number 2k Besides improving stability, the modification of a(-, ) should preserve consis- tency, i.e., the solution of the variational problem should reproduce the solution 64 of (4.2) if Vh is replaced by Hq(Q). For instance, consistency is preserved if one adds to o(-, ) a term that depends on the element residual acAUh + a V|, / of the solution, as in the SUPG (Streamline-Upwind/Petrov-Galerkin) method [23] nel a{uh, Vh) + ^2 (~K^uh + a vuh f, Te(a Vvh))ne = (/, vh) V vh e=l (4.3) and the GLS (Galerkin/least-squares [61]) method nel a{uk, Vh) + ^2 + a Vuh /, re(-tzAvh + a Vvh))ne e=l = (/, vh) Vvhevh. (4.4) These methods generalize the artificial diffusion and upwind finite difference schemes, which motivated a choice for the stabilization parameter [23]: = 2^ (th(Fe) + P?) ' Several variants of methods (4.3) and (4.4), and other choices for the stabi- lization parameter have been proposed (see [19] and references therein). 4.1.2 Stabilized Baumann-Oden Method The Baumann-Oden method consists of the DG method (3.10) plus the normal flux terms from method (2.58): nel a{uh,vh) (uh ul, vha w)^ e=l nel + ^2 ( (Uh> n)r n)Te ) (4.5) e=l = F(v) (ff, vha n)r_ + {g, {Vvfc} n)r VuG V* . 65 Engel et al [33] proposed a stabilized version of (4.5), which we refer to as SBO. This method includes an SUPG stabilization to a(*, ), such as in (3.11); a parameter a multiplying the second boundary integral in the LHS; and an interior penalty term rd ([/uj, [w/J)p. Note that the latter modification renders stability to (4.5) in the diffusive regime. If 7 = P fl P, i > j, we have (KJ,M)7 = (< = (ui - K)7 + K, v{)j = (uh-Kxtivh)inTi + {uh-uehxt,vh)inTj . The method SBO can be written as nel a(uh, vh) + ^2 ^ra (-kAu + a Vu, a Vv)Qe {uh u^, vha n)ri e=l + a (uft, K{Vvh} n)re (vh, k{S/uh} n)re + rd (uh ue^\ vh)re j nel = F(v) + ^2 To- (/, a Vn)ne (g, vha n)T_ e=l + a (g, {VtJfc} n)r + rd (g, vh)r V v e Vf . (4.6) Engel et al propose rd of order (9(/he), which is consistent with the interior penalty method [76]. 66 4.2 Discontinuous Enrichment Method We have developed discontinuous enrichment methods for advection and diffusion equations. The DEM for (4.1) benefits from both approaches, notably from the design a space of Lagrange multipliers. The DEM for the Poisson and advection problems weakly imposes the conti- nuity of Uh and Uhd-n respectively, across element interfaces. This suggests that one should try to balance these constraints for the advection-diffusion equation. For instance, we could weakly impose the continuity of kuh + u^a n. A first extension of the space Mh in (3.15) would be = JJspan{(/vu + va-ri) |nurg, v G Qi(Qe)} . e However, this space does not account for the domain outflow boundary. In particular, we add constant multipliers defined on T+: Wh = w; U W+ W+ = IJspanlw + va n|r^nr+, u G Qo(^e)} (4-7) e The enrichment function used for the advection equation is well suited for the advective regime, although it may not capture some features of second-order equations, such as boundary layers. For instance, consider Problem 4.6.1 below. Let us add the term kAu with k = 10-6 and outflow boundary conditions. Fig. 4.1 shows the DEM solution with homogeneous Neumann and Dirichlet outflow boundary conditions. The latter presented spurious oscillations in the neighborhood of the outflow boundary. We seek extra enrichment functions that can improve the DEM solution near boundary layers. The starting point is the asymptotic approximation described 67 Figure 4.1: DEM solutions for Prob. 4.6.1 (30) with Neumann (left) and Dirich- let (right) outflow boundary conditions in [62], which is reviewed for a particular problem. 4.3 Matched Asymptotic Expansion in the Unit Square Let a be a homogeneous velocity field a = (01,02)* with oi, 02 > 0. If Q C K has a boundary T of class C and /, g are also of class C, then one can construct a family of asymptotic series of the form n u(x,y) = lim Un(x,y) Un(x,y) = 'S^Ktui(x,y) (4.8) n-Aoo * ^ i=0 satisfying ||CUn f\\ = (P(/cn+1), where Â£ is the advection-diffusion operator [62]. Substituting (4.8) into the equation Cu = / and collecting the terms of same order in k yields the recurrence equations a Vu0 = / , (4.9) a Vui = Auj_i % > 1 . (4-10) Since the terms Ui(x7y) are solutions of first-order problems, we cannot impose boundary conditions on Ui along the entire boundary. Therefore, the asymptotic approximation (4.8)-(4.9) to (4.1) may not satisfy u = g. 68 The method of matching asymptotic expansions ([62, 73, 86] and the refer- ' ences therein) proposes the construction of additional asymptotic series defined on local coordinate systems near the outflow boundary. These series are re- quired to satisfy the differential equation under the new coordinates and some matching boundary conditions satisfied by the superposition of the series. We consider a zeroth-order asymptotic approximation u0(x, y). In this case, we can relax the smoothness of T. Given <5 > 0, we consider in particular the domain shown in Fig. 4.2, which excludes the corners of K. Figure 4.2: A domain with (^-boundary Let us add to (4.9) an inflow boundary condition: a Vo = / in Uq = g on (4.11) Problem (4.11) is referred to as the reduced problem [86]. Since u0 ^ g on r+ in general, we have to introduce other asymptotic expansions near T+. For instance, let Q1 = {(x,y) G ; |1 x\ < J and \y\ < 1 5} (Fig. 4.2) and vg G C1^1) such that uff(l 8,y) = 0 and vg(l,y) = g(l,y) uQ(l,y). We 69 seek a solution v for the problem / r1 kAv + a Vv = 0 in Q1 , v vg on T1 . (4.12) Note that u = uo + v is continuous at x = 1 5 and satisfies u = g at nr+. Consider the local coordinates C= (1 -x)/k , < = y , which maps I21 onto = (0, 2/k) x (1+5,1 5). Define u;(Â£, (j>) = u(x(C), y{4>)) and = vg(x(Q,y((j))). From (4.12), the function w satisfies ~ KWM ~ + = 0 rv r\i W,CC + aiw,C + !^wM ~ Ka2Wt4> = 0. (4.13) Therefore, w is the solution of the following problem: w, w = g on T . In general, substituting OO wiC, 4>) = ^KlWi{c, i=0 into (4.13) yields r ^o,cc + i'u;o,c = 0 j < wljCC + aiWi^ ozwoj , wi,CC + alwiÂ£ = ~wi-2, + 02^1-1,^ > %>2. (4.14) 70 Since we restricted ourselves to terms of order zero, we consider w = Wo, which is approximated as follows: For each 4> G (1 + <5,1 5), w0 satisfies w0,cc + io,c =0 , < w0(0 , limty0(C,^) =0; C->o that is, wo((,(j)) = <7(0,4>)e~ai^. Transforming back to (x, y), we approximate v in fT by vi{x,y) = vg(l,y)eaiix~1)/K = (g{l,y) u0(l,y))eai{x~1)/K . (4.15) Analogously, consider Q2 as in Fig. 4.2. We have that v2(x, y) = {g{x, 1) u0(x, l))ea2(2/_1)/K (4.16) constitutes an approximation in Q2. Extending the domains of vi>2 to Cl, we have that the function Ksixiy) = uo(x,V) + (0(1,1/) o(l,2/))eai(l-1)/K + (g(x, 1) u0(x, l))ea2{y~1)/K . (4.17) is an approximation of order 0(k1) of the solution of (4.1) on Cl excluding the regions near the corners. Indeed, ||Â£ti*s /|| = k ||Au0|| < Ck for some C > 0 independent of k and u*as g on T_ U (r1 fl T+) U (r2 fl T+). Remark 4.1 : Asymptotic expansions can be adapted for domains with corner points [62, Sec. IV.2]. In particular, the solution can be improved by adding a corner layer correction: uas(x, y) = Ksix, V) ~ (0(1,1) o(l, l))e(ai(*1)+0a(y_1))/,s . (4.18) 71 4.4 Boundary-layer Functions Let us denote iV| by NE1 from here on. Motivated by the former section, we propose the following enrichment functions for rectangular elements: . ne,3(x,v) = exP where Xq = (xq, yffi satisfies a Xq> a x V x E Qe [35]. Roos et al [86] refer to N% 2 ari(l as boundary-layer functions. In exponential-based finite-difference schemes [3, 36, 55], these functions solve the corresponding one-dimensional equations in the x and ydirections. Indeed, N% 2 and AT| 3 are solutions of the homogeneous equation kAu + a V = 0 , (4-19) if we assume u(x,y) = u(x) and u(x,y) = u(y), respectively. In the context of the discontinuous enrichment method, the boundary-layer functions are also free-space solutions of the advection-diffusion equation. We expect the polynomial field and the first enrichment N%^ to approximate the reduced solution uo(x, y), while the new enrichment functions should approx- imate the extra terms in (4.17). Since the functions v\ and u2 in (4.15)-(4.16) are approximations in the neighborhood of the outflow boundary, we restrict N% 2 and JV|; 3 to elements Qe along the outflow boundary. One could try to capture the corner layer correction by adding the function f(x,y) = exp ^ a<)j . (4.20) NeE>2{x,y) exp ai(x K 72 However, the L2-norm of the trace of f(x, y) over Cle rapidly converges to zero. Furthermore, the extra term in (4.18) is a punctual correction. Neverthe- less, / has been used to derive exponentially-fitted shape functions (L-splines, [78, 55]) and a stabilization parameter which depends on the direction of a [54]. 4.5 A Modified Gaussian Quadrature We would like to employ tensor-product quadrature rules to compute inte- grals such as I = j J d^dr] , where q is a polynomial of degree < 3. Let fji, Wi be the 2-point Gauss-Legendre quadrature weights and nodes, respectively. We have 2 I = f dÂ£ . i=l The integration with respect to f cannot be accurately performed with the Gauss-Legendre rule, specially if a\/K is large. We need a quadrature rule for /I _ . K /2ai/re p{0e~~^ dÂ£ = / p( 1 Kt/ai)t i fli Jo e t dt , where p is a polynomial of degree < 3. One can find t\$, wi;2 such that
r*2ai/K
p(t)e_t dt = wip(ti) + w2p(t2) ,
/
which implies
K [ / ' h\ (- h
lexp Wip 1 K + W2p 1 K
ai L V aj i, ai
= wTPp{iTP) + Â£*?(Â£*) ,
CexP SI,2 = 1- 3? 1 5*+ I'tO >
a i
. exp Wl,2 W 1,2 = K .
ai
73

If ai/rc > oo, t1)2 and iuij2 are the Gauss-Laguerre parameters [30]
d,2
2^pV2 and w^
2 \/2
We can approximate and in the advective regime by
Â£exp 2=F\/2 j ~exp 2y/2
f i 9 = 1 /c---- and w,, = k-------.
M oi 4ai
Fig. 4.3 compares the values of ^xp and E^xp for 1 < cli/k < 103.
Figure 4.3: |fcp and Â£{xp for 1 < a\/K < 103
4.6 Numerical Experiments on Rectangular Domains
Four problems are tested in this section. The first two are extensions of
the examples from the previous chapter, followed by two standard test problems
[39]. The approximation errors are calculated on a mesh of n x n rectangular
elements, with n = 20, 40, 80.
DEM is compared with SBO (Sec. 4.1.2) and the GLS method with the
stabilization parameter defined in [44]. Figure 4.4 shows the average CPU time
74

spent on each method. The computational cost of DEM is the highest because
of extra degrees of freedom associated with Lagrange multipliers. Also, the
Figure 4.4: Average CPU time of DEM, SBO, and GLS versus n
4.6.1 Advection Skew to the Mesh
We modify Problem 3.3.1 by adding the term kAu with k = 10_m,
m = 0,3,6 to (3.1). We set homogeneous Dirichlet conditions in the outflow
boundary, which leads to boundary layers as k decreases.
Figures 4.5-4.7 show how each method resolves the boundary layer when
9 = 30 and h = 1/40. In particular, they highlight a major deficiency of dis-
continuous Galerkin methods regarding outflow boundary conditions. Baumann
and Oden [10, Sec. 4.1] reported similar results for a one-dimensional problem.
The GLS method was able to approximate the boundary layer and preserve
stability. However, the accuracy drops from the interior to the boundary of the
domain. DEM, which was initially defective (Fig. 4.1), achieves higher accuracy
with the aid of the boundary-layer functions.
75

y y y
Figure 4.5: DEM elevation plots for Prob. 4.6.1, n = 40, 6 = 30, k = 10-m,
m = 0,3,6 (left to right)
Figure 4.6: SBO elevation plots for Prob. 4.6.1, n = 40, 6 = 30, k = 10 m,
m = 0,3,6 (left to right)
Figure 4.7: GLS elevation plots for Prob. 4.6.1, n 40, 0 = 30, k = 10 m,
m = 0,3,6 (left to right)
76

When k = 10-6, the diffusion effect is negligible (except at the outflow
boundary) and we can consider the reference solution u of Prob. 3.3.1 with the
additional constraint u |r+= 0. Figs. 4.8-4.9 display the errors of the solution
of each method with respect to this reference solution. Although the errors of
DEM and SBO were close, one should note that the error of SBO at the outflow
boundary (see Fig. 4.6 (right)) is difficult to detect in the L2 norm.
4-6.2 Unit Force
Problem 3.3.2 is modified in the same fashion as in Problem 4.6.1. Figs.
4.11-4.13 confirm the difficulties of the SOB and GLS methods in resolving the
boundary layer. We consider again the errors with respect to a reference solution
As in Prob. 3.3.2, the errors for 6 = 30 and 6 = 60 are identical, so the
latter is omitted. Notice that DEM achieved superconvergence for 9 = 45, as
in Prob. 4.6.1. The L2-norm of the GLS error was dominated by the error at
the strips of elements along the outflow boundary, which is of order 0(h^2).
77

h
h
Figure 4.8: Error plots of Prob. 4.6.1, 9 = 30, 45
h
Figure 4.9: Error plots of Prob. 4.6.1, 9 = 60
Figure 4.10: Error plots of Prob. 4.6.2, 9 30, 45
78

Figure 4.11: DEM elevation plots for Prob. 4.6.2, n = 40, 6 = 30, k = 10 m,
m = 0,3,6 (left to right)

Figure 4.12: SBO elevation plots for Prob. 4.6.2, n = 40, 6 30, k = 10 m,
m = 0,3,6 (left to right)
Figure 4.13: GLS elevation plots for Prob. 4.6.2, n = 40, 6 = 30, k = 10 m,
m = 0,3,6 (left to right)
79

4.6.3 Advection in a Rotating Flow Field
Consider a homogeneous Dirichlet Problem in Q = [0.5,0.5]2 without dis-
tributed sources (/ = 0), k = 10-6 and a (y,xY (Fig. 4.6.3).
Figure 4.14: Statement of Problem 4.6.3
There is an internal boundary along the negative y-axis, with the boundary
condition u(0,y) = w{y) = [cos(47ry + 7r) + l]/2. However, we weakly enforce
u(0,y~) = w(y) and [u(0, y)] = 0 in DEM and SBO.
We employed the reference solution from [54]. Figs. 4.15-4.16 indicate that
all methods had a similar performance, which is probably due to the regularity
of the solution. In this situation, the fastest method (GLS) is more appropriate.
Figure 4.15: DEM, SBO and GLS elevation plots for Prob. 4.6.3 (left to right)
80

Figure 4.16: Error plot of Prob. 4.6.3
4.6.4 Thermal Boundary Layer
Consider the boundary value problem outlined in Fig. 4.17, with k = 7 x
1CT4 and / = 0. The main difficulties are the formation of an outflow boundary
layer at x = 1 and a parabolic layer [86] along y = 0 (Fig. 4.18).
y u = 1
J a=(2y,0)f 0.5 u = 2y Li .i
X

u ~ 0 i ,
Figure 4.17: Statement and mesh pattern for Prob. 4.6.4
The errors shown in Fig. 4.19 were calculated with respect to the GLS
solution in a mesh of 200 x 200 elements. The drop of convergence rates for
SBO and DEM might be related to difficulties with the parabolic layer.
We generalize the results of Sec. 4.3 to a quadrilateral element Qe. Consider
the mapping F : K > Vte defined in Sec. 2.4 and suppose that Â£ = 1 is mapped
81

Figure 4.18: DEM, SBO and GLS elevation plots for Prob. 4.6.4 (left to right)
Figure 4.19: Error plot of Prob. 4.6.4
to an outflow side of fi. Let D and Q1 be the images of f2g and S71 (Fig. 4.2)
under the mapping F.
Let u0 be the solution of the reduced problem (4.11), r* = Te>1 D r+, Tf C
T61 \ rÂ£ and vg 6 C1^1) such that
vg |ri= 9 ~ uo and nff|ri=0.
As in (4.12), we define a local correction based on the solution v of the
problem
{kAv + a Vn = 0 in fP ,
v = vg on T1 .
82

^9(,n) = vg(x(,r}),y(tri))wdv(,Ti) = v(x(,Tj),y(,ri)). Prom (2.33)
and (2.46), v is the solution of
+ KAzVfy + kAsV^-'c
{B1 + kB3)vÂ£ + (B2 + kB4)v>T] = 0 in K , (4.21)
v = g on f ,
where
AM.n) =-&\$ Ai((,v) = 2x*x+?4V" M(,n)
b^v) =aiVx,j?2X", m,v) =a2X*
d + vl
i2
B,(lv) =

|J|
l
xxxAn + VxVAn + ^jf(xAx,v + VÂ£V,n) + ^Jjf(xl + y,i) )
\J\
21 xaxav + vaVav + Tjf Ow*? + VAV*) + T/f + vl) )
Let us consider the local coordinates
C=(i-0/ ,
= v-
Since the coefficients Ai, Bi of equation (4.21) depend on Â£ = 1 k(, we
must consider the expansion of these coefficients with respect to k. For instance,
the Taylor series for A\ (Â£, rj) centered at Â£ = 1 is given by
i= 0
If Ai^, i > 0 is defined as

d
83

then Ai(C, (f>) = -Ai (Â£(Â£), r}{4>)) can be written as
i=0
The functions A2,..., B2 are represented in a similar way. Let g((, (f>) =
5(Â£(C)>7The zeroth-order expansion for to(C, (f>) = w(Â£(C)i ^(0)) satisfies
-A.1,0{i>)wo,cc + Bito((f>)w0^ 0 .
Proceeding as in Sec. 2.4, we approximate w0 by the solution of
/ _
-4i,o(0)wo,cc + B1)0(^)wo,c = 0 ,
< w0(0,(f>) =g[0,) ,
lim wo((, ) =0 ,
C->
whose solution is
wo(C, 4>) = -5(0,4>) exp ( -
The corresponding function in the (Â£, fj) variables is
The approximations near the remaining sides of K can be found in a similar
way. These boundary-layer functions are used as enrichment functions in the
following section.
vi&v) = 9(1, V) exp
-a23v>(M)
x?v(1,v)+ y^v)
-l)
4.8 Generalized Boundary-layer Functions
Let us generalize the enrichment field to non-rectangular elements. For each
Si E re_ n r_, we add the enrichment function N% k(x,y) Nl(Â£(x,y)1r](x,y)),
where
84

exp
\J(l,ri)\a1yjt](l,Tj) a2xiTI(l,rj)
(CTl)
K z*(i ,v) + y%(hv)
i % = 3 4-1 ,
Al(t,v) = {
(4.23)
exp
[J(Â£i)l Q2^,e(^,i) Qi2/,e(Â£,i)
zf (Â£, 1) + y%(Â£, 1)
(v T 1)
, t = 2l
The functions N^ k, k > 1, can be seen as asymptotic approximations of
free-space solutions and are conveniently defined in the reference element.
To employ the quadrature rule suggested in Sec. 4.5, one has to evaluate the
coefficient of the exponential at the quadrature points of the alternate direction.
For instance, we can write as exp(a(Â£ 1)/k), where the constant a
depends on the quadrature point ry in the direction of rj.
4.9 Numerical Experiments on Quadrilateral Meshes
Let us consider Prob. 4.6.2 (k = 106) in O2 = [0, l]2 with the mesh M2
from Sec. 3.3 and in the domain 03 (Fig. 4.20) defined by the union of the sets
{(, y); -y/4 < x < 1 + y/A, y e (o, l)} ,
{(x,y) ; 1 < y < 1 + (Ax + 1)(5 4x)/144, x G (-1/4,5/4)} .
The domain f23 tests the performance of the boundary-layer functions near
boundaries which are not aligned with the coordinate axes. For convenience,
we refer to the meshes associated with f23 as M3. Figs. 4.21-4.26 display the
elevation and contour plots in this domain for each method.
85

Figure 4.20: Domains Q2 and 03 with respective meshes when n = 20
Ideally, the approximation errors on mesh M2 should be the same as in Prob.
4.6.2. We compare the errors on mesh M2 with those from Fig. 4.10 (unmarked
lines on Figs 4.27 4.29). Except for the loss of superconvergence of DEM at
9 = 45, the results are not too sensitive to the perturbation of the mesh. Note
also that the convergence rates on f23 were close to those on a square domain.
4.10 Limitations of DEM
As seen in Prob. 4.6.4, DEM does not perform well on problems with
boundary layers that are aligned with the velocity field.
Another adverse situation is when the inflow boundary of an element Vte is
a single edge, which is not covered by assumptions 3.1-3.3. In this case, the
enrichment N%^{x,y) is not included (see Sec. 3.5). If Vte is an interior element,
then no enrichments would be included at De.
These limitations are apparent in the problem described in Figs. 4.30 4.32.
The stabilization terms guarantee the robustness of the GLS and SBO methods
for k = 10m, m = 0,3,6, while the DEM solution becomes unstable as k
decreases (Figs. 4.33 4.41).
86

Figure 4.21: DEM elevation plots for Prob. 4.6.2 (Q3), 9 = 30, 45, 60 (left
to right)
Figure 4.22: SBO elevation plots for Prob. 4.6.2 (fi3), 9 = 30, 45, 60 (left to
right)
Figure 4.23: GLS elevation plots for Prob. 4.6.2 (D3), 9 = 30, 45, 60 (left to
right)
87