A NEW SUBMESH STRATEGY IN THE
TWOLEVEL FINITE ELEMENT METHOD FOR
THE ADVECTIVEDIFFUSIVE EQUATION
by
FengNan Hwang
B.S., FuJen Catholic University, Taiwan, 1995
A thesis submitted to the
University of Colorado at Denver
in partial fulfillment
of the requirements for the degree of
Master of Science
Applied Mathematics
1999
This thesis for the Master of Science
degree by
FengNan Hwang
has been approved
by
.ndrew Knyazev
ii/i 6/
Date
Hwang, FengNan (M.S., Applied Mathematics)
A New Submesh Strategy in the TwoLevel Finite Element Method for the
AdvectiveDifFusive Equation
Thesis directed by Professor Leopoldo P. Franca
ABSTRACT
The residualfree bubble (RFB) method for the general boundary value prob
lem is reviewed. The two level finite element method (TLFEM) derived from
the RFB method and its application to the advectivediffusive equation are
addressed. There is a boundary layer near outflow boundaries for the bubble
shape functions when the flow is advective dominated. In this work, we in
troduce a new submesh strategy depending on the direction of the flow in the
TLFEM and design an algorithm for generating this submesh. The numerical
results show that the new submesh is able to capture the boundary layer which
is caused by the choice of bubble functions. The effect of a more accurately
approximated residual free bubble function produces some improvements in
the solution of the advectivediffusive equation.
This abstract accurately represents the content of the candidates thesis. I
recommend its publication.
Signed
Leopoldo P. Franca
111
DEDICATION
To my parents and TzuChia
ACKNOWLEDGMENTS
First, I would like thank my advisor, Professor Leo Franca, for his in
structive guidance and wonderful patience during the preparation of this thesis.
Also thanks to the members of our research group, especially Ali Nesliturk and
Saulo P. Oliveira, for their technical support on programming. Finally, thanks
goes to Randy Chase for his linguistic assistance.
CONTENTS
Figures ............................................................ vii
Chapter
1. Introduction....................................................... 1
2. The review of RFB for boundary
value problems .................................................... 5
3. The TLFEM for the advectivediffusive problem..................... 11
4. The algorithm for nonuniform submesh generator.................... 24
5. Numerical results................................................. 31
5.1 A problem with discontinuous
boundary condition.............................................. 31
5.2 Thermal boundary layer problem.................................. 38
5.3 Bubble ramp problem ............................................ 47
6. Conclusion........................................................ 50
References........................................................... 52
vi
FIGURES
Figure
3.1 Obtaining x~ and x+ from any x G K ............................. 16
3.2 Problem statement............................................... 17
3.3 Nodal points and sides numbering................................ 17
3.4 Submesh refines on outflow boundaries........................... 19
3.5 The element parameter computation............................... 21
3.6 Submesh used for different cases: T indicates inflow boundary
and O indicates outflow boundary............................ 23
4.1 Mesh grading from left to right: a quadratic transformation . 25
4.2 Local node and side ordering.................................... 27
5.1 Domain for skew convection problems with two different velocity
fields.......................................................... 32
5.2 Two different 10x10 submeshes used for TLFEM: uniform (left)
and nonuniform (right).......................................... 33
5.3 Comparison of bubble shape functions by using two different
submeshes for the 45 problem: uniform (left) and nonuniform
(right)......................................................... 34
5.4 Comparison of bubble shape functions by using two different
submeshes for the 60 problem: uniform (left) and nonuniform
(right)......................................................... 35
vii
5.5 Comparison of TLFEM solutions by using two different sub
meshes for the 45 problem: uniform (left) and nonuniform (right) 36
5.6 Comparison of TLFEM solutions by using two different sub
meshes for the 60 problem: uniform (left) and nonuniform (right) 36
5.7 Comparison of GLS and TLFEM solutions for the 45 problem:
GLS (left) and TLFEM (right) ............................ 37
5.8 Comparison of GLS and TLFEM solutions for the 60 problem:
GLS (left) and TLFEM (right) ............................ 37
5.9 Problem statement for the thermal boundary layer problem ... 38
5.10 Mesh used for the thermal boundary layer problem......... 39
5.11 Submesh used for the thermal boundary layer problem...... 39
5.12 Comparison of TLFEM solutions using two different submeshes 40
5.13 Bubble shape functions for the diffusive dominated case: unifor
m submesh (left) and nonuniform submesh (right).......... 41
5.14 Bubble shape functions for the advective dominated case: uni
form submesh (left) and nonuniform submesh (right).......... 42
5.15 Bubble shape functions approximated by two different submesh
es for the bottom region: uniform (left) and nonuniform (right) 44
5.16 Bubble shape functions approximated by two different submesh
es for the top region: uniform (left) and nonuniform (right) ... 45
5.17 Comparison of the TLFEM solutions by using two different sub
meshes for the thermal boundary layer problem when k = 10~6:
uniform (left) and nonuniform (right).................. 46
viii
5.18 Comparison between the GLS method and the TLFEM method
for the thermal boundary layer problem when k = 10 6: GLS
(left) and TLFEM (right) ...................................... 46
5.19 Problem statement for the bubble ramp problem................. 47
5.20 Comparison of bubble shape function
ent submesh for the bubble ramp problem: uniform(left) and
nonuniform (right)...................................... 48
5.21 Comparison of the GLS and the TLFEM by using two different
submeshes for the bubble ramp problem..................... 49
IX
1. Introduction
We approximate the advectivediffusive equation in this work. There
are a lot of applications in fluid dynamics using this model as well as in the
semiconductor industry. For the advectivediffusive equation, the behavior of
the solution depends on the magnitude of the velocity field and the diffusivity
coefficient.
Let us define the meshPeclet number as Pe =
\a\h
2re
Here, a is the
velocity field, h is the mesh size and re is the diffusivity coefficient. When
this number is large, we say the model is advective dominated; otherwise, it
is diffusive dominated. In the finite element literature, it is well known that
the standard Galerkin method using piecewise linears performs poorly for the
advective dominated model. Spurious oscillations are frequently detected in
the solution.
In order to overcome this difficulty, stabilized finite element methods
have been introduced [3, 5, 10], and one version was denoted by the Galerkin
leastsquares method (GLS). The GLS method adds an artificial term to the
variational formulation. This additional term not only improves the numerical
stability of the Galerkin method but also preserves good accuracy.
1
On the other hand, it is not hopeless to apply the standard Galerkin
method to solve this problem. Herein we are interested in approximating the
residualfree bubble (RFB) method using the Galerkin method. The bubble
function is chosen so that the computed solution satisfies the original differen
tial equation in the interior of each element and vanishes on the boundary of
each element.
In [1] a relationship is established between the stabilized finite ele
ment method and the Galerkin method using piecewise linears enriched with
one bubble per element for the advectivediffusive model. Therein these two
methods are shown to be equivalent for diffusive dominated models and ap
propriately defined element parameters hk Later Brezzi, Franca and Russo
[2] proved the coercivity /stability inequality for a limiting case of the RFB
method. Further progress has been made by Franca, Nesliturk and Stynes in
[7], where they obtained the desired stability condition for positive but small
diffusivity parameters under some hypotheses that the flow velocity is constant,
the triangulation is regular, and the edges of triangulation are bounded away
from the direction of flow.
A challenge for the RFB method is to determine the residual free
bubble functions in a higher dimensional situation. The twolevel finite element
method (TLFEM) is a general framework to resolve this task. This method
2
was first successfully used for solving the Helmholtz equation [6]. In [7], the
numerical results for the advectiondiffusion problem showed that the TLFEM
performed as well as GLS and that there is no major qualitative difference
between these two methods. The TLFEM consists in partitioning the mesh into
submeshes, and an appropriate numerical method is then used to approximate
the PDEs governing the residualfree bubble basis functions instead of solving
analytically for the residualfree bubble functions. We can partition each mesh
into different submeshes arbitrarily. For simplicity of implementation, uniform
submeshes are considered first. However, due to the choice of bubble functions,
there might be boundary layers near outflow boundaries for the bubble shape
s
functions. Therefore, it seems more appropriate to choose a new submesh which
is able to capture the layer to produce more accurate approximate residualfree
bubble shape functions.
The purpose of this thesis is to introduce the new nonuniform sub
meshes depending on the direction of the flow, and compare numerical solutions
of the TLFEM for advectiondiffusion with uniform submeshes to see if we can
get any improvement on the numerical solutions under the same global mesh.
This thesis is organized as follows: In the next chapter, we review the
residualfree bubble method for the general boundary value problem. Then,
3
in chapter 3, we discuss the twolevel finite element method for the advective
diffusive equation derived from the RFB method ; In chapter 4, the algorithm
for generating nonuniform submesh is developed. Finally, numerical results
and conclusions are presented in chapter 5 and chapter 6 respectively.
4
2. The review of RFB for boundary
value problems
Let us consider Q to be an open bounded domain in R2 with boundary
For simplicity, we assume that dQ, is a polynomial curve in which case we
say that is a polynomial domain. (If 9Q is a curve, we can approximate it
with a polynomial).
First, we consider the general boundary value problem
Lu = f
in Q
(2.1)
I u = 0 on dVt
where L is a linear differential operator, e.g. L may be the advectivediffusive
operator, u is the unknown scalar function and / is a given source function.
We also assume that this problem is well posed.
The abstract variational formulation of (2.1) is as follows: Find a
scalar variable it V' such that
a(u,v) = {f,v) V v V, (2.2)
where V is a Hilbert space, a(, ) is a bilinear form from V x V to R and (, )
is the usual scalar product in L2(Q).
To specify the standard Galerkin finite element for (2.1), we partition the
5
domain fl into several pieces K (e.g. triangles, quadrilaterals etc.) in the stan
dard way, which forbids overlapping or any vertex on the edge of a neighboring
element, and so on. Thus,
U K = KlUK2U...uKm (2.3)
KeTh
where Th is a partition of f2.
We introduce the mesh parameter
h = max diam(i'f), diam(iir)= diameter of K (2.4)
KeTh
We now define 14 as a finitedimensional space, which is a subspace of V. Then
the standard Galerkin finite element method is: Find uh Â£ V* such that
a(uh,vh) = (f,vh) VvhÂ£Vh (2.5)
Now, we decompose the space 14 such that 14=14 + B, where 14 is the space
of continuous piecewise linear or bilinear polynomials and B is the space of
residualfree bubbles We will define the space B explicitly later. Then every
Uh Â£ Vh can be written in the form of Uh U\ + Ub, where ui Â£ 14 and
Ub Â£ B. For the residualfree bubble space, we require the bubble component
Ub of each vh to vanish on dK of each K and require each vh to satisfy the
original differential equations strongly, i.e.
L(i + Ub) = f in K (26)
6
or Lub = {Lui f)
in K
(2.7)
subject to zero Dirichlet boundary condition on the element boundary, i.e.
ub 0 on dK. (2.8)
By the classical static condensation procedure, first we set Vh = vbtK in K and
Vh = 0 elsewhere in (2.5) to have
a(ui + vbiK)K = (/, vbtfc)x VvbjK B (2.9)
where a(,)x and (,)k indicate that integration is restricted to the element
K.
Then, taking Vh = V\ in (2.5), we obtain
a(wfc,vi) = C/>i) (2.10)
(iq +Ub,Vi) = (j>i) (2.11)
The formulation (2.9) is automatically true due to our choice of bubble func
tions. In fact, this equation is also the variational formulation of (2.6) using
vb as test function restricted to each element K. Furthermore, the equation
(2.11) is the method to compute an improved bilinear or linear approximation
due to the residualfree bubbles effect. To find the residualfree bubble part
of the solution, we need to solve (2.5) which depends on the linear part of the
7
solution u\. Instead, bubble shape functions with i varying from one to the
number of element nodes (Nen) can be obtained by the following auxiliary
problems:
(i) For each i=1, 2, .. Nen, find fax such that
L(f>i,K = Llpi'K in K (2.12)
<;Kk = 0 on dK (2.13)
where the tpiiK are local basis function for U\ and
(ii) find
L4>f,K = / in K (2.14)
Vlk = 0 on dK (2.15)
Thus if Nen 1 = H ci,K<, i=1 h,K (2.16)
then
ub Nen ci,Ki>i,K + i>f,K (2.17)
i= 1
with the same coefficient Cij<.
Furthermore, we have the following representation:
Nen
Uh,K X] ci,K{4>i,K +
i=1
(2.18)
8
Hence, we can define the residualfree bubble space Bk as
Bk = span{pK, : faenjc, f} (2.19)
and
B= Â£ Bk (2.20)
Kerh
Therefore, we can restate the RFB method for general boundary value problem
as the following:
/
Find Uh = U\ + U(, E. Vh = Vi + B such that
' a{uh,vi) = (/, vi) Vuj e Hi (221)
a{'U'h,Vb,K)K = (/) ViF G Th and v^k Bk
\
or we can eliminate Ub.K in (2.21) and obtain an equation which only involves
up.
*
Find ui Vi such that
< (2.22)
a(i, ui) + vi)k = (/, Ui) V! Vi
To get Ub.K as a function of u\ and / requires solving equations (2.12)(2.15),
which is as complicated as solving the original differential equation, unless we
have special cases such as rectangular elements such that we can employ clas
sical analytical tools to get an exact solution within each element. Henceforth
our strategy is to approximate the bubble shape functions ipijK and by an
other appropriate finite element method. In other words, at the global level we
9
use the standard Galerkin method with piecewise linears for the original prob
lem. At the element level, we partition each element into a finer submesh and
then utilize nonstandard finite element methods to solve the bubble problems.
This is why we called this method a twolevel finite element method (TLFEM).
In this work, instead of GLS, we apply the improved Unusual Stabilized Finite
Element Method (the improved USFEM) [8] to approximate the bubble shape
functions. We discuss TLFEM in more detail in the next section.
10
3. The TLFEM for the advectivediffusive problem
In the present section we develop the TLFEM based on the RFB
method and discuss its application to the advectivediffusive problem. This
approach is general, and can be applied to an arbitrarily shaped domain with
out any difficulty. For the advectivediffusive problem, we set
L = kA + a V (3.1)
Here a is the given velocity field, assumed to be constant in each element, and
k is the given positive constant diffusivity coefficient. We are interested in the
advectivedominated case, i.e. when k C a The associated bilinear form is
a(u, v) = k(Vu, Vv) + (a V, v) (3.2)
At the global level, we rewrite equation (2.11) for the advectivediffusive prob
lem:
k(Vui, Vi>i) + K(Vub, Vvi) + (a Vui.vi) f (a Vub, tq) = (/, Vi) (3.3)
11
Then, let us consider the second term above:
K(Vit6,Vvi) = Vui)*
K
= ^2 k J Vub Viq dx
= WftVux nds k / ubAvi dx),
JdK JK
where n is an unit normal vector. The last equality is obtained by integration
byparts. Since the residualfree bubble function ub s are zero on element
boundaries and v\ is bilinear inside rectangular elements, we conclude that
K(Vub, Vvi) = 0.
Therefore, equation (3.3) simplifies to
k(Vui, Vui) + (aWui, iq) + (a Vub, tq) = (/, tq) (3.4)
Substituting (2.18) into (3.4) and setting iq = ipj enables us to write the ma
trix problem: Find the coefficient Cjs such that
^ci[(/cV'i/ji, V%) + (a V'ibi,'ipj) + (a Vfc.ipj)] = (/,^) (aV^/,^) (3.5)
i runs over all unknown interior nodes in the elements, say through N.
12
At the local element level, (2.12)(2.15) for the advectivediffusive problem
reads:
(i) for i = 1,  , Nen,
(ii)
a V4>i>K K,A0ijK = aVip^K in K (3.6)
4>lk = o on dK (3.7)
a K.A
= 0 on dK (3.9)
For each element K, let dK~ = { x e dK : a n(x) < 0} be its
inflow boundary and dK+ = { x G dK : an(x) > 0} be its outflow boundary,
where n is the outward normal unit vector to dK. Assume that a n(x) is
bounded away from zero; then we have dK = dK~ U dK+.
For small diffusivity we wish to consider the corresponding reduced problems
obtained by setting k =0:
(i) for i = 1, , Nen,
aVq>i,K = in K
on dK'
(3.10)
(3.11)
13
a V cpf>K f
in K
(3.12)
4>},k o
on dK
(3.13)
The characteristics of the reduced problems (3.10) and (3.12) are straight lines
parallel to the velocity field a.
Let us state some basic facts concerning the exact solutions of (3.6), (3.7), (3.8)
and (3.9):
(i) If we consider that the inflow boundary data of the reduced problem is a
discontinuous function g instead of the zero function, then the solution of the
reduced problem may be discontinuous with a jump across the characteristic.
In the full problem, the solution is continuous in K and the jump will be spread
out in a small region around the characteristic.
(ii) If the values attained by the reduced problems on dK+ do not coincide
with the boundary values described in the full problems, then solutions of the
latter problems will have a boundary layer at dK+: i.e. at a very narrow region
near the outflow boundary of the element, where the solution and its derivative
change abruptly.
For a given unit vector v with components (t^,^), the directional derivative
14
of u is defined by
vVu = Dvu (314)
In other word, Dvu is the component of Vit is the direction of v Hence,
a y can be viewed as the directional derivative of 4>^k in the streamline
direction multiplied by a similarly for a and a Therefore,
the analytical solution of (3.10) and (3.11) can be obtained simply by taking
integrals from x~ to cc at both sides of the equations along the streamline di
rection and then applying the fundamental theorem of calculus. The solution
of (3.10) and (3.11) is given by:
(3.15)
and similarly the solution of (3.12) and (3.13) is:
(3,16)
where x is the intersection of a line parallel to a passing through the point
x in (3.15) with the inflow boundary dK~. In addition, we define x+ as the
15
coordinate of the point on the outflow boundary dK+, aligned with the point
x+ in the streamline direction. (See Figure 3.1)
Figure 3.1. Obtaining x and x+ from any x K
To get the general idea of the behavior of the residualfree bubble shape func
tion, let us consider a simple example:
Let K be an unit square element defined in [0,1] x [0,1] with an uniform veloc
ity field of size one forming a 45 with the horizontal axis. (See Figure 3.2)
We assume the nodal points are labeled in ascending order in the counterclock
wise direction and let Side(l) be the edge joining the nodal points 1 and 2, the
Side(2) is the edge joining the nodal point 2 and 3 and so on. (See Figue 3.3).
16
Figure 3.2. Problem statement
Figure 3.3. Nodal points and sides numbering
17
It can be easily seen that Side(l) and Side(4) are inflow boundaries and Side(2)
and Side(3) are outflow boundaries for this example. Then, we can define the
local shape function tpi{xys explicitly.
$i{x,y) = (iz)(i y)
ifafav) = (i y)
y) = xy
= (i ~x)y
By formula (3.15), the values of each bubble shape function on the outflow
boundary dK+ can be evaluated by:
Mx.v) = < 1 X : on Side(2)
i y K : on Side(3)
o3{x,y) = < y  on Side(2)
x : on Side(3)
*(x,y) = 0
For bubble functions
18
boundary values prescribed in the full problems. It turns out that we are
expecting there may be an outflow boundary layer for these bubble shape
functions. This example gives us a general idea of how to choose the optimal
submesh (a mesh defined for each element) which is able to capture a jump
discontinuity of the exact solution in a thin numerical layer. Our strategy is to
use the nonuniform submesh which is more refined near the outflow boundaries
depending on the direction of the velocity field. Figure 3.4 illustrates this idea.
Figure 3.4. Submesh refines on outflow boundaries
We begin to approximate residualfree bubble shape functions by partitioning
each element K into the coarse submesh K* (a mesh defined for each element)
with diameter h* and denote by ipf the basis function for a piecewise linear
interpolation on the submesh.
19
Therefore, our unknown bubble basis function can be approximated by
^ = Â£4% O17)
i
Here N* is the number of all unknown interior node in the submesh. The index
i refers to specific bubble function that we are trying to compute.
Similarly, under the presence of a source term /, the other bubble shape func
tion is given by:
i
Then we can formulate the improved USFEM [8] for the case in which the
reactive term is equal to zero in the matrix formulation for (3.6) and (3.7) as:
for each % (from 1 to Nen) find cf \ such that
K W,*,0 + + (O v#,T0 Ol
l
= (a V*,*, > + TO. VO (3.19)
In addition, the matrix formulation for (3.8)(3.9) is given by: find c{, such
that
E 4 [( v**. o+w;. vo+( . to oi = u,r+ v
i
(3.20)
We use the stability r suggested in [8] as
r(x,PeK(x)) =
h*K2
6k + 6K^(PeK(x))
(3.21)
20
PeK{x)
C(z)
\a{x)\2h*K
3 k(x)
/
1 , if 0 < x < 1,
<
x if x > 1,
(3.22)
(3.23)
M*)2 = BM*)2!)1'2. (3.24)
Z=1
And the element parameter h*K is computed by using the largest streamline
distance of elements. See Figure 3.5.
Figure 3.5. The element parameter computation
Once the constants and c{ are determined, we substitute them in (3.17)
and (3.18) respectively to get the approximate residual basis function and
For a practical problem such as the flow over an airfoil or an automobile, the
direction of the flow varies elementwise. Therefore, it is necessary for us to
design a subroutine in our computer program to generate the nonuniform sub
mesh automatically. Before proceeding to discuss the algorithm for a submesh
21
generator, let us list all combinations of outflow and inflow boundary segments
of elements and indicate the submeshes which are used for each different case.
Other permutations are not listed, since they are covered by all combinations
listed in Figure 3.6.
22
I
I
II
o
I
o
o o
o
I
o
I
o
o
o o
o
o
Figure 3.6. Submesh used for different cases: I indicates inflow boundary
and O indicates outflow boundary.
23
4. The algorithm for nonuniform submesh generator
We begin our study of this subject with the algorithm for very sim
ple meshes in one dimension and proceed to generalize to more sophisticated
nonuniform irregular mesh generation schemes for the two dimensional case.
An algorithm for an uniform mesh on interval (a, b) is easily constructed:
STEP 1 Set h = n=number of elements
STEP 2 For i = 1,  , n:
Set oci = Then & = a+(ba)gi(di) defines the node location,
where gi(a) = a.
For our nonuniform mesh, we can define two mappings:
g2{oi) = a2 (4,1)
gz(a) = (1 a)2. (4.2)
Replacing g\{a) by g2(a) in the STEP 2, the mapping Â£2(q0 will take
the uniform mesh in a to a quadratically graded mesh from left to right in Â£.
(See figure 4.1) Conversely, the mapping g3 (a) will take the uniform mesh in
a to a quadratically graded mesh from right to left in Â£, if we rewrite STEP
2 as:
24
For i = 1, , n:
Set cti = (i 1 )h. Then Â£* = b (b a)g3(tti) defined the node
location.
Or if we wish to produce a mesh that is graded into both ends Â£ = a and Â£ = b,
we can combine these two mappings (/2(a) and gs(a) and define
94{a) = <
(1 a)2
0 < a < 0.5
0.5 < a < 1
(4.3)
Figure 4.1. Mesh grading from left to right: a quadratic transformation
Now. let us generalize this idea to a two dimensional irregular mesh.
Let If be a quadrilateral element defined by locations of its four nodal points
25
x%, a = 1 , 4 in the physical domain R2 and f2e = [0,1] x [0. l] be a corre
sponding biunit square in the computational domain. The nodal points and
sides of elements either in the physical or the computational domain are labeled
in ascending counterclockwise direction. See Figure 4.2.
The coordinates of a point (Â£,77) in the biunit square Qe are related to the
coordinates of a point (x, y) in K by the transformations:
4
x(^v) = I]Nq(Â£,77K a=l (4.4)
y(,v) = 'Z/Na{^y)yea, a=l (4.5)
where ^(Â£,77) = ^(1 Â£)(1 77), ^(Â£,77) = ^(1 + Â£)(1 77),
^3(^77) = J(1+ Â£)(! +77), /V4(Â£, t?) = ^(1Â£)(1+ 77).
Let us introduce some notations used in the algorithm.
1. Nsd: number of space dimensions (Here Nsd=2),
2. N*es: number of elements for submesh (N*e.s = n x m),
3. N*np: number of nodal points for submesh.
4. a: local node number for mesh (1 < a < 4)
5. a*: local node number for submesh (1 < a* < 4).
6. A* global node number for submesh (1 < 4* < N*np).
26
Figure 4.2. Local node and side ordering
7. i: spatial index (1 < i < Nsd),
8. e* : element number for submesh (1 < i* < N*es),
9. ID*(A*): destination matrix for submesh.
10. IEN*(a*,e*): location matrix for submesh.
11. xK (a, i): coordinates of nodal points for mesh
12. xK (A*, i) coordinates of global node number for submesh
Similarly to the global level, three data processing arrays, IEN*, ID* and
xK'(A*,i*) are needed as input data. For more detailed discussion see [9].
Here is the algorithm for the irregular mesh generator:
INPUT: xK(a, i), m,n
OUTPUT: ID*(A*), LM*(a*,e*), xK'{A*,i)
27
STEP 1 Determine whether each boundary of element K is inflow or out
flow by the definition and then properly label each corresponding segment of
the biunit square Â£}e.
STEP 2 Set h = , k = .
STEP 3 Compute the horizontal coordinates of the submesh for the biunit
square fie:
Define the following mapping:
/l(z) = z3 (4.6)
fi(x) = (l^)3 (4.7)
= X (4.8)
Note: The domains of these three functions are [0,1].
Case 1: Side(2) and Side(4) both are inflow boundary segments
For i = 1, ,n + 1:
Set ai = (i 1 )h. Then Â£ = 1 + 2/3(0;*);
Case 2: Side(2) and Side(4) both are outflow boundary segments
For i = 1, , [n/2]:
Set a.i = (i 1 )h. Then f* = 1 + 2/1(oj) and
for i = [n/2 + 1], , n + 1:
28
Set a.i = (i 1 )h. Then Â£ = 1 2/2(aj);
Case 3: Side(2) is inflow and Side(4) is outflow boundary segment
For i = 1,, n:
Set ai = (i 1 )h. Then ^ = 1 + 2
Case 4: Side(2) is outflow and side(4) is inflow boundary segment
For i = 1, , n:
Set ai = (i 1 )h. Then ^ = 1 2/2(ofz);
STEP 4 Compute the vertical coordinates of the submesh for biunit square
fte.
Similar to STEP 3, now we consider Side(l) and Side(3) of Qe.
STEP 5 Set ID*(j4*)=0 on the boundary nodes and eg=l;
STEP 6 For i = 1, , n; j = 1,..., m: do STEP 7 STEP 10
STEP 7 Construct ID*(A*) for submesh
Set A* = (i 1) (n + 1) + j, A* =
If ID*(A*) 7^ 0 then ID(A*)=eg and set eq=eq+1;
STEP 8 Construct IEN*(a*,e*) for submesh:
Set e* = (i 1) (n + 1) + j;
Then IEN*(l,e*)=A*, IEN*(2,e*)=A*+l,
IEN*(3,e*)=A* + n+2, IEN*(4,e*)=A* + n+1.
29
STEP 9 Compute the coordinate of node points (node .A*):
61 = Â£
Va = Vj (49)
STEP 10 Map the coordinates of node points {^aiVa) from the computa
tional domain to the physical domain by the transformations (4.4) and (4.5).
30
5. Numerical results
In this chapter, we will report three series of experiments for the
advectiondiffusion problem with TLFEM using the new nonuniform submesh
introduced in the previous section. In the following numerical results, we il
lustrate the applicability of the method for singularly perturbed problems, i.e.
for small values of the diffusivitv k compared with the advection field .
5.1 A problem with discontinuous
boundary condition
In the first two experiments, let us consider a unit square domain
Q with discontinuous boundary values at (0.5,0) and (0,1). The diffusivity is
k = 106, and we will test our methods for two different angles of the uniform
velocity field of size one with the horizontal axis. The first case is 45 and the
second one is 60. (see Figure 5.1 for problem statements).
In these problems, a discontinuous data at the inflow boundary is propagat
ed into the domain which causes an internal layer along the characteristic of
the problem starting at point(0.5,0). In addition, the problems are subjected
to homogeneous Dirichlet boundary conditions at the outflow boundary which
create the outflow boundaries.
31
0 u= ^ 0.5 u=o 1 x
Figure 5.1. Domain for skew convection problems with two different velocity
fields
32
For both cases, we employ an 20 x 20 uniform mesh for the TLFEM and the
GLS methods. First, we compare the different submesh partition strategies
for the TLFEM, say 10x10 uniform submesh and nonuniform submesh( See
Figure 5.2) For bubble shape function analysis, we pick one element to plot
Figure 5.2. Two different 10x10 submeshes used for TLFEM: uniform (left)
and nonuniform (right).
four shape bubble functions. The comparison of bubble shape functions by us
ing uniform and nonuniform submesh for 45 case is shown in Figure 5.3, and
for 60 is shown in Figure 5.4. For both cases, it is obvious that nonuniform
submeshes are successful to capture the outflow boundary layers, and produce
more accurate residual free shape functions. Meanwhile, the nonuniform sub
mesh yields slightly better performance in the global solutions for both cases,
(see Figure 5.5 and 5.6.) In Figure 5.7 and Figure 5.8, both the TLFEM and
the GLS perform similarly except in the crosswind internal boundary layer.
TLFEM performs better than GLS therein.
33
'1
'1
Figure 5.3. Comparison of bubble shape functions by using two different
submeshes for the 45 problem: uniform (left) and nonuniform (right)
34
'1
l
Figure 5.4. Comparison of bubble shape functions by using two different
submeshes for the 60 problem: uniform (left) and nonuniform (right)
35
uniform submesh
nonuniform submesh
Figure 5.5. Comparison of TLFEM solutions by using two different submeshes
for the 45 problem: uniform (left) and nonuniform (right)
uniform submesh
nonuniform submesh
1.5 r
1 1 Xaxis
Yaxis
Figure 5.6. Comparison of TLFEM solutions by using two different submeshes
for the 60 problem: uniform (left) and nonuniform (right)
36
GLS
TLFEM
Figure 5.7. Comparison of GLS and TLFEM solutions for the 45 problem
GLS (left) and TLFEM (right)
GLS TLFEM
Figure 5.8. Comparison of GLS and TLFEM solutions for the 60 problem
GLS (left) and TLFEM (right)
37
5.2 Thermal boundary layer problem
Let us consider a rectangular domain of sides 1.0 and 0.5. subject to
the boundary conditions presented in Figure 5.9. The velocity field is given by
a = (2y,0) and diffusivity k = 7 x 104. This problem can be viewed as the
simulation of the development of a thermal boundary layer on a fully devel
oped flow between two parallel plates, where the top plate is moving with the
velocity equal to one and the bottom plate is fixed.
Figure 5.9. Problem statement for the thermal boundary layer problem
The nonhomogeneous mesh for GLS and TLFEM methods consists 21 equally
spaced nodes in the xdirection, 11 nodes uniformly distributed in the interval
[0,0.1] and 11 nodes equally spaced on [0.1,0.5] in the ydirection. (see Figure
5.10)
We wish to compare TLFEM by using the different submeshes, say 10 x 10
uniform and nonuniform submesh (see Figure 5.11).
38
0.5
0.45
0.4
0.35
0.3
0.25
>
0.2
0.15
0.1
0.05
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Xaxis
Figure 5.10. Mesh used for the thermal boundary layer problem
Figure 5.11. Submesh used for the thermal boundary layer problem
39
Because the flow velocity increases along the ydirection, the characteristics of
solutions in the domain can be divided into two regions: for the bottom of the
domain, we have a diffusive dominated case ; for the top of the domain, we
have the advective dominated case.
We find out that the uniform submesh performs better than the nonuniform
mesh for the diffusive dominated case. On the other hand, the nonuniform
submesh still does a good job for the advective dominated case (See Figures
5.13 and 5.14). The global solutions for the two different submeshes are almost
the same (See Figure 5.12).
uniform submesh
nonuniform submesh
Xaxis
1 0
Yaxis
Xaxis
1 0
Yaxis
Figure 5.12. Comparison of TLFEM solutions using two different submeshes
40
1
*1
0.5
0.5
0.5
0.5
Figure 5.13. Bubble shape functions for the diffusive dominated case: uniform
submesh (left) and nonuniform submesh (right)
41
*1
0.5
0.5
Figure 5.14. Bubble shape functions for the advective dominated case:
form submesh (left) and nonuniform submesh (right)
uni
42
Next we change to a smaller diffusivity k = 106 and rerun TLFEM by using
the same meshes and submeshes and compare it with the GLS method (See
Figures 5.15 and 5.16 for the bubble shape functions comparison). Taking the
top plate velocity as the characteristic flow velocity, we have the meshPeclet
number Pe =
\a\h
2k
= 25000. This means that the entire domain is advec
tively dominated. The numerical results show that the nonuniform submesh
performs slightly better than the uniform one (See Figure 5.17). Meanwhile,
the TLFEM performs better than the GLS (See Figure 5.18). There are some
oscillations near the outflow boundary.
43
'1
'1
Figure 5.15. Bubble shape functions approximated by two different submesh
es for the bottom region: uniform (left) and nonuniform (right)
44
'1
1
Figure 5.16. Bubble shape functions approximated by two different submesh
es for the top region: uniform (left) and nonuniform (right)
45
uniform submesh
nonuniform submesh
Figure 5.17. Comparison of the TLFEM solutions by using two different
submeshes for the thermal boundary layer problem when k = 10~6: uniform
(left) and nonuniform (right)
GLS
TLFEM
Xaxis
Yaxis
0.5
Figure 5.18. Comparison between the GLS method and the TLFEM method
for the thermal boundary layer problem when k = 10~6: GLS (left) and
TLFEM (right)
46
5.3 Bubble ramp problem
Now. let us consider an Lshaped domain with external source /=1,
k = 10~6, and homogeneous Dirichlet conditions depicted in Figure 5.16
y
1
ii=0
ii
u=0 _______ Velocity:
a, = 1.0 a. = 0.0
0.5 M=0
: M=0
i u=0 /= l.O
0.........0.5 . 15 x
u = 0
Figure 5.19. Problem statement for the bubble ramp problem
A uniform partition into 300 elements are employed for the global mesh and
two kinds of 10x10 submeshes are used for TLFEM. Because of the presence of
the external force, the additional residual free bubble shape function
to be determined. The comparison of these approximated solutions using two
different submeshes are shown in Figure 5.20, where the nonuniform submesh
performs better than the uniform one.
ii=0
Velocitv:
a o o li S3' O II
K=0
u=0 /= io
47
uniform
nonuniform
0.1
0.1
0
0.05
0.05
Yaxis 0 5
Xaxis
Yaxis
0 5
Xaxis
Figure 5.20. Comparison of bubble shape function f by using two different
submesh for the bubble ramp problem: uniform (left) and nonuniform (right)
The solution exhibits a strong outflow boundary layer along x=1.5, two cross
wind boundary layers along y=l and y=0, and a crosswind internal layer along
y=0.5. Thus, it is one of the most stringent tests for the advective diffusive
problem. We find out that there are some improvements near outflow bound
aries in the numerical solutions for the TLFEM due to the new submesh we
choose. Meanwhile, the numerical results of the GLS and the TLFEM using
nonuniform submeshes are almost indistinguishable (See Figure 5.21).
48
TLFEM(nonuniform) TLFEM(uniform)
Yaxis 0 o xaxis
Figure 5.21. Comparison of the GLS and the TLFEM by using two different
submeshes for the bubble ramp problem
49
6. Conclusion
In this work, we formulated the twolevel finite element method based
on the standard Galerkin method using piecewise linears in the original mesh,
and the usual stabilized finite element method was used to approximate the
partial differential equations governing the residual free bubble functions. Once
these residual free bubble functions are determined, we can substitute them
into the Galerkin formulation to improve the accuracy of the global numerical
solution. The main advantage of this method is that we do not have to solve
these partial differential equations analytically. Therefore, it is suitable for the
finite element computation in a practical problem.
Due to the choice of bubble functions, the boundary layers occur
at outflow boundaries for the advectivedominated case. The new submesh
strategy w^as introduced to cure this problem. The numerical experiments
confirm our idea that nonuniform submeshes are able to capture the layer of
residual free bubble shape functions. The result display some improvements
in the solution of the advectivediffusive problem, and the TLFEM solutions
perform better than the GLS method in some cases. We note that the idea of
the new submesh strategy can be applied to other problems in fluid dynamics,
50
such as the incompressible NaiverStokes equations. Nonuniform submesh
depending on the direction of the flow can approximate the solution of the
pressure term more accurately. See [12] for its application.
51
References
[1] F. Brezzi, M.O. Bristeau, L.P. Franca, M. Mallet, and G. Roge. A relation
ship between stabilized finite element methods and the Galerkin method
with bubble functions. Comput. Methods Appl. Mech. Engrg., 96: 117
129, 1992.
[2] F. Brezzi, L.P. Franca, A. Russo. Further considerations on residualfree
bubbles for advectivediffusive equations. Methods Appl. Mech. Engrg.,
166: 2533, 1998.
[3] A.N. Brooks and T.J.R. Hughes. Streamline upwind/PetrovGalerkin for
mulation for convection dominated flows with particular emphasis on the
imcompressible NavierStokes equations. Comput. Methods. Appl. Mech.
Engrg., 32: 199259, 1982.
[4] G.F. Carey and J.T. Oden. Finite Elements: Computational aspects.
PrenticeHall, 1984.
[5] L.P. Franca, S.L Frey and T.J.R. Hughes. Stabilized finite element meth
ods: I. Application to advectivediffusive model. Comput. Methods. Appl.
Engrg., 95: 253276, 1992.
[6] L.P. Franca and A.P. Macedo, A twolevel finite element method and its
application to Helmholtz equation. Int. J. Numer. Methods Engrg., 43:
2332, 1998.
[7] L. P. Franca, A. Nesliturk, and M. Stynes. On the satibility of residualfree
bubbles for convectiondiffusion problems and their approximation by a
twolevel finite element method. Comput. Methods Appl. Mech. Engrg.,
166: 3549, 1998.
[8] L. P. Franca and F. Valentin. On an improved ususual stabilized finite
element method for the advectivereactivediffusive equation, To appear
[9] T.J.R. Hughes. The Finite Element Method: Linear Static and Dynamic
Finite Element Analysis. PrenticeHall, 1987.
[10] T.J.R. Hughes, L.P. Franca, and G.M. Hulbert, A new finite element
52
formulation for computational fluid dynamics: VIII. The Galerkin/least
squares method for advectivediffusive equations Comput.Methods Appl.
Mech. Engrg., 73: 173189, 1989.
[11] C. Johnson. Numerical Solution of Partial Differential Equations by the
Finite Element Method. Cambridge, 1987.
[12] A. Nesliturk. Approximating the incompressible Navier Stokes equations
by the twolevel finite element method. Ph.D thesis, University of Colorado
at Denver, 1999.
53
