On the Design of an Arti cial Di usion Model for the Lagrange-Galerkin Method on Unstructured ...

30 Pages·1996·1.26 MB·English

Report no. 96/07
On the Design of an Arti(cid:12)cial Di(cid:11)usion Model
for the Lagrange-Galerkin Method on
Unstructured Triangular Grids
Paul Houston Endre Su(cid:127)li
In this paper we propose the addition of an arti(cid:12)cial di(cid:11)usion
model to the Lagrange-Galerkin (cid:12)nite element method, which is de-
pendent on both the mesh function h and the numerical solution
uh. The purpose of this additional term is to produce smooth ap-
proximations to sharp features of the solution, such as internal and
boundary layers. Moreover, the added di(cid:11)usion will help to stabilise
the numerical scheme.
Further, we extend the a posteriori error analysis presented in [12]
to include this arti(cid:12)cial di(cid:11)usion model. Based on this a posteriori
estimate, we design an adaptive algorithm to ensure global control of
2 2
the error in the L (L ) norm with respect to a pre-determined toler-
ance. The performance of this numerical algorithm is demonstrated
by some numerical experiments.
Key words and phrases: Arti(cid:12)cial di(cid:11)usion, Lagrange-Galerkin (cid:12)nite element
methods, a posteriori error analysis, adaptive algorithms
The (cid:12)rst author would like to acknowledge the (cid:12)nancial support of the EPSRC.
The work reported here forms part of the research programme of the
Oxford{Reading Institute for Computational Fluid Dynamics.
Oxford University Computing Laboratory
Numerical Analysis Group
Wolfson Building
Parks Road
Oxford, England OX1 3QD May, 1996
2
Contents
1 Introduction 3
2 Notation and basic de(cid:12)nitions 4
3 Model problem and discretisation 5
3.1 The Lagrange-Galerkin method : : : : : : : : : : : : : : : : : : : 5
3.2 Design of the arti(cid:12)cial di(cid:11)usion model : : : : : : : : : : : : : : : 8
4 A posteriori error analysis 9
4.1 Error representation : : : : : : : : : : : : : : : : : : : : : : : : : 9
4.2 Interpolation/projection estimates for the dual problem : : : : : : 11
4.3 Strong stability of the dual problem : : : : : : : : : : : : : : : : : 15
4.4 Completion of the proof of the a posteriori error estimate : : : : : 16
5 Adaptive algorithm 20
6 Numerical Experiments 22
6.1 Experiment 1: Boundary and internal layer problem : : : : : : : : 22
6.2 Experiment 2: Rotating cylinder problem : : : : : : : : : : : : : : 25
7 Conclusions 28
References 29
3
1 Introduction
Over the last decade the Lagrange-Galerkin (cid:12)nite element method has been suc-
cessfully used for the numerical solution of unsteady convection-dominated dif-
fusion problems, see Bercovier et al. [2,3], Douglas & Russell [6], Lesaint [15],
Pironneau [17] and Su(cid:127)li [21,22], for example. This numerical scheme is based on
combininga(cid:12)niteelementprocedurewiththemethodofcharacteristicstoexploit
the `nearly-hyperbolic' nature of these problems. While the Lagrange-Galerkin
method gives a dramatic improvement over the standard Galerkin method in
terms of stability, the scheme may still produce oscillatory approximations in
the vicinity of sharp transitions in the numerical solution, such as internal or
boundary layers. Further, the integrals that arise in the formulation of the
Lagrange-Galerkin method involve the product of (cid:12)nite element basis functions
with their shifted counterparts. In general, these integrals cannot be evaluated
exactly and some form of approximation must be made. It is well known that
if quadrature formulae are employed, then the method (which is unconditionally
stable when exact integration is performed) becomes only conditionally stable,
see Morton et al. [16] and Priestley [18], for example.
To reduce these problems, Priestley [19,20] introduced the positive Lagrange-
Galerkinmethodbasedonthe(cid:13)uxcorrectedtransportalgorithmofBoris&Book
[4]. In this paper we shall consider a di(cid:11)erent approach: we propose the addition
of an arti(cid:12)cial di(cid:11)usion model to the basic Lagrange-Galerkin method dependent
on both the mesh function h and the numerical solution uh. This will not only
damp spurious oscillations generated by the method, but the added di(cid:11)usion will
help to stabilise the numerical scheme, cf. Su(cid:127)li [23]. Additionally, in the context
of mesh adaptation, we noted in [12] that for a convection-dominated di(cid:11)usion
problem, the spatial mesh re(cid:12)nement algorithm is governed by the `hyperbolic
part' of the residual of the underlying partial di(cid:11)erential equation. Further, it
was shown in [12] that this hyperbolic residual will concentrate the spatial mesh
near sharp features of the solution that are not orthogonal to the (cid:13)ow direction.
Indeed, internal layers that are aligned with the (cid:13)ow direction were less re(cid:12)ned,
despite the error being locally large in these regions. Thus, the addition of the
arti(cid:12)cial di(cid:11)usion into the numerical scheme will ensure that the spatial adaptive
algorithm is no longer dominated by the hyperbolic part of the residual, and
consequently that all sharp features of the solution will be appropriately re(cid:12)ned.
The outline of the rest of this paper is as follows: in Section 2 we summarise
some of the notational conventions that we shall use. In Section 3 we state the
model problem to be considered and formulate the Lagrange-Galerkin method
for this problem. Furthermore, in this section we design an appropriate arti(cid:12)cial
di(cid:11)usion model and formulate the `discontinuity capturing' Lagrange-Galerkin
method. In Section 4 we extend the a posteriori error analysis presented in
[12] for the discontinuity capturing Lagrange-Galerkin method using the general
approach developed by C. Johnson and his co-workers (see [7,8,9,10,11,14], for
4
example). InSection 5wedesign anadaptive algorithmbased onthisaposteriori
2 2
estimate to ensure global control of the error in the L (L )-norm with respect to
a pre-determined tolerance. In Section 6 we present some numerical experiments
to illustrate the performance of our adaptive strategy. Finally, in Section 7 we
summarise the work presented in this paper.
2 Notation and basic de(cid:12)nitions
Let Z denote the set of integers, N the set of positive integers, N0 the set of
+
non-negative integers, R the set of real numbers and R the set of positive real
numbers.
d
Let ! be a bounded open subset of R , where d 2 N. For 1 (cid:20) p (cid:20) 1,
p
let L (!) denote the usual Lebesgue space of real-valued functions with norm
2 2
k (cid:1) kLp(!). For p = 2 and for u;v 2 L (!) we denote by ((cid:1);(cid:1))! the L (!) inner
product de(cid:12)ned as
Z
(u;v)! := u(x)v(x)dx:
!
For ! = (cid:10), we denote k(cid:1)kL2((cid:10)) by k(cid:1)k, and ((cid:1);(cid:1))! by ((cid:1);(cid:1)).
If (cid:11) = ((cid:11)1;:::;(cid:11)d), with each (cid:11)i 2 N; i = 1;:::;d, we write j(cid:11)j := (cid:11)1+:::+
(cid:11)d and (cid:11)! := (cid:11)1!(cid:11)2!:::(cid:11)d!.
(cid:11) (cid:11)1 (cid:11)d
Let D := D1 :::Dd and Dj = @=@xj for 1 (cid:20) j (cid:20) d. For m 2 N0, we
m
denote by C (!) the set of all continuous real-valued functions de(cid:12)ned on !
(cid:11) m
such that D u is continuous on ! for all j(cid:11)j (cid:20) m. C (!(cid:22)) will denote the set of
m (cid:11)
all u in C (!) such that D u can be extended from ! to a continuous function
on !(cid:22) for all j(cid:11)j (cid:20) m. In particular, when m = 0 we simply write C(!(cid:22)) instead of
0 m m
C (!(cid:22)). C0 (!(cid:22)) will denote the set of functions in C (!(cid:22)) with compact support
m;1 m
in !(cid:22). Also, C (!(cid:22)) will denote the set of functions in C (!(cid:22)) for which, for
(cid:11)
0 (cid:20) j(cid:11)j (cid:20) m, D u is Lipschitz-continuous in !(cid:22).
m;p
Further, for m 2 N0, let W (!) denote the classical Sobolev space endowed
with the norm k (cid:1) kWm;p(!) and the semi-norm j (cid:1) jWm;p(!) (cf. Adams [1]). For
m m;2 m
p = 2 we write H (!) for W (!). H0 (!) will denote the closure of the space
m
of in(cid:12)nitely smooth functions with compact support in ! in the norm of H (!).
m (cid:0)m
The dual space of H0 (!) will be denoted by H (!).
2
Let X be any of the spaces just de(cid:12)ned. Then X will denote the topological
product X (cid:2)X.
5
3 Model problem and discretisation
Given a (cid:12)nal time T > 0, we shall consider the following unsteady convec-
2 (cid:0)1 2
tion-di(cid:11)usion problem: given that f 2 L (0;T;H ((cid:10))) and u0 2 L ((cid:10)), (cid:12)nd
1 2
u(t) 2 H0((cid:10))\H ((cid:10)) such that:
ut +a(cid:1)ru(cid:0)(cid:15)(cid:1)u = f; x 2 (cid:10); t 2 I; (3.1a)
u((cid:1);t) = 0; x 2 @(cid:10); t 2 I(cid:22); (3.1b)
u(x;0) = u0(x); x 2 (cid:10); (3.1c)
2
where (cid:10) is a bounded convex polygonal domain in R with boundary @(cid:10) and
I = (0;T]. Further, we assume that the di(cid:11)usion coe(cid:14)cient (cid:15) > 0, the velocity
(cid:12)eld a 2 C(0;T;C01((cid:10)(cid:22))2) and that a is incompressible, i.e.
r(cid:1)a = 0 8x 2 (cid:10); t 2 I: (3.2)
3.1 The Lagrange-Galerkin method
In this section we shall formulate the standard Lagrange-Galerkin method for
(3.1) without the additional arti(cid:12)cial di(cid:11)usion term. However, we (cid:12)rst need to
introduce the following notation.
0 1 M M+1
Let 0 = t < t < ::: < t < t = T be a subdivision (not necessarily
n n(cid:0)1 n
uniform) of I, with corresponding time intervals I = (t ;t ] and timesteps
n n(cid:0)1 n
kn = t (cid:0)t . For each n, let T = fKg be a partition of (cid:10) into closed triangles
K, with corresponding mesh function hn 2 C0;1((cid:10)(cid:22)) satisfying
n
jrhn(x)j (cid:20) (cid:22) 8x 2 K 8K 2 T ; (3.3a)
2 n
c1hK (cid:20) meas(K) 8K 2 T ; (3.3b)
n
c2hK (cid:20) hn(x) (cid:20) hK 8x 2 K 8K 2 T ; (3.3c)
where hK is the diameter of K (= longest side of K) and c1, c2 and (cid:22) are
positive constants. Further, h is de(cid:12)ned to be the global mesh function de(cid:12)ned
n
by h(x;t) = hn(x), for (x;t) 2 (cid:10)(cid:2)I and we de(cid:12)ne the corresponding time step
n
function k = k(t) by k(t) = kn; t 2 I .
n n
Let S = (cid:10)(cid:2)I ; for p;q 2 N0 we de(cid:12)ne
hn 1 n
S = fv 2 H0((cid:10)) : vjK 2 Pp(K) 8K 2 T g;
q
X
hn j hn
V = fv : v(x;t)jSn = t vj; vj 2 S g;
j=0
h hn
V = fv : v(x;t)jSn 2 V ; n = 1;:::;M +1g;
where Pp(K) denotes the set of polynomials of degree at most p over K.
In the following, we shall assume that p = 1 and q = 0. We note that if
hn
v 2 V for n = 1;:::;M + 1, then v is continuous in space at any time, but
6
n
may be discontinuous in time at the discrete time levels t . To account for this,
we introduce the notation
n n
v(cid:6) := sl!im0+v(t (cid:6)s);
and
n n n
[v ] := v+ (cid:0)v(cid:0):
n n
For some n 2 N0, we associate with T the set E = f(cid:28)g consisting of those
2 n
line segments in R which appear as an edge of some K 2 T . We also denote
n n
by Ei , those (cid:28) in E which are interior to (cid:10) (i.e. not part of @(cid:10)).
n
For each (cid:28) 2 Ei , let n(cid:28) denote the unit normal to (cid:28) in the outward direction
hn
to K, and de(cid:12)ne for w 2 S (for some n 2 N0),
" #
@w
= lim (rw(x+sn(cid:28))(cid:0)rw(x(cid:0)sn(cid:28)))(cid:1)n(cid:28); x 2 (cid:28);
@n(cid:28) s!0+
that is, [@w=@n(cid:28)] is the jump across (cid:28) in the normal component of rw. Finally,
we introduce the discrete second derivatives
(cid:12)" #(cid:12)
2 X (cid:12)(cid:12) @w (cid:12)(cid:12) 1 n
DhwjK = (cid:28)2@K\Ein(cid:12)(cid:12) @n(cid:28) (cid:12)(cid:12) hK; K 2 T ; (3.4)
(cid:12)" #(cid:12)
(cid:14)2 X (cid:12)(cid:12) @w (cid:12)(cid:12) 1 n
DhwjK = (cid:28)2@K\En(cid:12)(cid:12) @n(cid:28) (cid:12)(cid:12) hK; K 2 T ; (3.5)
where, in (3.5) we de(cid:12)ne rw(x+sn(cid:28)) = 0 if (cid:28) 2 @(cid:10).
To construct the Lagrange-Galerkin method, we (cid:12)rst de(cid:12)ne the particle tra-
jectories (or characteristics) associated with problem (3.1): the path of a particle
located at position x 2 (cid:10) at time s 2 I is de(cid:12)ned as the solution of the initial
value problem
d
X(x;s;t) = a(X(x;s;t);t); (3.6a)
dt
X(x;s;s) = x: (3.6b)
1 (cid:22) 2
We note that for a 2 C(0;T;C0((cid:10)) ) there exists a unique solution to (3.6).
The Lagrange-Galerkin method makes use of the material derivative, Dtu,
which is de(cid:12)ned, for u smooth enough, as follows:
d
Dtu(x;t) := u(X(x;s;t);t) js=t
dt
@
= u(x;t)+a(x;t)(cid:1)ru(x;t) 8x 2 (cid:10); t 2 I: (3.7)
@t
7
Hence, using the material derivative, equation (3.1) may be rewritten in the
following (weak) form: (cid:12)nd u(t) 2 V, such that
(Dtu((cid:1);t);v)+((cid:15)ru((cid:1);t);rv) = (f((cid:1);t);v) 8v 2 V; (3.8a)
(u((cid:1);0);v) = (u0;v) 8v 2 V; (3.8b)
1
where V = H0((cid:10)). The Lagrange-Galerkin time-discretisation involves approx-
imating the material derivative by a divided di(cid:11)erence operator along the par-
ticle trajectories. The simplest appropriate discretisation is the backward Euler
method, giving for n = 0;:::;M:
!
n+1 n+1 n n
u((cid:1);t )(cid:0)u(X((cid:1);t ;t );t )
;v
kn+1
n+1 n+1
+((cid:15)ru((cid:1);t );rv) (cid:25) (f((cid:1);t );v) 8v 2 V; (3.9a)
(u((cid:1);0);v) = (u0((cid:1));v) 8v 2 V: (3.9b)
n n
If we now de(cid:12)ne uh to be the Galerkin (cid:12)nite element approximation to u((cid:1);t )
n
at time t ; then applying the (cid:12)nite element method to (3.9) yields the Lagrange-
n+1 hn+1
Galerkin discretisation of (3.1) as follows: (cid:12)nd uh 2 S for 0 (cid:20) n (cid:20) M such
that
!
n+1 n n+1 n
uh (cid:0)uh(X((cid:1);t ;t ))
;v
kn+1
n+1 (cid:22) hn+1
+((cid:15)ruh ;rv) = (f;v) 8v 2 S ; (3.10a)
0 h0
(uh;v) = (u0;v) 8v 2 S ; (3.10b)
where f(cid:22)jSn+1 := f((cid:1);tn+1). This is the same approach as that used by Bercovier
et al. [2,3], Douglas & Russell [6], Lesaint [15], Pironneau [17] and Su(cid:127)li [21,22].
n+1
Alternatively, if we integrate (3.10a) with respect to t over I , we obtain
the following equivalent formulation: (cid:12)nd uh such that, for n = 0;1;:::;M,
uhjSn+1 2 Vhn+1 and satis(cid:12)es
h (cid:22) hn+1
(Dtuh;v)n+1 +((cid:15)ruh;rv)n+1 = (f;v)n+1 8v 2 V ; (3.11a)
0 h0
(uh(cid:0);v) = (u0;v) 8v 2 V ; (3.11b)
where
h
DtuhjSn+1
n+1 n+1 n+1 n+1 n n
:= (uh(cid:0)(X(x;t ;t );t )(cid:0)uh(cid:0)(X(x;t ;t );t ))=kn+1:(3.12)
2 n+1 2
Further, for v;w 2 L (I ;L ((cid:10))), we have used
Z tn+1
(v;w)n+1 := (v;w)dt:
tn
8
3.2 Design of the arti(cid:12)cial di(cid:11)usion model
In this section we discuss the design of an appropriate arti(cid:12)cial di(cid:11)usion model
for the Lagrange-Galerkin method. One approach is to modify the test functions
v in (3.10), (3.11); for example, if we consider the steady convection di(cid:11)usion
problem
a(cid:1)ru(cid:0)(cid:15)(cid:1)u = f; x 2 (cid:10); (3.13)
then the streamline di(cid:11)usion method applied to (3.13) is a generalisation of the
standard Galerkin (cid:12)nite element method with the test functions v replaced by
v +(cid:14)a(cid:1)rv;
where
1
2
(cid:14) = min((cid:13)1h =(cid:15);(cid:13)2h=jaj); (3.14)
2
and (cid:13)1, (cid:13)2 are positive constants, cf. Hansbo & Johnson [11]. This modi(cid:12)cation
of the test functions corresponds to adding arti(cid:12)cial di(cid:11)usion in the streamwise
direction; Hughes et al. [13] proposed additional modi(cid:12)cations of the test func-
tions to introduce arti(cid:12)cial di(cid:11)usion in other directions.
A second alternative is to explicitly add a term of the form
((cid:15)^ruh;rv)
to the numerical scheme, where (cid:15)^ = F(uh;h); e.g. for the streamline di(cid:11)usion
method applied to the time-dependent problem (3.1), Hansbo & Johnson [11]
proposed the following de(cid:12)nition for F:
2 3=2
F(uh;h) = max((cid:13)3h R(uh)(cid:0)(cid:15);(cid:13)4h (cid:0)(cid:15);0); (3.15)
where
n+1 n
R(uh)jSn+1 := ja(cid:1)ruh (cid:0)fj+j[uh]=kn+1j
and (cid:13)3, (cid:13)4 are positive constants.
Here, we adopt a similar approach, however, we shall propose a di(cid:11)erent
de(cid:12)nition for F. Our main criticism of (3.15) is two-fold. Firstly, the arti(cid:12)cial
di(cid:11)usion model (3.15) will add more dissipation in elements where sharp features
of the solution are not orthogonal to the (cid:13)ow direction a, i.e. where a(cid:1)ruh is
large. However, less dissipation will be added where a(cid:1)ruh (cid:24) 0, i.e. in layers
that are orthogonal to the (cid:13)ow direction, cf. [12]. Secondly, the dependence of
n+1
(3.15) on the numerical solution uh leads to a nonlinear system of equations
even when the underlying partial di(cid:11)erential equation is linear.
9
n+1
Instead, we propose to de(cid:12)ne FjK for K 2 T , as follows
(cid:14) (cid:14)
(cid:15)^ 3 2 n (cid:15)^ 3 2 n
F(uh;h)jK = max(C1h Dhu(cid:22)h (cid:0)(cid:15);C2h Dhuh (cid:0)(cid:15);0); (3.16)
n n+1 n n (cid:15)^ (cid:15)^
where u(cid:22)h = uh(X(x;t ;t );t ) and C1, C2 are positive constants. The advan-
tages of this arti(cid:12)cial di(cid:11)usion model are as follows. Firstly, F is independent
of a (cf. [12]), so that sharp features of the solution will be treated equally,
irrespective of their orientation with respect to the (cid:13)ow direction. Secondly, F
n n+1
depends on uh, not uh , hence we only need to solve a linear system of equa-
tions when the problem under consideration is linear. Finally, the calculation of
n
F is relatively cheap since u(cid:22)h has already been calculated in the time-stepping
procedure of the Lagrange-Galerkin method.
We may now formulate the Lagrange-Galerkin method with the arti(cid:12)cial dif-
fusion model (3.16) as follows: (cid:12)nd uh such that, for n = 0;1;:::;M, uhjSn+1 2
hn+1
V and satis(cid:12)es
h (cid:22) hn+1
(Dtuh;v)n+1+(((cid:15)+^(cid:15))ruh;rv)n+1 = (f;v)n+1 8v 2 V ; (3.17a)
0 h0
(uh(cid:0);v) = (u0;v) 8v 2 V ; (3.17b)
n+1
where(cid:15)^= F(uh;h)andFjK forK 2 T isde(cid:12)nedby(3.16). Inthefollowingwe
shall refer to (3.17) as the `discontinuity capturing' Lagrange-Galerkin method.
4 A posteriori error analysis
In this section we shall extend the a posteriori error estimate presented in [12] to
the discontinuity capturing Lagrange-Galerkin method (3.17). However, before
2 2
we proceed, let us introduce the following notation: for v;w 2 L (0;T;L ((cid:10)))
we de(cid:12)ne
XM Z tn+1
(v;w)Q := (v;w)dt;
tn
n=0
1=2
kvkQ := ((v;v)Q) ;
where Q := (cid:10)(cid:2)I.
4.1 Error representation
The (backward) dual problem takes the form: (cid:12)nd (cid:30) such that
(cid:0)(cid:30)t (cid:0)r(cid:1)(a(cid:30))(cid:0)(cid:15)(cid:1)(cid:30) = e (cid:17) u(cid:0)uh; x 2 (cid:10); t 2 I; (4.1a)
(cid:30)((cid:1);t) = 0; x 2 @(cid:10); t 2 I(cid:22); (4.1b)
(cid:30)(x;T) = 0; x 2 (cid:10): (4.1c)
10
In the following theorem, we give an existence and uniqueness result for problem
(4.1), which guarantees su(cid:14)cient regularity of the solution (cid:30) to enable us to
perform the proceeding a posteriori error analysis (see [12]).
2 2
Theorem 4.1 Suppose that e 2 L (0;T;L ((cid:10))), where (cid:10) is assumed to be con-
1 2 2 1
vex; then (4.1) has a unique (weak) solution (cid:30) 2 H (0;T;L ((cid:10)))\L (0;T;H0((cid:10))
2
\H ((cid:10))).
Multiplying (4.1a) by e and integrating by parts in both space and time, we
get the following error representation formula:
XM Z tn+1
2
kekQ = (e;(cid:0)(cid:30)t (cid:0)r(cid:1)(a(cid:30))(cid:0)(cid:15)(cid:1)(cid:30))dt
tn
n=0
XM Z tn+1 XM Z tn+1
= (et +a(cid:1)re;(cid:30))dt+ ((cid:15)re;r(cid:30))dt
tn tn
n=0 n=0
XM
n n 0
(cid:0) ([uh];(cid:30)(t ))+(u0 (cid:0)uh(cid:0);(cid:30)(0))
n=0
XM Z tn+1 XM Z tn+1
= (f (cid:0)a(cid:1)ruh;(cid:30))dt(cid:0) ((cid:15)ruh;r(cid:30))dt
tn tn
n=0 n=0
XM
n n 0
(cid:0) ([uh];(cid:30)(t ))+(u0 (cid:0)uh(cid:0);(cid:30)(0));
n=0
h
where we have used (3.8a). If we now let (cid:8) 2 V , then using (3.17) we have
XM Z tn+1
2 n
kekQ = ([uh]=kn+1 +a(cid:1)ruh (cid:0)f;(cid:8)(cid:0)(cid:30))dt
tn
n=0
XM Z tn+1
+ ((cid:15)ruh;r((cid:8)(cid:0)(cid:30)))dt
tn
n=0
XM Z tn+1
h n
+ (Dtuh(cid:0)([uh]=kn+1 +a(cid:1)ruh);(cid:8))dt
tn
n=0
XM Z tn+1 XM Z tn+1
n n (cid:22)
+ ([uh]=kn+1;(cid:30)(cid:0)(cid:30)(t ))dt+ (f (cid:0)f;(cid:8))dt
tn tn
n=0 n=0
XM Z tn+1
0
+(u0 (cid:0)uh(cid:0);(cid:30)(0))+ ((cid:15)^ruh;r(cid:8))dt
tn
n=0
h
= I+II+III+IV+V+VI+VII 8(cid:8) 2 V : (4.2)