Stabilized Galerkin methods and layer-adapted grids for elliptic problems

Computer methods in applied mechanics and engineering ELSEVIER

Comput.

Methods

Appl. Mech. Engrg.

166 (1998)

165-182

Stabilized Galerkin methods and layer-adapted problems R. Hangleiter, University

of GBftingen, Mathematics Received

12 November

grids for elliptic

G. Lube*

Department, 1997 revised

D-37083 1 January

Giirtingen, Germany 1998

Abstract Global error estimates are considered for stabilized Galerkin finite element methods to second-order elliptic boundary value problems with emphasis on the singularly perturbed case. The critical point is the appropriate choice of numerical damping parameters. Furthermore, we address the resolution of exponentially decaying boundary layers with anisotropic finite elements on layer-adapted grids. The numerical analysis is based on both anisotropic interpolation estimates and sharp estimates of derivatives 0 1998 Elsevier Science S.A. All rights reserved.

1. Introduction In this paper, stabilized Galerkin methods are applied to elliptic boundary value problems of second order (modelling advection-diffusion-reaction problems), with emphasis on singularly perturbed problems. These problems are basic models in fluid mechanics. In particular, they may appear within the iterative solution of coupled incompressible Navier-Stokes problems, e.g. as the energy or mass transport equation or (as a system) in the transport step of projection type methods [5,16]. Standard Galerkin finite element methods may suffer from numerical instabilities generated by dominant advection or reaction. Hughes et al. introduced the concept of stabilized Galerkin methods (cf. Section 2 for details) which are already well accepted (cf. e.g. the monographs [ 18,191). These methods combine, as compared to the basic Galerkin method, improved stability and accuracy due to a residual formulation. The numerical analysis of such methods is more or less restricted to quasi-uniform or isotropic meshes so far. On the other hand, in singularly perturbed cases, it is often desirable to resolve boundary or interior layers. The goal of the present paper is two-fold: (i) to give a critical review of standard stabilized Galerkin methods on quasi-uniform meshes, in particular of the choice of the inherent damping parameters; (ii) to resolve boundary layer regions using anisotropic mesh rejnement (at least in simplified situations) together with an appropriate design of the damping parameters. These results are essential tools to derive a-priori error estimates (in a stabilized energy norm) which are uniformly valid w.r.t. the singular perturbation parameter. An outline of the paper is as follows: In Section 2 we review a class of stabilized Galerkin methods, derive the basic quasi-optimal estimates (cf. Theorem 2.3) and give a priori error estimates on quasi-uniform meshes (cf. Theorem 2.4). Section 3 contains a framework of a priori generation of layer-adapted grids and (anisotropic) interpolation estimates for exponentially decaying boundary layers. This approach is discussed in more detail in Section 4 for typical two-dimensional problems. Some numerical results are reported in Section 4 as well. For a subdomain G c fi we denote by Wk.“(G) the standard Sobolev space of functions with derivatives of

* Corresponding

author.

00457825/98/$19.00 0 1998 Elsevier Science S.A. All rights reserved PII: SOO45-7825(98)00089-9

order Sk belonging to Z,“(G). The norm resp. seminorm on W”“‘(G) are denoted by /I. /li,p,G resp. 1. ll,,,,(;. (., .)(; is the inner product in L’(G). In case of G = R we usually omit the index R. C denotes a generic constant not depending on singular perturbation and discretization parameters.

2. Stabilized

Galerkin

2.1. Continuous

methods

for elliptic

problems

problem

Let R CR”, d S 3, be a bounded domain with a Lipschitz Consider the following elliptic boundunly value problem

continuous

boundary

?)R and the outer normal ri.

L,u:=-eAu+h’.Vu+cu=f

ini

(I)

u=o

on df1.

(2)

Of particular interest is the singularly perturbed case 0 < E << 1 where the solution of (l)-(2) is mainly characterized by the solution u(, of the limit problem for E = 0. Different kinds of interior resp. boundary layers of the solution of (l)-(2) can appear in subregions where LQ,is not smooth Lesp. where ug does not satisfy (2). In particular, boundary layers may appear at outflow parts r+ of aR where b . n’ > 0 or at characteristic parts & with c.Fi= 0. Throughout this paper, we state the following assumptions for problem (l)-(2): (H.la)

O
1,

GE w’qn)“,

(H.lb)

(i) c:=c-iV.g20,

c EL==(R), (ii)

3~>0:

A weak solution u E V := W,#j.2(fl) of (l)-(2)

f’E L’(O) Ic-V.btl~~(:;

a.e.

in R

satisfies

Find u E V such that for all u E V: B,(u, u) = (.f; u) , B,(u, u) := (al,

(using integration

2.2. Stabilized

Vu) +:

((i+%L, u) - (i+4b,u)) +

by parts of the advective

Gulerkin

(3)

term). The standard

(4) energy norm reads

methods

Y!! = {K} with We consider Lagrangian elements on simplices KC Rd of an admissible triangulation 0 = U, K. II;(K) is the set of polynomials of maximal degree 1 2 1 on K. The Lagrangian interpolant of a continuous function u is uniquely determined by (I~‘u)(x,) = u(x,) for all nodal points of K. Let & C V be the finite-dimensional subspace of conforming finite elements such that V,,l, C 14(K). The standard Gulerkin method reads Find all E V,#such that for all uh E V,: B,(u,,

u,,) = (,j u,,)

(6)

For later extension on possibly different sizes of the element in different directions, we introduce the following notation: For KC R2 let EK be the longest edge of K. Then, we denote by h,,, := meas its length analogy. and by h,,, := 2 meas(K)/h,,, the diameter of K perpendicular to E,. In the 3D-case we proceedby Again, let EK be the longest edge of K, and let FK be the larger of the two faces of K with EE:CF,. Then, we the length of E,, by h2,K : = 2 meas /h , ,K the diameter of FK perpendicular to denote by h,,, : = meas /(h, ,Kh2,K) the diameter of K perpendicular to FK. Then, hold h , ,K s . 2 h,,,. EK’ and by h,., : = 6 meas Later on, an element K will be called isotropic if /z,,~ = h,,, resp. anisotropic if h ,,K << h,,,. Throughout this paper, for interpolation estimates on an element K we require the

R. Hangkiter,

G. Lube

Maximal angle condition There is a constant y* < r (independent angle yK of any element K is bounded by y* : yK s y*. The 3D-condition can be formulated by analogy interpolation result on arbitrary elements:

of h and K E &)

0 cm

W))k

