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
University of Colorado at Denver
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
(DEM) for the advection equation, the advectiondiffusion equation and the
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 residualfree bubbles (RFB). The
enrichment functions derived herein were crucial in capturing boundarylayers on
advectiondiffusion problems and stabilizing a loworder velocitypressure 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
Farhat, Lutz Tobiska and Alexandre Madureira for helpful discussions/comments
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. Linear Advection Equation.......................................... 35
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. Advectiondiffusion Equation.................................... 62
4.1 Finite Element Methods......................................... 63
4.1.1 Stabilized Methods............................................. 64
4.1.2 Stabilized BaumannOden Method................................. 65
4.2 Discontinuous Enrichment Method................................ 67
4.3 Matched Asymptotic Expansion in the Unit Square................ 68
4.4 Boundarylayer 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 Boundarylayer 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 BaumannOden Method..................................... 97
5.3 Discontinuous Enrichment Method.................................... 99
5.3.1 FreeSpace 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 Velocitypressure Field.................................. 119
5.4.2 Driven Cavity Flow.............................................. 121
5.4.3 The Stickslip 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 nonmortar 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.13.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 = 10m,
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 = 10m,
m = 0,3,6 (left to right)............................................ 79
4.13 GLS elevation plots for Prob. 4.6.2, n = 40, 9 = 30, k 10m,
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, steadystate NavierStokes 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 103......................... 92
4.37 SBO elevation plot for Prob. 4.10, k = 103.......................... 92
4.38 GLS elevation plot for Prob. 4.10, k = 103.......................... 92
4.39 DEM elevation plot for Prob. 4.10, k = 106......................... 93
4.40 SBO elevation plot for Prob. 4.10, k = 106.......................... 93
4.41 GLS elevation plot for Prob. 4.10, k = 106.......................... 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 stickslip 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 NavierStokes 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 NavierStokes 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 NavierStokes
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 firstorder 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:
Residualfree 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
Twolevel FEM. An approximate RFB method where the residualfree 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 underresolved fine scales, such as
Stabilized methods. The Galerkin variational formulation is modified with
meshdependent 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:
Linear advection equation;
Linear advectiondiffusion equation;
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 secondorder 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 adhoc
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, firstorder 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 advectiondiffusion in Chapter 4. The
method is enhanced with exponential boundarylayer functions formerly
derived by the method of asymptotic expansions [62].
The objective of Chapter 5 is to stabilize the bilinearvelocitypressure
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.12.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 nonhybrid 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 Lipschitzcontinuous boundary [27]. We
represent the boundary of Cl by T. This thesis is particularly concerned with
twodimensional problems (n = 2) and realvalued 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 = supu() .
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^ = n0 K and Wu^ = little n.
All definitions above may be extended to vectorvalued 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 seminorm 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 + IV'I2)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'lli/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 R1/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) (27)
e=l
Given a fixed order of the elements Qe, we define the jump {v} of v G X as
[I(a5) = v()nt)(aj)n# ccGfTnfV i>j. (2.8)
We also define
nel nel
= , fJre , = .
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 saddlepoint 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 nonhomogeneous 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 Py
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 finitedimensional 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 ~ Pfcg\Kerfl < C(/f Vh\\V + IIP ~ ftllg )
(2.22)
Condition (2.22) is known as infsup or BabuskaBrezzi 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 Py
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
interpretation of the Lagrange multiplier (see also [18, Prop. II.l.l]):
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{llgldiv ;getf(div;0) qn = 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 finitedimensional 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 ith side
and the ith 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 rV2yi
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 firstorder 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)
Secondorder 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,
dJ1 dJ1
B , C x .
dr]
21
Therefore,
^,ij WJtlAnAkj + AuBkjU^k + A^iCkj'Qjk (2.35)
Let J1^ be the ftth row of J L We write (2.35) in vector notation as
H(u) = J'iJ(w)J~l+J1(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,,kiAiiAki 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?
VMfirFi + x>(x*
\J\
+
\J\,v \ ( \JU
y,z 1 vtn y,srf\ + 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^voF1 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 infsup condition was proven
[81, Lemma 10] for meshdependent 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 nonconforming method. In fact, several nonconforming methods
can be written as conforming, primal hybrid methods [81].
Another approach is to approximate the saddlepoint 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 nonconforming domain decomposi
tion method that weakly imposes continuity between subdomain interfaces by
introducing a finitedimensional 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 nonoverlapping subdomains Dk, 1 < k < Np\
that is,
~Nd
n = (J Dk.
k=1
Let Xh be a finitedimensional 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 onedimensional piecewisepolynomial 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 (259)
rfcJes
Consider the following bilinear forms:
nd
aD(u,v) = ]T(Vu,V?V , (2.60)
k=1
bD(v,fj) = (M>/*)rw (261)
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; ne P,(f2e) V Qe G D*, n U*nr= 0} ,
Mh = {v> G Wfc' ; ^l7e Ai(t) if e 7 or G 7> ,
W? = 0* e C(f{H}); /x7G 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 finitedimensional 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 finitedimensional 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 finitedimensional
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 freespace 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 nondiffusive 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 firstorder, 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 righthand 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 welldefined traces of restrictions of Vh to each element of % These
approximate traces are known as numerical fluxes in the discontinuousGalerkin
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
advectionreaction 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) <39)
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) elementby
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(uhuh, 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 highorder 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
advectiondiffusion 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 nonmortars. 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
nonmortars
Figure 3.3: Mortar and nonmortar 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 secondorder problems [35]. Since we are
interested in advectiondominated problems, it is desirable that the method
be able to solve firstorder problems as well. In particular, problem (3.1) is
the hyperbolic limit of the advectiondiffusion 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 advectiondominated advectiondiffusion
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 tdt0 ,
&(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.13.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 advectionreaction 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
aVueh
+ (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 that $2 = 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) =
xxp yyp
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 = yyA
(3.26)
By Cramers rule, if
XB XA di x_xA cii xB xA XXA
At 1 As
Vb VA Cl2 yyA 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
nonsmooth 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) = (xsx1)(yyi)(y5yi)(xxi)=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.13.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 \an\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 [{uhulf+ {u\(ul)2)]andT .
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] / aVuhuh>
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 piecewiselinear 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(tl)(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\rtWh\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 F1 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 nonrectangular 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.153.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.183.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.223.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
rN
/ 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.13.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. Advectiondiffusion Equation
We now consider the extension of (3.1) to the secondorder problem of
advectiondiffusion:
kAu + a Vn = f in Q ,
(41)
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, secondorder 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 wellposed 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
advectiondominated 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 boundaryvalue 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
advectiondiffusion 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 meshdependent
residual terms to the variational formulation (which we refer to as stabilized
methods, see Sec. 4.1.1), enrichment of the polynomial field with residualfree
bubbles [22, 20] or Lsplines [55, 78], and the use of parameterdependent 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 advectiondiffusion equation in a similar manner.
Stabilizing terms have been proposed to nonconforming [63, 72, 88], domain
decomposition [1, 2, 4, 94], and the BaumannOden 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 finitedimensional 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 (StreamlineUpwind/PetrovGalerkin) 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/leastsquares [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 BaumannOden Method
The BaumannOden 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
= (uhKxtivh)inTi + {uhuehxt,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 Uhdn respectively, across element interfaces. This suggests that
one should try to balance these constraints for the advectiondiffusion 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 + vari) 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 nr^nr+, u G Qo(^e)} (47)
e
The enrichment function used for the advection equation is well suited for
the advective regime, although it may not capture some features of secondorder
equations, such as boundary layers. For instance, consider Problem 4.6.1 below.
Let us add the term kAu with k = 106 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)
nAoo * ^
i=0
satisfying CUn f\\ = (P(/cn+1), where Â£ is the advectiondiffusion 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 . (410)
Since the terms Ui(x7y) are solutions of firstorder 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 zerothorder 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Â£ = ~wi2, + 02^11,^ > %>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(l1)/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 Boundarylayer 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 boundarylayer functions. In
exponentialbased finitedifference schemes [3, 36, 55], these functions solve the
corresponding onedimensional equations in the x and ydirections. Indeed,
N% 2 and AT 3 are solutions of the homogeneous equation
kAu + a V = 0 , (419)
if we assume u(x,y) = u(x) and u(x,y) = u(y), respectively. In the context
of the discontinuous enrichment method, the boundarylayer functions are also
freespace solutions of the advectiondiffusion 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 L2norm 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 exponentiallyfitted shape functions (Lsplines,
[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 tensorproduct 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 2point GaussLegendre
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
GaussLegendre 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 GaussLaguerre 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
enrichment functions require additional quadrature rules.
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.54.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 onedimensional 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 boundarylayer functions.
75
y y y
Figure 4.5: DEM elevation plots for Prob. 4.6.1, n = 40, 6 = 30, k = 10m,
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 = 106, 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.84.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.
46.2 Unit Force
Problem 3.3.2 is modified in the same fashion as in Problem 4.6.1. Figs.
4.114.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
adapted from the advective problem for k = 106 (Fig. 4.10).
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 L2norm 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 = 106 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 yaxis, 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.154.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.
4.7 Asymptotic Expansion in Quadrilaterals
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 nffri=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=(i0/ ,
= 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 zerothorder 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 boundarylayer 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 Boundarylayer Functions
Let us generalize the enrichment field to nonrectangular 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 41 ,
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
freespace 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 boundarylayer 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.214.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.13.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