= (A u)

methods

of residual

interior

the following

s 1 + 1, v K E & .

In the singularly perturbed case 0 < E << 1, it is often useful to stabilize unphysical Galerkin solution. Under assumptions (H.l), a solution u E V satisfies L,u =,f Q(u) EL’(a) the residual formulation (L,u -A

s.t. the maximal

(cf. [3]). Then, we have with 7, := u - Ir’u

V u f W’+‘.‘(K) : \jq>d~jm,z,K s Ch;l’,l -“‘lu/,+ ,,2,k

B,(u, u) + c

167

I Comput. Methods Appl. Mech. En,qrg. 166 (1998) 165-182

(7) oscillations of the discrete in L*(a) and thus with

(8)

v uE V

K

Stabilized

Galerkin

Find U,, E V,, such that B,,(lJ,, B,&

w) := B,(n, w) + c

type start from (8):

v) = FsG(u)

(L,u, $(w)),

V u E V,

(9) (10)

3

K

&,cw)

:=

.ftw)

+

c

(A

@tw))K

(11)

.

K

We consider

the following

J,@J,:=

G,(&+V)u, + yK(-eAu/,

with a unique (H.2)

class of methods

with

fcu,),

sign of all yK and the following

3 0 E (0, 0):

(i)

lS,Cf

Oa 1x1 aa,

(12)

minimal design properties 0 (ii) 6,~’ s 2 C a.e. in fi

with C, from ( 16);

d z h:,,

where 0 = 3 if inf, yK * 0 and 0 = i if supK yK d 0. (H.2) is valid for the streamline upwind method (SUPG) with SK 2 0, yK = 0, the Galerkin/Least-squares method (GLS) with SK = yK 2 0 and the Douglas-Wang method with SK = - yK 2 0. We recover the Galerkin approach with SK = yK = 0. The critical point is the choice of the sets {SK} and {yK} such that (9)-( 12) yields a stable and accurate method. 2.3. Quasi-optimal We introduce

estimates

on arbitrary

the following

]~~u~~~~~ := l]lull1~ + 7

stabilized

sKll”vuII&:K

grids energy norm l/l . lllSG defined by + c

max(O;

sgn(yK))ll

-

E Au

+

(13)

cuIIi.2.K

K

and T:(u)

:= c

8,/l&+. VUII(:,2;K ;

K

T&(u)

:= 2

/YK/ /i--E AU + cuii;,,;,

.

(14)

K

Additional stability of the skew-symmetric part of the operator (or weighted control of the streamline is the main effect of the stabilization. We need some auxiliary results. LEMMA

2.1. For u, E V, we obtain under the assumptions

B,,(u,,,

derivative)

(H. l), (H.2):

u,,) 2 Co(Wliu&

with C,(0) = 1 - G/2 if inf, yK SO and C,(B)=1 -a(0), o(0)=i(t3+-) Consequently, the stabilized schemes (9)-(12) are uniquely solvable. PROOF. (i) The inverse estimate (with C, = 0 if 1 = 1)

(15) if SUP,Y,~O.

168

R. Hangleiter.

llA4tw

condition

zz

G. Lube

I Comput.

Methods

Appl.

Mech.

Engrg.

166 (1998)

16.5-182

v?iN I.2.K 3

(16)

IyKyK( s 6, and (H.2) imply

T;,,(u) =S;

&II - E Au + c&,.,

1< 2 T (~%K(IAbll;,~.~ + %llcr&.~)

d ~lii~lli~

(17)

Now, set w = u in (9)-( 12): B,,(u,

u) 2 111~111; + c

K

S,~~~~VU~~;,,., + c

+ c Y&-E K

K Au + cul/,:,,,

&-&Au

+ cu, ~VU),

+ (/T-VU, -E Au + CU)~)

(ii) In case of inf, yK 2 0 we find using ( I3), ( 14) ( 17) and Young’s

B,,(u, u) 3 -(l

-v

lllulll:.~;

-

&

inequality

-&JT;~,(U).

Equilibrating the constants with CY* = (Y,/0 = 1 /v/l (iii) In case of supK yK S 0 we similarly obtain

+ 0 yields the assertion

of constants

Now, we introduce

LEMMA that

with (Y= CX,= CY~and a” ~ &Y ~ H = 0 yields the assertion.

0

the notation

1

B, := II&z.,

(with a; > 0)

111~111~; -&c&d

-+)Tli(al+&

Equilibration

-&)lllull,;+(l

ad

(18)

c, :=

I14,r,K9

2.2. For II E V with Au E L’(K),

G := IlCll”.X.K >

Z, := min(6;‘;

v K E 9jT and u E V,, we

obtain

B&

under the

~‘) assumptions

(19) (H. I), (H.2)

(20) Qs,W

:=

~~~~~~l~ + %d + T’;,,W+ T Z&&I:

with K = max(3, K’) PROOF (i) We estimate

(21)

if inf, yK 2 0 vesp. K = 3 max(2, K’) (f supK yK d 0.

the different

terms of

I~,,(~, rJ)l s III + III1 + III11

(22)

separately using standard inequalities with appropriate via integration by parts, (H. 1b), (14) and (19) that III:= IB,(u, u)/ = Ie(Vu, Vu) + ((c -v. =

ElUl &I,.~ + ‘d

constants

L8 > 0: For the Galerkin

part, we obtain

6)u, u) ~ (u, 6vu)l

+fW ~~~~o.,~~~~~~~o,, 11~11”.,.,~111~111~ + ew

c max(l;KNI~II,III~III, +~~L44

+ T?(u)

R. Hangleiter,

G. Lube

I Comput. Methods Appl. Mech. Engrg.

169

166 (1998) 165- 182

(23) For the stabilized

terms, we find (with the obvious

(II\ := 1c

A u + cu, lhv),

S,{(-E

K

+ (6Vu,

definition

of TsS, cf. (14)):

6*Vv),}

(24)

1c yjJ-

lIIIl:=

E Au+cu,-•Au+cu),+(~.V~,-•Au+cu),}

K

(25) Summarizing

(22)-(25),

we arrive at

(26) > 0. We obtain the assertion (ii) Let inf, yK _supK yK s 0. Using (17), we have L, + L, + B(L, + L,)

IB,,(~~ 4 =s

Ill4ll: +

2

where the dots summarize the remaining i=l,... ,6 and using d<+. Cl For the stabilized THEOREM Galerkin

2.3.

methods

(9)-(12),

The assumptions

methods

with L, = L, = 2L, = 2L, = L, = L, = iC,,(B).

4 + 4 + L4T;$) 2

+ *- *

(27)

in (26). The desired result follows with L, = +Co(0),

u-terms

then we have the following

(H. I)-(H.2)

Let now be

variant of Cea’s Lemma.

imply the quasi-optimal

a priori error estimate

of the stabilized

(9)-(I2)

GCQ,,(zd - Z:(‘u) lllu- UhIIlSG

(28)

where C denotes a constant which is independent

ifs,=&-0, PROOF.

of E, h, 8,. Note that in particular

Qsc(v) := lllulllc and C = 1

G=G. Let

e,,:=C/h-U=(Uh-Z:j’U)+(I~‘L1-U)--Xh+rlk

where Zl:’ : V + V, denotes consistency of the methods,

the Lagrangian we obtain

interpolation

operator.

colllx,IIl%G d &LYh~x/J= 4&I - 77hl XJ = -&&7h~x/J

The triangle

inequality

and (21) conclude

the proof.

0

By Lemma

2.1, Lemma

2.2 and

the

We are now prepared to give an error estimate on quasi-u@rm grids where h, = II ,,K = h,,, and to determine the parameters SK and yK. Note that we introduce (w.r.t. the anisotropic case in Section 3) a modified definition of the mesh P&let number Pe, : = h,,,B,e

’

(29)

THEOREM 2.4. Let the solution of (l)-(2) be smooth mccording to u E V n W ‘+“‘(0) und the triangulation
error estimate of the stabilized

s CT

From Theorem

(30)

h:(c + B,h,

+ &h;)lulf_,.z

2.3 and the approximation

Gulerkin

methods reads

/,. .

(31)

result (7) we conclude

using (H.2)(ii)

+ 7’%,) + T;,I(T,,)+ T Z,li&,) lll~,,lll.L c C(llIv,lllZ~ s Balancing

C c hf;Je K

+ C,hf,

+ 6,BZ + Z,hf,,)lulf,

,,2,K

(32)

the terms

6KBi -Z,hf,,

= ht., min(6.l;

Bi.6

we arrive with the mesh Peclet number

6, - h,,,Bi’

as in (29) at

8, - h:,Kem’

if Pe, 3 1 ;

Note that this definition

‘) ,

gives no contradiction

if Pe, s 1 to (H.2). This concludes

REMARK 2.5. 6) In the weakly anisotropic cuse h, ,K/hr,,K = : &. s r mesh-refinement, we have the modified estimate

the proof.

for all K E F,,, s.t. r *

0

1 which allows for local

(ii) In contrast to other papers (e.g. [7,15]), we did not include the reaction term in the design of the parameters 6, and yK. The main point is that, in the symmetric ca.se the quasi-optimal estimate in the energy norm Ilj . /IIc with the optimal constant C = 1 of the Galerkin method cannot be improved. For linear elements, i.e. 1 = 1, it is even wrong to use the GLS method, as numerical dissipation is then subtracted instead of being added to the Galerkin method. On the other hand, we did not find an improvement with the Douglas-Wang variant in numerical experiments 121 for such problems. 0 2.5. The singularly

perturbed

cuse

In the singularly perturbed case 0 < E << I, the estimates of Theorem 2.4 are, in general, not satisfactory. The interpolation estimate in (28) requires information on local Sobolev norms of the solution which are in general not uniformly bounded w.r.t. E + +0: /U ~ I;,“~ll,,,.~:K c Ch;‘;m”‘Iul,+ ,,2,K = Ch;+,lm”‘?‘), One advantage

of stabilized

Galerkin

methods

s(1) >O.

is the high-order

accuracy

away from boundary

and interior

R. Hangkiter,

G. Lube I Comput. Methods Appl. Mech. Engrg.

166 (1998) 165-182

171

layers where the solution is smooth. This has been proven for simplified problems in [20] and e.g. in [ 131 for reaction-diffusion resp. advection-diffusion problems. Different methods have been proposed to remedy (restricted) oscillations of the discrete solutions of stabilized Galerkin methods in boundary or interior layers. For methods of shock-capturing type cf. [4] and the references given there. On the other hand, we are not aware of any rigorous convergence analysis of such methods. In numerical experiments we sometimes did not observe convergence with h + +0 on isotropic meshes. Other methods using stabilizing terms of higher order [8,22], require regularity of the continuous solution in an unrealistic way. So we decided to use stabilized Galerkin methods as defined in Section 2 together with an appropriate resolution of boundary and interior layers.

3. A framework

of anisotropic

a priori layer refinement

In the remainder of this paper, we combine in the singularly perturbed case stabilized Galerkin methods with the resolution of layers. This is often an important point in practical computations. An isotropic resolution is, in general, too expensive due to the lower-dimensional character of layer phenomena, an anisotropic approach is much better suited. This approach is rather new in the literature and far away from being complete [19,6]. In Section 3 we consider basic aspects of such an approach in the special case of exponential boundary layers occurring at a (d - 1)-dimensional edge resp. face 2 of a polyhedral domain 0 C Rd. The ingredients are the construction of luyer adapted grids and the (anisotropic) interpolation of boundary layer functions together with the design of the discretization parameters 8, and yK. In Section 4 we apply this approach in more detail to 2D-model problems. 3.1. Generation

of hybrid layer-adapted

grids

Let h be a user-given global mesh size which is a-priorily set to resolve the solution away from layers. The following steps of the a priori generation of hybrid grids are proposed: l Detection of boundary layers at a manifold 2 C a0 and of its width a, l Generation of a structured mesh in the layer strip Q(s) of thickness (~?(a,) with possibly large aspect ratio h, K >>hd.K Generation of an (possibly unstructured) isotropic mesh of size h ,,K - h,,, - O(h) away from layer regions. To be more precise, we discuss the following simplified situation to be considered below. Here, we assume that a, d h. l

EXAMPLE 3.1. Let 2 C an be an edge resp. a face in the 2D- resp. 3D-case. Consider an edge- (resp. face-) fitted Cartesian coordinate system (4, p) with p : = dist(x, 2) and tangential variables 4 : = (4,) . . , 4<,,,-_, ). A boundary fitted layer-adapted mesh is built in the following way in the strip Q(z)), i.e. 0 c p d uI (cf. Fig. 1):

Fig. 1. Adapted

grid in a boundary

layer strip

172

R. Hungkiter,

G. Lube

I Cornput.

Methods

Appl.

Mech.

Engr,y.

lb6 (1998)

165

182

Setting p=p,, i=O ,..., N+ 1 with p,,=O and pNI, --u:, we construct a quadrilateral resp. hexahedral mesh in the strips %“‘(s): p,_ , s p up, with corners in the lines resp. planes p = p, s.t. a maximal angle condition (similarly defined as in Section 2.2) is valid. Furthermore, suppose that the resulting grids at p = p, consist of isotropic elements of size O(h). Assuming p, - p, ~, << h (to be verified below), we subdivide each quadrilateral resp. hexahedral element into simplicial elements K in such a way that the maximal angle condition and the coordinate system condition (cf. Section 3.2 below) are fulfilled. Then, set h,,, := p, - p,_ , = ui’h with a’,’ to be determined later. REMARK 3.2. 6) The distribution of (anisotropic) strips parallel to -?’ IS not necessary. But it turns out from numerical experience, that the longest edge of the finite elements has to be carefully oriented w.r.t. 2, at least in the neighborhood of characteristic layers. (ii) The proposed u priori construction of layer-adapted grids simplifies the presentation but it may also be incorporated within an uduptive method where the global mesh size h can be modified based on appropriate a posteriori estimates or error indicators. (iii) Furthermore, the restrictions to straight manifolds generating the layers and to boundary layers are not necessary. But an adaptive method with anisotropic resolution of boundary and interior layers has to be implemented with (cf. e.g. [21]). 3.2. Anisotropic

interpolation

estimates

The following condition is useful to describe the orientution of the layer-adapted grid and to derive anisotropic interpolation estimates w.r.t. to an appropriate fixed Cartesian coordinate system (+,, , 4<, , , p). Coordinute appropriate

system condition (20): The angle r,& between the longest side of element fixed Cartesian coordinate system) is bounded by lsin &.I s Ch,,,/h,,,.

K and the &axis

(of an

The 3D-coordinate system condition can be formulated by analogy (cf. 131). It is satisfied in the situation of Example 3.1 w.r.t. the (4, p)-system. The following interpolation result takes advantage of different directions in K using the multi-index notation Ly=(a ,/...)

LY<,)EN;:,

IQI = c

a, ?

and derivatives D” w.r.t. the fixed coordinate interpolant in c, C V := W’*2(fl).

h;:=h;‘l.. system.

.h;;”

Furthermore,

we denote

by f:,“(. ) the Lagrangian

THEOREM 3.3 ([2]). Assume that a finite element K E y,-,, satisfies the maximal angle resp. coordinate system conditions. Furthermore, let be v E W”‘“(K), k E { 1, . ,I + l}, p E [ 1, ~1 s.t. p > 2/k. Then, we obtuin for fixed m E (0, , k - l} Iu - $~i,,.,:,

c C

c h;lD”vl,.,:, !nl:k-,,,

.

In view of Theorem 2.3, we derive an upper bound of the term Qsc(v - iz)v) as in (21) for a smooth function v on a mesh satisfying the coordinate system condition w.r.t. a fixed coordinate system. COROLLARY 3.4. Under the ussumptions of Theorem elements K E Y,,, and with v, := (I - iy’)v that

3.3. we obtain for u with vIK E W’+“x(K)

for all

(34) (35)

R. Hangkiter, PROOF.

Theorem

lrl,l~,,;,

s

3.3 implies

C , ,F,_ a

+

m

G. Lube

I Comput. Methods Appl.

Mech. Engrg. 166 (1998) 165-182

173

for v, := (Z- rl(‘)u with k = 1 + 1 and p = m that

hi? meas(fWX,m,K

9

hence with EK;p,v as defined in (35) follows

(36) REMARK

3.5. For anisotropic

3.3. Interpolation

interpolation

of exponential

estimates

on quadrilateral

elements

see Ill.

layer terms

Now, we apply the anisotropic estimates to the interpolation of exponential boundary layer terms on a layer-adapted grid at a (d - 1)-dimensional straight manifold as discussed in Example 3.1. To be more precise, consider the boundary-fitted local coordinate system (4, p) and let the following exponential decay condition be valid: (ED)

For given numbers

r > 0, (T E (0, 11, for Ial s 1 + 1 and each element

K E a(s)

holds

dist(K, 2)) .

]l~“u\l~,~_~ d ~~~~~~ exp(-rem”

(37)

This condition implies with u0 > 0 that ]ID”u~]~,~,~c CE V’1if dist(K, 2) 2 ((~(1 + 1) + uO}TP’e”llog E]. The r.h.s. value may be defined as the layer thickness az. Using an appropriate cut-off function, we suppose that the layer term u vanishes in n\%(s). Then, we have the following global interpolation estimate: LEMMA 3.6. Let a layer-adapted grid be constructed in %(x) fl R as in Example 3.1 with (possibly anisotropic) elements K with hj,k = O(h), j d d - 1 and h,,, = a,h, aK > em. Suppose th;; a boundary layer term u E W’+“= (0) satisfies (ED). Then, we obtain for the interpolation error 7, := (I - I, )u that (with EK:P,y defined in (35)) Q&J

6 C c

Fk := c

(38)

xK ak-2YdEK;p,y ,

IPl=l

PROOF.

h2’+dFKGK,

Corollary

Q;&,> c

G, := exp(-21’-”

,

dist(K, ~))(a,e~“)*“+”

(39)

IYl’l

C

3.4, condition

(ED) and the special mesh construction

imply

cK mm(K) c c c EK,p,yh~a+P)IIDa+P+yull~.~,K la/=l-I

lpl=l

Iyl=I

d-l

< Cc

2 c K lal=1-1 IPI=

c

EK;P,y G

XE -2u(ad+Pd+yd)exp(-2J’mc c C c K

h2’+d exp

h2~a~~P~~~‘(aKh)2~a~~Pd~i’

jyl=l

-2rdist(K, EC

The idea now is to use G, for an appropriate

dist(K, 2)) 2)

aK *(r+l) I(-)

Err

c

lPI=lvl=l

aK

choice of a, and to minimize

‘-2ydEK;p,u

.

cl

the expression

FK with respect to 8,

(resp. YK). DEFINITION 3.7. A layer-adapted mesh in the strip Q(x) supK G, S 1 (with G, as in (39)) is called layer-resolvent.

fl L! as constructed

in Example

3.1 and satisfying

174

R. Harqleiter,

G. Lube I Comput. Mrthod.s Appl. Mech. Eqr,q

EXAMPLE 3.8. Condition G, s I can be satisfied recursively strip “u”‘(x), i.e. p, , C p := dist(K, C) c p,, holds

a

h

11) ._

166 (1098) 16.5-182

(with increasing

i) if for all elements

K in the

~ P, - P, I h

K’-G-

(40)

(1) In particular, in the strip %’ ” nearest to 2, one has to take aK = aK = E(‘. One should of course set the values ai’ in (40) as large as possible, i.e. c,,._P,-P, aK .-

I

(41)

h

Setting formally r = 0, hence aK = ug’ = E”, we get a mesh of so-called Shish& type. The number of anisotropic layers %‘I’ is then given by N- (1 + 1)h ‘Ilog ~1, h ence the number of (anisotropic) elements in the layer a(c) is of order O(h-dllog ~1). Such meshes take no advantage of the exponential decay of layer functions. A mesh satisfying (4 1) with I’ > 0 leads to a logarithmically graded mesh of so-called Gartland type [ 191. The number of layer strips 02c(‘I is considerably smaller than for a Shishkin mesh. THEOREM the choice

qf Lemma 3.6 he valid. Furthermore,

3.9. Let the assumptions

IyKyKI s 8, = min(h,,,B,‘;

let the mesh be layer-resolvent.

h:,Kc-‘)

With

(42)

we ohfain Q&(ll,)

c C c

h”+” 2

K

Ch”llog

C

(e + BKh,,, + C,h;,,)

(431

d,K

Elmax (E’-” + B,h + C,hh,,,)

.

(44)

K

PROOF.

Starting

c l!a=lvl=~ This implies c

from Lemma 3.6, we estimate

1-2yc,_I.K aK

h h d.K ’

c

hza;

different

terms of FK separately:

‘vti - h,,,h,,,

IPI-lvl-1

,

c lPl=lvl=~

using (H.2)(i) EK,p,ya~p:-2y2S C 2

E, ;

E, := E + e,h;,,

+ 6,B:, + Z,h;,,

.

(45)

lPl=Irl=l

proceed now with the minimization of E, w.r.t. 8, as in the proof of Theorem 2.4. Finally, ratio holds h,,,/h,,, c E-“. The number of elements is bounded by h-“[log ~1. 0 We

for the aspect

REMARK 3.10. (i) The estimate of Theorem 3.9 remains valid if we add to a layer term u satisfying (ED) a function w E w/+I.r =ZC f C(E) for all IcyId 1 + 1, e.g. the global part U,,,, of an asymptotic W with II~“410.~.n expansion of solution u (cf. Section 4). (ii) The interpolation estimate of Theorem 3.9 is of order I and almost uniformly valid w.r.t. E. Note that there is a gap of Q((E~) in the first r.h.s. term as compared to estimates of smooth solutions away from the layer region (cf. Theorem 2.4). 1 as an upper bound of G, with obvious changes in (iii) As an alternative, one can choose (I + 1)llog l Example 3.8 and Theorem 3.9. In particular, the number of elements in a Shishkin type mesh reduces to 0(h -“) but the logarithmic term in (44) still remains.

R. Hangkiter,

3.4. Condensed

G. Lube

grids at geometrical

I Comput. Methods Appl. Mech. Engrg.

175

166 (1998) 165-182

singularities

A condensed mesh (not necessarily of isotropic type) appears if different layer adapted regions “zc(&) match at convex geometrical singularities Y = i-l:, , i5,, 1 d S cd (corners resp. edges and comers of a polyhedral domain in R2 resp. R’). Special layer terms can compensate perturbations arising on Q(Y) n a0 by such interacting layers terms. We neglect here a (possible) singular behaviour of the solution caused by incompatibility of the data at .!X For simplicity, let us assume that a Cartesian coordinate system (x,, . . . ,x4) is adapted to Y such that the edges resp. faces _&, . . . , xs are located at x,~= 0, s = 1, . . . , S. Furthermore, assume that a layer term z satisfies at Y the following exponential decay condition: (EDC)

For given numbers fI;=, %(&zl, holds

$ > 0, cryE (0, 11, 1 GsSSSd,

for lal
1 and each element

KEQ(Y):=

(EDC) implies certain data compatibility conditions at Y to guarantee z E W’+“=( %( 9)). Suppose again that, using suitable cut-off functions, the layer term z vanishes in n\‘%(Y). Note that h,,, has to be chosen according to the layer %(A,) with maximal u,. A straightforward calculation (using similar arguments as in the proof of Theorem 3.9) yields that the result of Theorem 3.9 remains valid (with obvious modifications). We skip details but give some remarks for simple examples in Section 4. REMARK 3. II. A detailed calculation shows that the mentioned modification of Theorem 3.9 for corner or edge layer terms with (46) remains valid if we add a smooth function (as in Remark 3.10 (i)) and exponentially decaying layer terms (according to Lemma 3.6 and Theorem 3.9) on a condensed mesh at the boundary of a hyperplane 2 as discussed in Example 3.1.

4. Application

to two-dimensional

model problems

The goal is now to apply the preceding results to model problems on the unit square containing different types of exponential boundary layers. Therefore, we take advantage of asymptotic expansions of the solution which are available in some cases. Furthermore, we support the theoretical results by numerical examples. 4.1. Asymptotic

expansion

of the continuous

solution

In the singularly perturbed case, the behaviour of the solution of (l)-(2) can be arbitrarily complicated, depending essentially on the behaviour of the characteristics of the limit operator L,. So, for brevity, we restrict ourselves to problems L,u:=

--EAu+g.Vu+cu=f

in 0=(0,

1)2

on a.0

U=O

(47) (48)

with sufficiently smooth data and with a highly organized flow field b’such that the asymptotic structure of the solution is relatively simple. Additionally, we assume regularity of the solution s.t. u E ,‘+““(a) for 0 c E c 1, hence interior layers are avoided. Moreover, let on any edge 2 hold either 2 c f, with lbt. c]lrn 2 p > 0 or .I5 c 4, i.e. b’. GIZ = 0. The latter assumption together with smooth data considerably simplifies the different possibilities of boundary layer behaviour at f+ and r,. We discuss three typical situations below. For the numerical analysis of stabilized Galerkin methods on hybrid meshes we propose a modified application of the quasi-optimal estimate of Theorem 2.3 using an asymptotic expansion u = uy + rM of the solution u with a sufficiently small remainder rM:

lb - ~,lllscc c{Q,&~

- Il(‘ut)

+ Q&,

- Ijj’r,,,)l.

(49)

176

R. Hangleitrr,

G. Lube

I Comput. Methods Appl. Mrch. Engrg.

166 (1998) 16.5%182

In general, it is much easier to find estimates of the derivatives of the asymptotic expansion. On the other hand, we require only a low order estimate of the remainder. A standard expansion of problems (47)-(48) under the assumptions stated above has the following structure [19,10]:

=

,;)E’U,(X)+ 2 E’%,(X)+ 2 lz,(x) ,=c

+ r&x,

E)

, -0

lJ, denotes the global expansion or regular part of the solution. Boundary layer expansions (of outflow or characteristic type) are collected in V,. At corners of the domain, corner layer expansions Z, occur typically if there holds V, + U, # 0. The global part U,,,,:= Et0 E’u,(x) solves recursively the system

1

f L,,u, := b--vu, + cu, =.f; := ‘Au

j=O j= l,...,

,~,

in .Q\r

M

with homogeneous Dirichlet conditions on the inflow boundary 11 only. Interpolation estimates of r/, can be derived from Theorems 2.4 resp. 3.9 (or Remarks 3.10, 3.1 I) for subdomains away from resp. in boundary layer regions, provided that one additionally has (H.3)

U,,,,E

W’+‘,“(fI) .

In particular, no interior layer originates at points of ? chosen away from layer regions as indicated by Theorem 4.2. Interpolation

of boundary

The global expansion characteristic boundaries I;,. Following Example neighborhood ?&EC(x)of one-to-one mapping (x,,

layer corrections.

Or+ U 4,. As a consequence, 2.4 on isotropic meshes.

Numerical

the sets {x}, {SK} are

results

U, does in general not satisfy the boundary conditions (48) at outflow resp. r, resp. r,. According to our assumptions, an edge 2 C aR\T belongs either to r+ or 3.1, we introduce a layer-adapted mesh and a local coordinate system (4, p) in a 2 where p(x) := dist(x, 2) and 4 denotes a tangential variable on 2. The smooth .x2)-(4, p) transforms the elliptic operator L, to

L,u : = - ELZU+ L,,u

with A(& p) := X2_, (ap/ax;)’ coefficients at p = 0 yields

> 0; B,(& p) := g-V& B,(4, p) := g.Vp, B,(+, p) := c. Taylor expansion

of the

K >

B,(~>P)-C

5 : = P/E“,

we determine

A(& P) = ,$) Al(+> 0)~’ + &I”) Introducing

the transformation

I =o

B,,,W)P’+O(~~~~L

j=O,

1,2.

(52)

g E [0, l] such that

(53)

4.2.1. Outjlow layers The simplest case appears if _X is part of the outjIow boundary r+ s.t. B2,“(& 0) = 6.Vp = -b’. $1, s p < 0. Then, (53) results in cr = 1 and the representation (with K := M) . h.o.t..

,

a2 +. L,, .= -A,, 2 + 4, al

a al

R. Hangkiter,

G. Lube

I Comput. Methods Appl. Mech. Engrg.

The formal expansion of the boundary recursively from the system

layer

L,:u,+ = - C Lk+u,Tpk,

L,+v,+ =o;

correction

j=

l,...,

V, := Z,co u,‘(@, l)d

M

171

166 (1998) 165-182

can

then

be determined

in %(.Z)

k=l

with

u,+

+

U,

=

0

on Z: We obtain sufficiently

smooth solutions

(via smoothness

of the data and of U,,,) of the

form

In conclusion, mind &,,I(+,

the boundary

layer term V, is of exponential type satisfying ence, application of Theorem 3.9 results in O)l. H

O)/A,(&

1x1 d 6, = min(h,,B,‘;

(ED) with (T = 1 and 0 < r <

h&6-‘).

In particular, in the strips nearest to the edge, one obtains 8, c C&*. On a Shishkin type mesh, this is even valid in the whole layer strip. Therefore, the artificial diffusion of stabilized Galerkin methods in outflow layers is much smaller as dictated by the standard choice on isotropic meshes. One can even set 6, = yK = 0 on a Shishkin mesh (cf. Example 4.1). EXAMPLE

4. I. Consider

u = 2 exp(-x, ,=I

problem

(47) with b’= (1, l)T, c = 0 and the exact solution

lc) - ,Ij exp(-xi

/E)

with vanishing global expansion U,,, and outflow layers at x, = 0 and x2 = 0. The corner layer term II:=, exp(-x, /E) can be resolved by the condensed mesh at the origin. The remainder of the asymptotic expansion vanishes. Application of Theorem 3.9 to the outflow layers and Remark 3.11 to the corner layer yield the global (and almost uniform w.r.t. E) estimate llll.4- U,l&

s Ch2Jlog El

We resolve the outflow layers by a Shishkin type mesh. In Fig. 2 we present the convergence history in the energy norm and L”-norms (p = 2 and p = 0~) for E = lop3 and different choices of 6, = yK in the layer elements. First resp. second-order accuracy are observed. On the other hand, for moderate h, the error is

1o-2

-A.O

0

eo +xqah’ 45h

h

Fig. 2. Convergence

h in the energy norm and LO-norms (p = 2, CD)for EL 4.1 on a Shishkin mesh.

considerably an isotropic

smaller for 8, = 0 and 8, = mesh.

lh2

in comparison

to 8, = Ch which would be the standard choice on

4.2.2. Characteristic layers Assume now that an edge 2 C aR is characteristic, i.e. z’ C I;,. This implies B,,,,($+ 0) = 6.v~ = -6. n’j2 =: 0. Further, let k, be the smallest index in (52) such that B,,, $0. Under the assumption c(x) = B,,(4, p) 3 y > 0, we conclude

from (53) that LT= 1 /2. This implies

with I?!= 2M the representation (54)

(i) Reaction-diffusion problems The simplest case appears in (54) if 6~ 6, hence k, = k? = x and L:: := -A, fl’/aJ’ + B,,,,,Z. The formal expansion of the boundary layer correction V, := Zff(, u:,~(& L)E”’ can then be determined recursively from second order ordinary differential equations

L$J::= 0, and vi,+u,=O,

L;;v:’ = - i: L:‘v:‘, ) i(=I j=O ,...,

K; v(:,_, =O, j=

j=

1,...,2M l,...,

in a(Z),

M on 2. We obtain solutions

of the form

type satisfying (ED) with LT= I /2 i.e. the layer term V, is of exponential Bo,(,(& O)/A,,($, 0). Let us recall that we set yK = 8, = 0 in the symmetric case. EXAMPLE

4.2. Consider

u = E exp(-x, ITI

problem

and

0 < r < min+

(47) with h’= (0, O)T, c = 1 and the exact solution

/ &) - ,Q exp(-x,

&)

The solution has a vanishing regular expansion CJ,,,,,two exponential boundary layer terms and a corner layer term. The latter term compensates the influence of the exponential edge layer terms. Application of Theorem 3.9 and Remark 3.11 result in

which is uniformly valid w.r.t. E. In Fig. 3 we present the convergence history in the energy norm (unscaled and scaled) for a Shishkin type mesh and different values of E. We observe similar results as in Example 4.1 and additionally robustness w.r.t. E. Furthermore, the scaled results (energy norm divided by dm) indicate that this factor in the theoretical estimate is sharp. Now, we compare the layer resolution with a Shishkin and a Gartland type mesh on moderate tine meshes with 65 X 65 = 4225 resp. 43 X 43 = 1849 grid points. We observe that the error is nearly of the same order on both meshes (at least for small E). The Gartland type mesh is clearly preferable due to the smaller number of required (anisotropic) elements. Table 1 Comparison E

IO

2

10 h 10 ‘I’

of Shishkin

and Gartland

type grids

Mesh type

Energy norm

L’-norm

Lx-norm

Shishkin Gartland Shishkin Gartland Shishkin Gartland

6.430E-003 I .S60E-002 1.7I I E-003 2.254E-003 2.842E-004 2.68 I E-004

2.234E-004 8.072E-004 I .500E-004 4.191E-004 2.639E-005 3.101E-005

I .413E-003 5.406E-003 8.659E-003 2.132E-002 2.043E-002 I .5I 1E-002

R. Hangkiter,

G. Lube I Comput. Methods Appl. Mech. Engrg.

179

166 (1998) 165-182

h

Fig. 3. Convergence

in the energy

norm (unscaled

and scaled) for Ex. 4.2 on a Shishkin

mesh

(ii) Advection-diffusion problems The approach in (54) is more complicated if b’Z 6, b’. n’ = 0 at f, and min , = ,,2 kj is finite. Then, the simplest caseappearsfork,=Oandk,~1,henceL:::=-A,a2/a12+B,.,a/a~+B,,,Z.Theformalexpansionofthe boundary layer correction V, := CT:,, uy,?(4,J)e”’ can then be determined recursively from second-order parabolic differential equations

L,;;u:’=fy := - i: L$f_,

L$Ji =.f:: := 0,

)

j=

l,...,

2M

in Q(X)

k=I and Dirichlet conditions on J$ (i.e. at p = 0) resp. at (48). A smooth transformation to

4 =0

such that U, + V’ satisfies the boundary

conditions

is given by

6, := u, exp The leading

.

term has the form

Corresponding expressions can be found for higher-order terms 6: of the layer expansion. are sufficiently smooth if the (rather restrictive) compatibility conditions

akU (H.4)

s=O,

Ikl c 1

at

4 =0

are valid. Therefore, the boundary dm. Application of Theorem IykyKI d 6, = min(h,,B,‘; EXAMPLE

4.3. Consider

(resp. at L? n K)

layer term vM satisfies condition 3.9 yields

(ED) with u = l/2

hi k~-‘). problem

The terms ~7: resp. VT

(47) with g= ( 1, O)T, c = 0 and the exact solution

and 0 < r < min,

1/

180

K. Hungkiter,

G. Lube

I Cornput.

Methodc

Appl.

Mwh.

Engrg.

166 (1998)

165182

=y$yexp (- 4,(IXL) )

l4

with a parabolic Theorem 3.9

layer term at x2 = 0. No global and corner layer terms appear. The error estimate follows from

II[U- u/J;,

d CP( & + h)llog El .

In Fig. 4 we present the convergence history in the energy norm and L”-norms (p = 2, x) for E = IO-” and different choices of the parameters 6, in the layer regions. Again, we obtain first resp. second-order convergence for the energy resp. the L”-norms. As opposed to the case of outflow layers, the choice of 8, in the layer strip has apparently no significant influence. REMARK 4.4. The structure of the characteristic boundary layer terms complicated in the general case, for a detailed discussion cf. [lo]. 4.3. Remarks

on interpolation

arising

from (54) is much

more

of corner layer terms and of‘ the remainder

The (anisotropic) error estimates-as given in Section 4 so far-are in general not sufficient to get error estimates which are uniformly valid w.r.t. E. The$r.yt problem is a possible interaction of different layer terms at geometrical singularities (here at corners) of the boundary. In general, one has to construct special corner layer expansions Z,,, as solutions of elliptic equations in appropriately stretched variables (i,, I&). They have to compensate the expansion U, + V, in such a way that ~1; = U, + VM + Z, satisfies the given (homogeneous) Dirichlet conditions on aR. The second problem is that there appear at corners the same geometrical singularities of the solution as for standard elliptic equations [ll], but on E-dependent scales. Z,M is only sufficiently smooth (and satisfies condition (EDC)) if data compatibility conditions (of higher-order) are satisfied. Let us briefly discuss the arising problems. The simplest case appears for d#usion-reaction problems. The influence of corner layers is then only local, i.e. restricted to a 0( Allog cl)-region of a corner 112,141. Singularities can, in general, be resolved only using an appropriate geometric refinement. A resolution of Shishkin type is not sufficient. The situation is more complicated for advection-diffusion-(reaction) problems. In general, the effect of corner layers is not local in corner regions. Singularities of the solution caused by geometrical reasons at points P of r are distributed along subcharacteristics of the limit operator L,, passing through P. Thus, interior layers appear which could be resolved within an iterative process inbedded in an adaptive refinement approach 1211. We conclude with remarks on estimates qf the remainder term. Recall that estimates of Q,s,((l - I)(‘)r,)

0.0015

OOM 0.0035 0.003 00025

0.002 P e y o.m15-

Fig. 4. Convergence

in the energy norm and L”-norms

(p = 2, p) for Ex.

4.3

R. Hangkiter,

G. Lube

I Comput. Methods Appl. Mech. Engrg.

166 (1998) 165-182

181

require control of low order derivatives only. For reaction-diffusion problems, such estimates can be found in [ 141 for the genera1 case of polygonal domain. The advection-diffusion problem is again more complicated. For the simplest case of outflow layers in the unit square only, we refer to the discussion in [6]. a more general result is open so far.

5. Summary

and concluding

remarks

In this paper, we consider error estimates (in the energy norm) for stabilized Galerkin methods. After a review of such methods on isotropic meshes we focus the discussion on the resolution of boundary layers. The main point is a refined design of critical discretization parameters which depend on the diameter of the largest ball inscribed in an element K. The analysis is based on anisotropic interpolation estimates of different layer parts of the asymptotic expansion of the solution. This has been done for the special (but important) case of exponentially decaying boundary layers appearing at (d - I)-dimensional edges or faces of an polyhedral domain L2 C Rd. Furthermore, we discuss different types of layer-resolvent grids. Meshes of Shishkin type, recently introduced in the literature, are covered. Grids which are more adapted to the exponential decay of the layer (e.g. of Gartland or Bachvalov type) require much less grid points in the layer and give essentially the same numerical results. Numerical results for different exponential layer types confirm the theoretical convergence rate for the energy norm 1 (here I= 1). A theoretical foundation of the observed second order convergence rate (for I= 1) in the L”-norms will be considered elsewhere. A remarkable fact is that the discretization error is almost the same for different choices of the sets yK and 8, provided that h is sufficiently small. However, it is preferable to stabilize the Galerkin method (at least away from layers) in order to get a robust discrete problem. Open problems are the discussion of more complicated layers, the extension to curved manifolds generating the layers and to interior layers. Other important topics are the resolution of geometrical singularities and estimates of the remainder of the asymptotic expansion which will be discussed elsewhere.

References [I] Th. Apel, A note on anisotropic interpolation error estimates for isoparametric quadrilateral finite elements, SFB 393, Preprint 96-10, TU Chemnitz-Zwickau, 1996. [2] Th. Ape1 and G. Lube, Anisotropic mesh refinement for singularly perturbed reaction-diffusion problems, Appl. Numer. Math. 26 (1998) 415-433. [3] Th. Ape1 and Cl. Lube, Anisotropic mesh refinement in stabilized Galerkin methods, Numer. Math. 74 (1996) 262-282. [4] R. Codina, A discontinuity-capturing crosswind-dissipation for the finite element solution of the convection-diffusion equation, Comput. Methods Appl. Mech. Engrg. 110 (1993) 325-342. [5] R. Codina and J. Blasco, Analysis of a finite element approximation of the stationary Navier-Stokes equations using equal order velocity-pressure interpolation, Preprint CIMNE No. 113, Barcelona 1997 (submitted to SINUM). [6] M. Dobrowolski and H.-G. Roos, A-priori estimates for the solution of convection-diffusion problems and interpolation on Shishkin meshes, ZA 16 (1997) 1-13. [7] C. Farhat and L.P. Franca, Bubble functions prompt unusual stabilized finite element methods, Comput. Methods Appl. Mech. Engrg. 123 (1995) 299-308. [8] L.P. Franca and E.G. Dutra do Carmo, The Galerkin gradient least-squares method, Comput. Methods Appl. Mech. Engrg. 74 (1989) 41-54. [9] L.P. Franca, S.L. Frey and T.J.R. Hughes, Stabilized finite element methods: I. Application to the advective-diffusive model, Comput. Methods Appl. Mech. Engrg. 95 (1992) 253-276. [ 101 H. Goering, A. Felgenhauer, G. Lube H.-G. Roos and L. Tobiska, Singularly perturbed differential equations, Math. Research, Vol. 13 (Akademie-Verlag, Berlin, 1983). [l l] P. Grisvard, Elliptic Problems in Non-smooth Domains (Pitman, Boston, 1985). [ 121 H. Han and R.B. Kellog, Differentiability properties of solutions of the equation -e2 Au + TU=f in a square, SIAM J. Math. Anal. 21 ( 1990) 394-408. [ 131 C. Johnson, A.H. Schatz and L.B. Wahlbin, Crosswind smear and pointwise errors in streamline diffusion finite element methods, Math. Comput. 49 (1984) 179, 25-38. [14] R.B. Kellog, Boundary layers and corner singularities for a self-adjoint problem, in: M. Costabel et al., eds., Boundary value problems and integral equations in non-smooth domains (Dekker, New York 1995) 121-149.

[ 151 G. Lube, Stabilized Galerkm finite element methods for convection dominated and incompressible flow problems, Banach Center Public., vol. 29. Inst. Math., Polish Acad. SC., Warsrawa, 1994. [ 161 G. Lube and F.-C. Otto. A nonoverlapping domain decomposition method for incompressible flow problems, appears in: Proc. Third Int. Conf. Numeric. Modelling in Contin. Mech., Prague, 1997. [ 171 G. Lube and D. WeiM, Stabilired finite element methods for singularly perturbed parabolic problems. Appl. Numer. Math. 17 (1995) 43 I-459. 1181A. Quarteroni and A. Valli, Numerical Approximation of Partial Differential Equations (Sprmger, lYY4). [ 191 H.G. Roos, M. Stynes and I_. Tobiska, Numerical Methods for Singularly Perturbed Differential Equations (Springer, Heidelberg, 1996). [201 A.H. Schatz and L.B. Wahlbin, On the finite element method for singularly perturbed reaction-dlffuvion problems in two and one dimensions, Math. Comput. 40 (1983) 47-89. [2l] T. Skalicky and H.-G. Roes, Anisotropic mesh refinement for problems with internal and boundary layers. Preprint Math-NM-09-97, TU Dresden, Faculty of Math., 1997. [22] F. Valentin and L.P. Franca. Combining stabilized tinite element methods, UCDiCCM Report No. 48, Denver. 1995.