- Email: [email protected]
Residual-based stabilized higher-order FEM for advection-dominated problems
Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138 www.elsevier.com/locate/cma
Residual-based stabilized higher-order FEM for advection-dominated problems Gert Lube *, Gerd Rapin Mathematics Department, NAM, University of Go¨ttingen, Gottingen D-37083, Germany Received 4 January 2005; received in revised form 7 July 2005; accepted 26 July 2005
Abstract We reconsider the numerical solution of linear(ized) advection–diffusion–reaction problems using higher-order finite elements together with stabilized Galerkin methods of streamline-diffusion type (SUPG) and with shock-capturing stabilization. The analysis improves the a priori analysis in our previous paper [T. Knopp, G. Lube, G. Rapin, Stabilized finite element methods with shock capturing for advection–diffusion problems, Comput. Methods Appl. Mech. Engrg. 191 (2002) 2997–3013]. The theoretical results are supported by some numerical experiments. Ó 2005 Elsevier B.V. All rights reserved. Keywords: Advection–diffusion; Incompressible flow; Stabilized finite elements; SUPG
1. Introduction The motivation of the present paper stems, e.g., from the finite element simulation of the non-isothermal and incompressible Navier–Stokes problem ot~ u r ðmr~ uÞ þ ð~ u rÞ~ u þ rp ¼ bh~ g; r ~ u ¼ 0; ot h þ ð~ u rÞh r ðarhÞ ¼ q_ V =cp
ð1Þ ð2Þ ð3Þ
for velocity ~ u, pressure p and temperature h in a polyhedral domain X Rd, d 6 3, with source terms bh~ g and q_ V =cp . This model describes, e.g., the air flow in buildings, etc. [17]. The momentum and continuity Eqs. (1) and (2) describe the fluid motion; the heat transfer is driven by the advection–diffusion Eq. (3). Turbulence may occur at high Rayleigh or Reynolds numbers. A standard approach is to consider the Reynolds averaged Navier–Stokes equations (RANS) together with, e.g., the k– turbulence model. Within a statistical turbulence model only averaged values are considered. An eddy viscosity ansatz for turbulent effects is modeled as an additional diffusion term with eddy viscosity mt. Then the averaged values for ~ u, p and h are determined by (1)–(3) with m and a replaced with (variable) viscosities me = m + mt and ae = a + mt/Prt. The eddy viscosity term mt is determined, e.g. in the k– turbulence model, by mt = clk2/ where the turbulent kinetic energy k and the dissipation rate of k are defined by additional (nonlinear) advection–diffusion–reaction equations. A standard algorithmic treatment of the coupled model is to semi-discretize, in an outer loop, in time (with possible step control) using an A-stable method and then, in an inner loop, to decouple and linearize the resulting system. A block *
Corresponding author. E-mail addresses: [email protected] (G. Lube), [email protected] (G. Rapin).
0045-7825/$ - see front matter Ó 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2005.07.017
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138
4125
Gauss–Seidel method with fixed point or Newton-type iteration per time step leads to linearized Navier–Stokes problems (of Oseen-type) and linearized advection–diffusion–reaction problems as auxiliary problems, see [17]. A proper numerical approach to the latter model Lu :¼ r ðaruÞ þ ~ b ru þ cu ¼ f
ð4Þ
in X
is an important ingredient of the approach and will be discussed in this paper. The streamline upwind Petrov–Galerkin (SUPG) method, proposed by Brooks and Hughes [1], was the first variationally consistent, stable and accurate finite element model for advection-dominated problems. It initiated the development of stabilization techniques for advection-dominated and related problems. For an overview, see, e.g., [14]. A first relevant analysis of the SUPG method can be found in Johnson et al. [15] in case of regular solutions. The theory has been refined over the years in several directions. Here we only mention the analysis of a hp-version by Houston and Su¨li [12]. Nevertheless, for non-smooth solutions, localized oscillations of the SUPG solution may still exist in the neighborhood of steep gradients. As a remedy, discontinuity- or shock-capturing terms can be added to enhance the stability. Linear (but non-consistent) schemes for low-order elements are considered, e.g., in [16,22]. Mizukami/Hughes [13] introduced the first nonlinear discontinuity-capturing schemes DC1 and DC2. The idea was to enhance, additionally to streamline upwinding, numerical viscosity in the direction of $uh. The consistent approximate upwind (CAU) in [9] provided a further refinement. For recent developments of the CAU scheme to higher-order elements, we refer to [10]. Moreover, Codina considered in [5,6] the discontinuity-capturing/crosswind-dissipation (DC/CD) scheme with additional anisotropic viscosity. The important question of low-order (nonlinear) schemes which satisfy a discrete maximum principle is discussed in the recent papers [2,3], see also the monograph [8]. A first theoretical result for such nonlinear schemes is seemingly due to Szepessy [23]. In our previous paper [18] we considered the a priori analysis of a rather general class of shock-capturing schemes. The goal of the present paper is an extension of the stabilized higher-order FE method of the recent paper [12] to the case of shock-capturing stabilization. In particular, we address the choice of the stabilization parameters, see Section 2, extend and refine the analysis of shockcapturing schemes given in [18], see Section 3, and provide some numerical experiments, see Section 4. 2. Stabilized FEM for advection–diffusion–reaction model Following basically [12], we describe and analyze the SUPG-stabilization of the advection–diffusion–reaction model. In contrast to [12], we give a refined definition of the stabilization parameters depending on all critical parameters. 2.1. Problem statement For the advection–diffusion–reaction scheme (4), we assume a; c 2 L1 ðXÞ, ~ b 2 ðH 1 ðXÞÞ \ ðL1 ðXÞÞ , f 2 L2(X) and d
ðr ~ bÞðxÞ ¼ 0;
cðxÞ P x P 0;
aðxÞ P a0 > 0;
a.e. in X.
d
ð5Þ
For simplicity only, we analyze the homogeneous Dirichlet problem u¼0
ð6Þ
on oX.
The basic variational formulation of (4)–(6) reads: Find u 2 V :¼ H 10 ðXÞ
s.t. Aðu; vÞ ¼ lðvÞ
8v 2 V
ð7Þ
with Aðu; vÞ ¼ ðaru; rvÞX þ ð~ b ru þ cu; vÞX ; lðvÞ ¼ ðf ; vÞX .
ð8Þ ð9Þ
2.2. Finite element discretization Suppose a family of admissible triangulations Th ¼ fT g of the polyhedral domain X where h is the piecewise constant mesh function with h(x) = hT = diam (T), x 2 T. We assume that Th is shape-regular, i.e. there exists a constant Cr 5 Cr(h) such that C r hdT 6 meas ðT Þ
8T 2 [h Th .
ð10Þ
Moreover, we assume that each element T 2 Th is a smooth bijective image of a given reference element Tb , i.e., T ¼ F T ð Tb Þ for all T 2 Th . Here, Tb is the (open) unit simplex or the (open) unit hypercube in Rd.
4126
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138
For p 2 N, we denote by Pp ð Tb Þ the set of polynomials of degree at most p on Tb . To each T 2 Th an integer pT P 1 is assigned; then we define the vector p ¼ fpT : T 2 Th g. We introduce a conforming finite element space on Th by V ph ¼ X ph \ H 10 ðXÞ;
b X ph ¼ fv 2 CðXÞjvj T F T 2 P pT ð T Þ 8T 2 Th g.
ð11Þ
A key point in the subsequent analysis are the following local inverse inequalities 8w 2 X ph
krwkðL2 ðT ÞÞd 6 linv p2T h1 T kwkL2 ðT Þ
ð12Þ
on T 2 Th . The constant linv depends only on the shape-regularity parameter Cr in (10), see [21]. Interpolation estimates will be given later on. The standard Galerkin-FEM is, in general, not appropriate for the numerical solution of problem (7)–(9) in the singularly perturbed case. Therefore we consider the well-known streamline upwinding (SUPG)-stabilization: Find uh 2 V ph : As ðuh ; vÞ ¼ Ls ðvÞ with As ðu; vÞ :¼ Aðu; vÞ þ Ls ðvÞ :¼ lðvÞ þ
X
X
8v 2 V ph
ð13Þ
dT ð b Lu; ~ b rvÞT ;
ð14Þ
T 2Th
dT ðf ; ~ b rvÞT ;
ð15Þ
T 2Th
b LujT :¼ r PT ðaruÞ þ ~ b ru þ cu 2
ð16Þ
d
d
and the orthogonal projection PT : ½L ðT Þ ! ½P pT ðT Þ . This modification of the (variable) diffusion term allows the application of inverse inequalities [12]. Remark. Other residual-based stabilization variants are the Galerkin/least-squares (GLS) method and the algebraic subgrid scale (ASGS) method where the test function dT ð~ b rvÞ is replaced with dTLv and dT(L v), respectively. All these P 1=2 dT k~ b rðÞk2 . Our numerical methods combine stability and accuracy and allow additional error control of T 2Th
0;T
experience shows that the additional test terms of the GLS and ASGS methods are less important. For this reason, we restrict the subsequent analysis to the SUPG-case. 2.3. Basic quasi-optimal error estimate The following estimates of the bilinear form As, using the stabilized norm ! 1=2 X pffiffiffi p ffiffi ffi 2 2 2 jjjvjjj :¼ k arvk 2 d þ k cvk 2 þ dT k~ b rvk 2 ; L ðT Þ
ðL ðT ÞÞ
ð17Þ
L ðT Þ
T 2Th
are slight modifications of results given in [12,18] and will be important in Section 3, too. For the convenience of the reader, we include the proofs in the Appendix. First of all, we obtain the following stability and continuity results. Lemma 2.1. Let the assumptions (5) and the following condition be valid: ! 1 h2T 1 0 6 dT 6 min 4 2 ; . 2 pT linv kakL1 ðT Þ kckL1 ðT Þ
ð18Þ
Then the bilinear form As in (14) satisfies the stability estimate 1 2 As ðv; vÞ P jjjvjjj 2
8v 2 V h .
ð19Þ
Moreover, for u 2 H 10 ðXÞ with ($ Æ (a$u))kT 2 L2(T) and vh 2 V ph , we obtain As ðu; vh Þ 6 Qs ðuÞjjjvh jjj Qs ðuÞ :¼ jjjujjj þ
X T 2Th
with aT := infx2Ta(x).
(
) !1=2 2 bkðL1 ðT ÞÞd 1 k~ 2 min ; þ kukL2 ðT Þ dT aT
X T 2Th
!1=2 dT k r PT ðaruÞ þ
2 cukL2 ðT Þ
ð20Þ
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138
4127
We combine both results to prove an a priori estimate for the SUPG scheme. Lemma 2.2. Let u 2 V be a solution of (7)–(9) with ($ Æ (a$u))kT 2 L2(T) for all T 2 Th and let the assumptions (5) and (18) be valid. This implies the quasi-optimal a priori error estimate of the SUPG scheme (13)–(16): !1=2 X 2 b LÞðuÞk 2 jjju uh jjj 6 C 1 inf Qs ðu vh Þ þ C 2 dT kð L ð21Þ L ðT Þ
p
vh 2V h
T 2Th
with data-independent constants Ci. In order to determine the stabilization parameter dT we prove an error estimate using appropriate approximation properties. Let be v 2 Hk(T) with k P r > d2. The standard Lagrange interpolation I Th;p v satisfies kv I Th;p vkH s ðT Þ 6 C I
hls T kvkH k ðT Þ ; pks T
0 6 s 6 l ¼ minðpT þ 1; kÞ.
ð22Þ
But we need a global interpolation operator. Unfortunately we cannot use the global Lagrange interpolation Ih,p for nonconstant pT, since then the Lagrange interpolation is in general not continuous and therefore it is not included in V ph . Nevertheless, in order to avoid additional technical problems, we require that I h;p v 2 V ph is satisfied, which is valid if p = pT for all T 2 Th . Then the global Lagrange interpolation Ih,pv for all v 2 H 10 ðXÞ \ H r ðXÞ; r > d2 with vjT 2 Hk(T), k P r, satisfies (22) with s 2 {0, 1} for all T 2 Th . Remark (i) The assumption p = pT on Th is quite restrictive. A remedy is a non-conforming domain decomposition method where a constant polynomial degree is assumed in each subdomain. The polynomial degree can differ from subdomain to subdomain, cf. [24]. (ii) The assumption r > d2 can be relaxed to r > 12 if the Scott–Zhang quasi-interpolant is used. But the right-hand side norm in (22) has to be replaced by the Hk-norm over a patch of elements containing T . The following result is a slight modification of Theorem 9 in [12]. Corollary 2.1. Let the assumptions (5), (18) be satisfied. Assume that the weak solution u 2 H 10 ðXÞ of (7)–(9) belongs locally to H kT ðT Þ, k T > d2 and a 2 W kT 1;1 ðT Þ for all T 2 Th . Moreover, we assume that I h;p u 2 V ph ðXÞ. Then the following error estimate for the solution uh of the SUPG scheme (13)–(16) is valid: X h2ðlT 1Þ T
2
jjju uh jjj 6 C
T 2Th
þ
2
2
kukH kT ðT Þ kakL1 ðT Þ þ 2ðk T 1Þ
kakW kT 1;1 ðT Þ
pT
kckL1 ðT Þ h2T p2T
kakL1 ðT Þ
(
þ
dT k~ bk2ðL1 ðT ÞÞd
2 bkðL1 ðT ÞÞd 1 k~ þ min ; dT aT
)
! h2T . p2T
ð23Þ
2.4. Parameter design Up to now, we have not specified the stabilization parameter dT. Here we revisit the analysis already considered in [12]. Corollary 2.2. Let the assumptions of Corollary 2.1 be valid. The choice of the stabilization parameters according to ( ) hT 1 h2T dT min ; ; ð24Þ kckL1 ðT Þ p4T l2inv kakL1 ðT Þ p k~ bk 1 T
L ðT Þ
is consistent with condition (18). Moreover, it leads to X h2ðlT 1Þ opt 2 2 T M T kukH kT ðT Þ jjju uh jjj 6 C 2ðk T 1Þ T 2Th pT
ð25Þ
with
2
M opt T
:¼ kakL1 ðT Þ 1 þ
kakW kT 1;1 ðT Þ 2
kakL1 ðT Þ
þ PeT þ CT þ min
kakL1 ðT Þ aT
Pe2T ; p2T l2inv
! ð26Þ
4128
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138
and PeT :¼
hT k~ bkðL1 ðT ÞÞd pT kakL1 ðT Þ
CT :¼
;
kckL1 ðT Þ h2T p2T kakL1 ðT Þ
ð27Þ
.
Proof. The idea is to equilibrate the dT-dependent right-hand side terms of the error estimate (23) under the side conditions (18) and continuity of dT. 1 dT
We consider the case are equilibrated if dT
k~ bk2
6
ðL1 ðT ÞÞd
aT
in the min-term of the right-hand side of (23). Then the dT -dependent terms in (23)
hT . ~ pT kbkðL1 ðT ÞÞd
ð28Þ
In order to fulfil the condition (18), we set dT according to (24). This leads to ( ) ( ( ) 2 ) 2 2 bkðL1 ðT ÞÞd h2T hT k~ bkðL1 ðT ÞÞd hT k~ bkðL1 ðT ÞÞd kckL1 ðT Þ h2T 2 2 1 k~ Z T :¼ min ; ¼ min max ; ; pT linv kakL1 ðT Þ ; dT aT p2T pT p2T p2T aT kakL1 ðT Þ 2 2 2 ¼ kakL1 ðT Þ min PeT ; max PeT ; CT ; pT linv ; aT hence M opt T :¼ kakL1 ðT Þ þ 6 kakL1 ðT Þ 6 kakL1 ðT Þ
kak2W kT 1;1 ðT Þ kakL1 ðT Þ
þ
kckL1 ðT Þ h2T p2T
2
þ dT k~ bkðL1 ðT ÞÞd þ Z T
! ZT 1þ þ PeT þ CT þ kakL1 ðT Þ kak2L1 ðT Þ ! 2 kakW kT 1;1 ðT Þ kakL1 ðT Þ 2 2 2 1þ þ PeT þ CT þ min PeT ; pT linv . 2 aT kakL1 ðT Þ kak2W kT 1;1 ðT Þ
This leads to the estimate (25)–(27).
h
Remark. Let us discuss the optimality of the result. For brevity, we consider the case a(x) a. The estimate (25)–(27) gives in the advection/reaction-dominated limit maxfPeT ; CT g P p2T l2inv the estimates X
hT
T 2Th
pT k~ bkðL1 ðT ÞÞd
X T 2Th
kð~ b rÞðu uh Þk2L2 ðT Þ 6 C
X h2lT 1 kuk2H kT ðT Þ ; 2k T 1 p T 2Th
X h2lT pffiffiffi 1 2 2 k cðu uh ÞkL2 ðT Þ 6 C kukH kT ðT Þ ; 2k T kckL1 ðT Þ p T 2Th
PeT P CT ;
CT P PeT .
These estimates are simultaneously optimal in h and p with respect to the (weighted) streamline-derivative in the advectiondominated limit PeT P CT P p2T l2inv , and in the (weighted) L2-norm in the reaction-dominated limit CT P PeT P p2T l2inv , respectively, see also [12], Remark 10a. Please note that both estimates are also valid for the complete error term jjj Æ jjj2. In the diffusion-dominated limit with PeT 6 1 ð6 p2T l2inv Þ, we obtain for the diffusion-relevant part of the jjj Æ jjj-norm 2
akrðu uh ÞkL2 ðXÞ 6 C
X T 2Th
a
h2ðlT 1Þ 2 kukH kT ðT Þ ; p2ðkT 1Þ
PeT 6 1
which is simultaneously optimal in h and p with respect to the H1-seminorm, see also [12], Remark 10a. Nevertheless, there remains a transition region between the diffusion-dominated limit and the advection/reaction-dominated limit with loss of optimality with respect to p by p1=2 (at least using the present technique). We want to point out that a similar analysis has been given in Remark 10 of [12], but the analysis is seemingly not complete in the transition region. Remark. We would like to point out that the dependence of the parameters dT on the spectral order p is specified. For fixed hT, it tends to zero with increasing p. This is in (formal) agreement with the behaviour of the artificial viscosity in the spectral vanishing viscosity method, see, e.g., [11].
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138
4129
3. Shock-capturing stabilization The SUPG method stabilizes most of the unphysical oscillations of the Galerkin discretization which are caused by ? dominating advection. Nevertheless, there may remain spurious localized oscillations, often in crosswind-directions ~ b . Such perturbations are very dangerous in the numerical solution of coupled nonlinear problems, see Section 1. A wellknown remedy is the shock-capturing variant of the SUPG method. For a review of such methods we refer to [22]. Versions which add linear crosswind-stabilization terms can be constructed for piecewise linear and bilinear elements. Energy estimates and even maximum norm estimates can be found e.g. in [16,22,26]. In practice, these methods suffer from smearing in crosswind-directions since, due to the linearity of the schemes, the stabilization effect cannot be restricted to certain subregions (e.g. close to shocks). An interesting new (but nonlinear) scheme in [2,3] for piecewise linear elements allows to satisfy a discrete maximum principle. 3.1. Framework of shock-capturing The extension of shock-capturing schemes to higher order elements requires a nonlinear construction. Different proposals in the literature [5,9,13,22] are covered by the following class of schemes [18]: Find U h 2 V ph : As ðU h ; vÞ þ Asc ðU h ; U h ; vÞ ¼ Ls ðvÞ with Asc ðw; u; vÞ ¼
X
8v 2 V ph
ðsT ðwÞDsc ru; rvÞT ;
ð29Þ
ð30Þ
T 2Th
where Dsc denotes a symmetric, positive semi-definite matrix function with kDsc kL1ðXÞd d 6 1. The non-negative limiter function sT(w) has to restrict the effect of shock-capturing to subregions where the (modified) residual b Lw f is too large. The common feature of different proposals is that sT is a function of the scaled (modified) residual, namely sT ðwÞ ¼ sT ðRT ðwÞÞ;
RT ðwÞ :¼
kb Lw f kL2 ðT Þ kwkH 1 ðT Þ þ jT
ð31Þ
;
where we introduced a regularization parameter jT > 0. Standard shock-capturing schemes [5,9,13,22] with jT = 0 lead to ill-posed nonlinear problems. For a possible choice of jT, see Example 3.1. Some results for the h-version of this scheme can be found in [18]. Here we refine and extend the analysis to the hp-version of the shock-capturing scheme (29) and (30). Lemma 2.1 allows to prove the solvability of the nonlinear scheme (29)– (31). Unfortunately, a (local) uniqueness result is still open. Theorem 3.1. Let the assumptions (5) and (18) be valid. Moreover, assume that the limiter function sT(w) in (30) is continuous w.r.t. w. The SUPG-scheme with shock-capturing stabilization (29) and (30) has at least one solution U h 2 V ph . The following a priori estimate is valid for each solution 2 X pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1=2 2 ð32Þ jjjU h jjj2 þ sT ðU h ÞDsc rU h 2 6 Cjjjf jjj L ðT Þ
T 2Th
with the dual norm jjjf jjj :¼ supv2V p ;jjjvjjj¼1 ðf ; vÞX . h
Proof. The proof relies on a variant of Brouwers fixed point theorem where Lemma 2.1 is used; for details see [18], Theorem 3.5. h 3.2. Shock-capturing with isotropic diffusion A priori error estimates of scheme (29) and (30) require more information about the shock-capturing term Asc. Let us first consider the case of isotropic diffusion with almost quadratical dependence of the function sT(w) on the scaled residual RT ðwÞ according to 2
Diso sc ¼ I;
2 iso siso T ðwÞ :¼ lT ðwÞ½RT ðwÞ ¼
liso T ðwÞ½RT ðwÞ
½kwkH 1 ðT Þ þ jT 2
ð33Þ
Lwh f kL2 ðT Þ . The class of schemes (29) and (30) with (33) covers, in particular, the consistent approximate with RT ðwh Þ ¼ k b upwind (CAU) scheme of [9] with jT = 0; it is also closely related to the original discontinuity-capturing schemes DC1 and
4130
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138
DC2 of [13]. The basic idea of these schemes was to introduce, additionally to streamline upwinding, numerical viscosity in the direction of $Uh. The following global a priori estimate shows that the convergence properties of the method are at least as good as for the SUPG-FEM. Moreover, we have stronger control of the gradient, including the crosswind-direction(s). Let again vh ¼ I h;p u 2 V ph be an appropriate interpolant of the solution u 2 V of (7)–(9). In the forthcoming analysis we introduce the quantities iso iso F iso h ¼ F h ðI h;p u; U h Þ :¼ sup F T ðI h;p u; U h Þ; T
F iso T ðvh ; U h Þ :¼
jvh jH 1 ðT Þ jU h jH 1 ðT Þ þ jT
.
ð34Þ
Assume, for simplicity, that Ih,pu can be taken as the Lagrange interpolant of the (sufficiently smooth) solution u. Then it is shown in [18], Remark 4.6 that the constant F iso h is uniformly bounded w.r.t. the mesh size h and viscosity a. See also next remark. p Theorem 3.2. Let the assumptions (5), (18), and (33) with 0 6 liso T ðwÞ 6 qdT for all w 2 V h and a sufficiently small q > 0 be 2 satisfied. Assume that the weak solution u 2 V of (7)–(9) satisfies locally ($ Æ (a$u))jT 2 L (T) and belongs to H kT ðT Þ; k T > d2 for all T 2 Th . Moreover, we assume that vh :¼ I h;p u 2 V ph . Then the following error estimate for the solution Uh of the stabilized discrete problem (29), (30) and (33) is valid: X X
2 2 2 2 siso þ K2 dT kð b L LÞðuÞkL2 ðT Þ ð35Þ jjju U h jjj þ T ðU h Þju U h jH 1 ðT Þ 6 K 1 Qs ðu I h;p uÞ T 2Th
T 2Th
with data-independent constants Ki. Proof. We set e := Uh vh. and g := u vh. Lemma 2.1 and the definition of Aiso sc imply 1 iso iso iso jjjejjj2 þ Aiso sc ðU h ; e; eÞ 6 As ðe; eÞ þ Asc ðU h ; e; eÞ ¼ As ðU h ; eÞ þ Asc ðU h ; U h ; eÞ As ðvh ; eÞ Asc ðU h ; vh ; eÞ 2 X ¼ As ðg; eÞ Aiso dT ðð b L LÞu; ~ b reÞT I þ II þ III. sc ðU h ; vh ; eÞ
ð36Þ
T 2Th
We estimate the right-hand side terms separately. Lemma 2.1 implies 1 2 2 I ¼ As ðg; eÞ 6 jjjejjj þ 2½Qs ðgÞ . 8
ð37Þ
The estimate of term II takes advantage of the structure (33) of siso T ðwÞ and of the terms defined in (34): X
II 6
T
6
!1=2 2 siso T ðU h ÞjejH 1 ðT Þ
X
!1=2 2 siso T ðU h Þjvh jH 1 ðT Þ
T
3 X iso 1 X iso 2 2 2 sT ðU h ÞjejH 1 ðT Þ þ s ðU h Þ½F iso T ðjU h jH 1 ðT Þ þ jT Þ 4 T 3 T T
3 1 X iso 2 2 6 Aiso l ðU h Þ½F iso sc ðU h ; e; eÞ þ T RT ðU h Þ 4 3 T T X 3 1 iso 2 6 Aiso dT R2T ðU h Þ sc ðU h ; e; eÞ þ q½F h 4 3 T X 3 iso 2 2 2 2 b b b 6 Asc ðU h ; e; eÞ þ q½F iso d k Lek þ k Lgk þ kð L LÞuk 2 2 2 T L ðT Þ L ðT Þ L ðT Þ . h 4 T 2Th
ð38Þ
In the last substep we used the Young inequality together with LU h f kL2 ðT Þ 6 k b LðU h vh ÞkL2 ðT Þ þ k b Lðvh uÞkL2 ðT Þ þ kð b L LÞukL2 ðT Þ . RT ðU h Þ ¼ k b For the last right-hand side term in (36) we find III 6
X 1 X 2 2 dT k~ b rekL2 ðT Þ þ 2 dT kð b L LÞukL2 ðT Þ . 8 T 2Th T 2Th
ð39Þ
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138
4131
Summarizing (36)–(39) we obtain X X 1 1 2 2 2 2 2 iso 2 iso 2 b b ½ jjjejjj þ Aiso ðU ; e; eÞ 6 2 Q ðgÞ þ q½F d k Lek þ k Lgk dT kð b L LÞukL2 ðT Þ . 2 2 h T s L ðT Þ L ðT Þ þ ðq½F h þ 2Þ h sc 4 4 T 2Th T 2Th ð40Þ Inverse inequalities (12) and assumption (18) imply X pffiffiffi 2 1 X pffiffiffi 2 2 dT k r PT ðareÞ þ cekL2 ðT Þ 6 k arvkðL2 ðT ÞÞd þ k cvkL2 ðT Þ . 2 T 2Th T 2Th 2
1 Using the triangle inequality together with the latter estimate, we find with q½F iso h 6 16 that X X 2 2 2 2 2 q½F iso dT k b LekL2 ðT Þ 6 2q½F iso dT k~ b rekL2 ðT Þ þ kr PT ðareÞ cekL2 ðT Þ h h T 2Th
T 2Th
1 2 2 2 6 2q½F iso h jjjejjj 6 jjjejjj . 8
ð41Þ
Please note that we take advantage of the assumption that vh = Ih,pu is the Lagrange interpolant of u 2 V in order to have a proper bound of F iso h . Furthermore, we obtain via triangle inequality and the definition of Qs(g) that 1 X X 2 2 2 2 2 iso 2 ~ b q½F iso d k Lgk 6 2q½F d k b rgk þ kr P ðargÞ cgk ð42Þ 2 2 2 T T T L ðT Þ L ðT Þ L ðT Þ 6 ½Qs ðgÞ . h h 8 T 2Th T 2Th We collect (40)–(42) to X 1 2 2 2 ~ ~ jjjejjj þ Aiso dT kð b L LÞukL2 ðT Þ sc ðU h ; e; eÞ 6 K 1 ½Qs ðu vh Þ þ K 2 2 T 2Th ~ i . The triangle inequality concludes the proof, using that the term Aiso ðU h ; g; gÞ can be estimated with appropriate positive K SC similarly as the term II. h Remark. In order to relax the regularity assumptions on u 2 V in Theorem 3.2, let vh ¼ I h;p u 2 V ph be chosen as a quasiinterpolant of u 2 V with the local stability estimate jvh jH 1 ðT Þ 6 CjujH 1 ðV ðT ÞÞ on a suitable neighborhood of V(T) of element T. Then F iso h is uniformly bounded w.r.t. h, see [18], Remark 4.6. The present analysis yields so far only an upper bound for the limiter function liso T ðU h Þ 6 qdT ; it gives no hint to its explicit construction. Example 3.1. For the special case of a = const. and of piecewise linear elements (i.e. p = 1), the CAU-scheme in [9] is determined by ! ~
hT j~ bk j d0 hT cau lT ðU h Þ ¼ max 0; dk dT ; dk ¼ min 1; ; 6a j~ bk j ~ brU h rU h if j$Uhj 5 0 and ~ bk ¼ 0 else. Numerical results, see, e.g. [7,9,13], show that this scheme is suitable for where ~ bk ¼ jrU j2 h
solutions with strong gradients in layers. For problems with smooth solutions, however, it suffers from smearing in crosswind-direction(s). An improved version of the CAU scheme, even for higher-order elements, has been discussed in the recent paper [10]. The analysis (with some open ends) is based on an linearized iterative scheme. The present analysis allows to set lcau T 6 qdT . Other analytical tools are necessary to close the gap between theoretical and practical results. 3.3. Shock-capturing with crosswind-diffusion The crosswind-smearing of the DC- and CAU-schemes lead Codina in [5,6] to the construction of a shock-capturing scheme with anisotropic diffusion in crosswind-direction(s). This discontinuity-capturing crosswind-dissipation (DC/CD) scheme can be seen as a special case of the general schemes (29) and (30) with crosswind-diffusion and almost linear dependence on RT ðwÞ: 8 ~ ~ > jwj : H 1 ðT Þ þ jT ~ 0; b¼0
4132
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138
Codinas scheme is introduced in [6] with jT = 0 and the limiter function 2kakL1 ðT Þ dc=cd cd lT ðwÞ ¼ lT ðwÞ :¼ l0 hT max 0; b hT RT ðwÞ
ð44Þ
h R ðwÞ
T T can be seen as a pseudo mesh Peclet number. Numerical experiments show that b is a critical with l0 ¼ 12. The term 2kak L1 ðT Þ parameter. Codina proposed b = 0.7 for p = 1 and b = 0.35 for p = 2. Here we try to mimic the proof of Theorem 3.2. The following result shows that this is possible under the following restriction on the limiter function
0 6 lcd T ðwÞ 6 qdT RT ðwÞ;
w 2 V ph
ð45Þ
with a suitable parameter q > 0. This means that both variants of the shock-capturing scheme with isotropic and with anisotropic diffusion have the same upper bounds of the additional shock-capturing term. Theorem 3.3. Let the assumptions of Theorem 3.2 be valid. Then we obtain for the shock-capturing scheme (29), (30) and (43) with assumption (45) the estimate: X X
2 2 2 2 cd 1=2 scd rðu U h ÞkL2 ðT Þ 6 K 3 Qs ðu I h;p uÞ þ K 4 dT kð b L LÞukL2 ðT Þ ð46Þ jjju U h jjj þ T ðU h Þk½Dsc T 2Th
T 2Th
with data-independent constants Ki. Proof. We modify the proof of Theorem 3.2 and introduce the quantities cd F cd h ðI h;p u; U h Þ :¼ sup F T ðI h;p u; U h Þ; T
F cd T ðI h;p u; U h Þ :¼
1=2 k½Dcd rI h;p ukL2 ðT Þ sc
jU h jH 1 ðT Þ þ jT
ð47Þ
.
The constant F cd h is uniformly bounded w.r.t. h and a. We start again from (36), but the estimate of term II is different depending on the structure of scd T ðwÞ. Proceeding similarly as in (38) and using (47), (43) and (45) we obtain
II 6
X
!1=2 2 cd 1=2 scd rI h;p ukL2 ðT Þ T ðU h Þk½Dsc
X
T
6
!1=2 2 cd 1=2 scd rekL2 ðT Þ T ðU h Þk½Dsc
T
3 X cd 1 X cd 1=2 2 2 2 sT ðU h Þk½Dcd rekL2 ðT Þ þ s ðU h Þ½F cd sc T ðjU h jH 1 ðT Þ þ jT Þ 4 T 3 T T
3 1 X cd 2 6 Acd l ðU h Þ½F cd sc ðU h ; e; eÞ þ T RT ðU h ÞðjU h jH 1 ðT Þ þ jT Þ 4 3 T T X 3 1 cd 2 6 Asc dT R2T ðU h Þ. sc ðU h ; e; eÞ þ q½F h 4 3 T Assuming a sufficiently small q, we proceed as in the proof of Theorem 3.2.
ð48Þ h
As in the previous case (cf. discussion of Theorems 3.2), 3.3 provides only an upper bound of the limiter function lcd T ðwÞ 6 qdT RT ðwÞ. As already discussed for the CAU-scheme, other analytical tools are necessary to refine the design of the limiter functions lT(Æ). p Remark. We propose to choose the parameter jT in such a way that lcd T ðwÞ remains bounded. For all w 2 V h , the inverse estimate (12) leads to
RT ðwÞ ¼ 6
k r PT ðarwÞ þ ~ b rw þ cw f kL2 ðT Þ jwjH 1 ðT Þ þ jT linv p2T hT
1 kakL1 ðT Þ þ k~ bkðL1 ðT ÞÞd þ kckL1 ðT Þ kwkL2 ðT Þ þ kf kL2 ðT Þ . jT
Choosing jT ¼ kckL1 ðT Þ kwkL2 ðT Þ þ kf kL2 ðT Þ
ð49Þ
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138
we obtain
(
dT RT ðwÞ
1 h2 ; ; 4 2 T kckL1 ðT Þ pT linv kakL1 ðT Þ L1 ðT Þ
hT 6 min pT k~ bk
)
4133
linv p2T kakL1 ðT Þ þ k~ bkðL1 ðT ÞÞd þ 1 . hT
This means that the artificial viscosity in Acd sc is at most of order hT/pT. 4. Numerical experiments Let us verify the error estimates of Sections 2 and 3, starting with a problem with a smooth solution. Example 4.1. Consider problem (4)–(6) with ~ b ¼ ðx2 ; x1 ÞT , a 2 {103,106}, and c 2 {0,103} in the domain X = (0, 1)2. The source term f is chosen in such a way that uðxÞ ¼ sinðpx1 Þ sinðpx2 Þex1 x2 is the (smooth) exact solution of (4)–(6). The numerical solution of the SUPG-scheme (13)–(16) with the parameter choice according to (24) is calculated with FEMLAB 2.3 on quasi-uniform meshes for Pp -elements with p 2 {1, 2, . . ., 6}. The additional shock-capturing stabilization is not required here. pffiffiffi 2 In Fig. 1 we present the p-convergence of the numerical error measured in the energy norm jjjvjjjE :¼ ðk arvkL2 ðXÞ þ 2 1=2 b cover all cases kvkL2 ðXÞ Þ for different values of h. The chosen values of c, a and h, p together with the variable flow field ~ of the parameter design for dT according to (24). The convergence rates, including the h-convergence for fixed p (not shown here), confirm the theoretical result of Corollary 2.2. In general, we observed for problems with a smooth solution a remarkable robustness of the results with respect to the design of the stabilization parameters. This is a consequence of the residual-based stabilization terms. We consider now the shock-capturing method (29)–(31) in case of rough solutions where typically unphysical oscillations appear around layers of the continuous solution. Interior error estimates, i.e. in subdomains X 0 X away from 0
0
10
10 h=1/4 h=1/8 h=1/16 h=1/32
–2
10
10
–4
–4
10
energy
energy
10
–6
10
–8
–6
10
–8
10
10
–10
–10
10
10
–12
10
h=1/4 h=1/8 h=1/16 h=1/32
–2
1
–12
2
3
p
4
5
10
6
1
2
3
0
5
6
0
10
10 h=1/4 h=1/8 h=1/16 h=1/32
–2
10
h=1/4 h=1/8 h=1/16 h=1/32
–2
10
–4
–4
10
10
–6
energy
energy
4
—
—
10
–8
10
–10
–8
10
10
–12
–12
10
10
–14
1
–6
10
–10
10
10
p
–14
2
3
p
—
4
5
6
10
1
2
3
p
4
5
6
—
Fig. 1. p-convergence for Example 4.1 with smooth solution and different values of a 2 {103, 106}, c 2 {0, 103}, and h 2 f14 ; 18 ; 161 ; 321 g.
4134
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138
boundary and interior layers, were given, e.g., by [20,25]. It is well-known that the nonlinear discrete problem (29)–(31) is very ill-posed for jT = 0, cf. [2,22]. Here we solve the nonlinear problem via simple iteration n nþ1 s n 2 N0 : As ðU nþ1 h ; vÞ þ Asc ðU h ; U h ; vÞ ¼ L ðvÞ
8v 2 V ph .
ð50Þ
Using the initial guess U 0h ¼ 0, the first iterate is given by the solution of the SUPG scheme (13)–(16). For a related approach for the CAU scheme, we refer to [10]. The iteration (50) with the rather strong stopping criterion kU hnþ1 U nh k1 6 104 is robust and fast for p = 1. For p P 2 we sometimes observed that the simple iteration did not converge. Then the application of Newton-type methods is recommended, see, e.g., [22]. Example 4.2. Consider problem (4) and (6) in the domain X = (0, 1)2 with a = 106, ~ bðxÞ ¼ p1ffiffi5 ð1; 2ÞT , c = 0. The source 1 2x x 2 4 ffi Þ is the exact solution. It is characterized by a characteristic term f is chosen in such a way that uðxÞ ¼ 12 ð1 tanh 1pffiffiffi 5a pffiffiffi interior layer of thickness Oð aj ln ajÞ around 2x1 x2 ¼ 14. We consider the solutions of the stabilized discrete scheme (29)–(31) without and with shock-capturing with different Ppelements for p 2 {1, 2, 3, 4} on a quasi-uniform mesh with h ¼ 641 . More precisely, we applied the DC/CD variant (44) with b = 0.7, j = 104, l0 = 1 for p 5 1 and l0 = 2 for p = 2. We stopped the simple iteration scheme (50) after 10 iterations if it had not reached the stopping criterion. In Table 1 we compare the error in different norms for the SUPG and the DC/CD scheme for Example 4.2 for p 2 {1, . . ., 4}. First of all, we observe the error reduction with increasing p. Although the SUPG scheme always shows a slightly smaller error, one observes that the errors are of the same magnitude. The larger errors of the shock-capturing scheme can be explained by the slightly increased crosswind-diffusion. We want to study the crosswind-diffusion effects more precisely. In Fig. 2 the crosswind-diffusion parameter s is presented. One can observe very well, that the additional diffusion is located around the layer. As expected, away from the layer no additional crosswind-diffusion is added. Moreover, we present in Fig. 3 contour-plots of the SUPG- and DC/ CD-solutions for p = 1 and p = 4. Indeed we observe a slight additional crosswind-effect of the shock-capturing. Finally, in Fig. 4 we plot the cross-section of the SUPG- and of the DC/CD-solutions in the crosswind-direction at x1 + 2x2 = 1. We observe the classical Gibbs phenomenon of the SUPG solution in the neighborhood of the boundary layer, i.e. there appear significant over- and undershoots of the discrete solution. The oscillations of the SUPG scheme in the characteristic layer are clearly damped with shock-capturing. It is important to observe that the strong gradient
Table 1 The error for the SUPG scheme without and with shock capturing for Example 4.2 with a = 106, c = 0, j = 104, b = 0.7, l0 = 1, h = 1/64 in different norms p
1 2 3 4
k kL2 ðXÞ
jjj Æ jjjE
k kL1 ðXÞ
SUPG
SC/CD
SUPG
SC/CD
SUPG
SC/CD
0.0559 0.0310 0.0216 0.0160
0.0587 0.0315 0.0234 0.0181
0.0589 0.0354 0.0266 0.0211
0.0615 0.0358 0.0280 0.0227
0.4421 0.3734 0.3098 0.2582
0.4985 0.3773 0.3515 0.3174
Fig. 2. The crosswind-diffusion parameter s for Example 4.2 with a = 106, c = 0, j = 104, b = 0.7, l0 = 1, h = 1/64, p = {1, 4}.
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138 SUPG
SC
1
1
0.9
0.9
0.8
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0 0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0 0.1
0.2
0.3
0.4
SUPG 1
0.9
0.9
0.8
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.5
0.6
0.7
0.8
SC
1
0 0.1
4135
0.5
0.6
6
0.7
0.8
0 0.1
0.2
0.3
0.4
4
Fig. 3. Contour plot for Example 4.2 with a = 10 , c = 0, j = 10 , b = 0.7, l0 = 1, h = 1/64, p = {1, 4}. On the left-hand side without shock capturing and on the right-hand side with shock capturing.
of the SUPG solution in the layer is still preserved by the shock-capturing scheme. Although not reflected in the present analysis of Section 3, this seems to be a result of the proper construction of the limiter function of the DC/CD-scheme with (44). Moreover, we observe that, for fixed h, the resolution of the strong gradient improves with increasing p. 5. Summary In this paper we extended the a priori analysis in [12] of a hp-version of stabilized Galerkin methods of streamline-diffusion type (SUPG) for diffusion–advection–reaction problems to the case of shock-capturing stabilization. This additional stabilization in crosswind-directions is an important ingredient within the solution of coupled nonlinear flow models, e.g., for turbulence models. In particular, we extended our previous results in [18] by studying carefully the dependence of stabilization parameters on the spectral order and on characteristic local quantities (like the mesh Peclet number). The analysis is supported by some numerical experiments. We improved our previous results in [18] for certain shock-capturing variants of the SUPG-scheme. In particular, we give better upper bounds of the limiter function which are in better agreement with schemes presented in the literature. Nevertheless, there are still open problems: Besides the uniqueness problem, a refined control of the crosswind-error (via the limiter function) would be desirable. Moreover, the convergence of iteration schemes for the arising nonlinear discrete problem is open. It seems to be useful to provide a careful and fair comparison of the present residual-based stabilization approach to the concept of edge-stabilization [2,4]. Moreover, it is necessary to extend the present concept to discontinuous Galerkin methods (eventually patchwise) in order to simplify the application of hp-adaptivity and domain decomposition techniques. In paper [19], we extend the present approach to the linearized incompressible Navier–Stokes problem, thus completing the analysis of auxiliary problems appearing within the iterative solution of coupled problems as shown in Section 1.
4136
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138 1.2
SC SUPG
1
1.2 1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0
–0.2 0
0.05
0.1
0.15
0.2
0.25
1.2
–0.2
0
0.05
0.1
0.15
0.2
0.25
1.2
SC SUPG
1
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0
0.05
0.1
0.15
0.2
0.25
SC SUPG
1
0.8
–0.2 0
SC SUPG
–0.2 0
0.05
0.1
0.15
0.2
0.25
6
Fig. 4. SUPG-FEM without and with SC-stabilization for Example 4.2 with a = 10 , c = 0 and h ¼ 641 ; p 2 f1; 2; 3; 4g.
Appendix Here we summarize some basic estimates for the SUPG-scheme applied to the linear advection–diffusion–reaction problems, see Section 2. Proof of Lemma 2.1. Setting u = v in (14) and integrating the advective term by parts together with (5), we obtain X pffiffiffi pffiffiffi 2 2 2 ðk arvkðL2 ðT ÞÞd þ k cvkL2 ðT Þ þ dT k~ b rvkL2 ðT Þ þ dT ðr PT ðarvÞ þ cv; ~ b rvÞT Þ. As ðv; vÞ ¼ T 2Th
We need the following auxiliary estimates, which use the local inverse estimates (12) together with the assumptions (18): pffiffiffi l p2 1=2 dT ðr PT ðarvÞ; ~ b rvÞT 6 dT inv T kakL1 ðT Þ k arvkðL2 ðT ÞÞd k~ b rvkL2 ðT Þ hT pffiffiffiffiffi 1 pffiffiffi b rvkL2 ðT Þ ; 6 pffiffiffi k arvkðL2 ðT ÞÞd dT k~ 2 pffiffiffi 1=2 b rvÞT 6 dT kckL1 ðT Þ k cvkL2 ðT Þ k~ b rvkL2 ðT Þ dT ðcv; ~ p ffiffiffiffiffi 1 pffiffiffi 6 pffiffiffi k cvkL2 ðT Þ dT k~ b rvkL2 ðT Þ . 2 We combine both estimates by using Youngs inequality pffiffiffi pffiffiffi 2 1 X ~ 1 2 2 2 2 As ðv; vÞ P jjjvjjj dT kb rvkL2 ðT Þ þ k arvkðL2 ðT ÞÞd þ k cvkL2 ðT Þ P jjjvjjj . 2 T 2Th 2
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138
4137
We now prove the continuity result: The symmetric part of As(u, vh) is bounded by jjjujjj jjjvhjjj. The remaining terms can be bounded as follows: ð~ b ru; vh Þ ¼ ðu; ~ b rvh Þ ( ) !1=2 2 X bkðL1 ðT ÞÞd 1 k~ 2 6 min ; kukL2 ðT Þ dT aT T 2Th X
dT ðr PT ðaruÞ þ cu; ~ b rvh ÞT 6
T 2Th
X
! 1=2 X pffiffiffi 2 2 k arvh kL2 ðT Þ þ dT k~ b rvh kL2 ðT Þ ; T 2Th
!1=2 dT k r PT ðaruÞ þ
cuk2L2 ðT Þ
T 2Th
X
!1=2 dT k~ b rvh k2L2 ðT Þ
.
T 2Th
Summarizing all terms and dividing by jjjvhjjj, we find the assertion.
h
Proof of Lemma 2.2. For arbitrary vh 2 V ph we set f := uh vh and g := u vh. Lemmas 2.1 and 2.2 imply 0 " #1=2 1 X X 1 2 2 Ajjjfjjj. jjjfjjj 6 As ðf; fÞ ¼ As ðg; fÞ þ dT ððL b LÞu; ~ b rfÞT 6 @Qs ðgÞ þ dT kðL b LÞukL2 ðT Þ 2 T 2Th T 2Th This estimate together with the triangle inequality conclude the proof.
h
Proof of Corollary 2.1. According to Lemma 2.2, the right-hand side of (21) is bounded by ( ) 2 X bkðL1 ðT ÞÞd 1 k~ 2 2 2 jjju U h jjj ¼ C jjju I h;p ujjj þ min ; ku I h;p ukL2 ðT Þ d a T T T 2Th þ
X
dT k r PT ðarðu I h;p uÞÞ þ cðu
2 I h;p uÞkL2 ðT Þ
T 2Th
þ
X
! 2 b dT kð L LÞukL2 ðT Þ .
ð51Þ
T 2Th
Now we estimate term by term using (18), and the interpolation estimate (22). The first one can be estimated by 2lT 2 X h2 hT kakL1 ðT Þ þ kckL1 ðT Þ T2 þ dT k~ bk2ðL1 ðT ÞÞd kuk2H kT ðT Þ . jjju I h;p ujjj2 6 3C 2I 2k T 2 p p T T T 2Th
ð52Þ
The second term in (51) is bounded in a similar way. Using the inverse inequality (12) and the definition of PT the third term of (51) can be bounded by ! T X X p4T hT2lT 2 h2l 2 2 2 2 2 T dT k r PT ðarðu I h;p uÞÞ þ cðu I h;p uÞkL2 ðT Þ 6 2C I dT linv kakL1 ðT Þ 2 2kT 2 þ kckL1 ðT Þ 2kT kuk2H kT ðT Þ . hT p T pT T 2Th T 2Th ð53Þ The last term of (51) can be estimated as follows: X X dT kð b L LÞuk2L2 ðT Þ ¼ dT kr ðaruÞ r PT ðaruÞk2L2 ðT Þ T 2Th
T 2Th
62
X
2 2 dT kr ðaruÞ r I Th;p ðaruÞkL2 ðT Þ þ kr I Th;p ðaruÞ r PT ðaruÞkL2 ðT Þ ;
T 2Th
I Th;p
is the local Lagrange interpolation. Using the approximation property (22) and the inverse inequality (12) yields where for the consistency error ! 2lT 4 4 X X 2 2 2 2 hT 2 pT T b dT kð L LÞukL2 ðT Þ 6 2 dT linv 2 kI h;p ðaruÞ PT ðaruÞkL2 ðT Þ þ C I 2kT 4 karukH kT 1 ðT Þ hT pT T 2Th T 2Th ! 2lT 4 T 2 X p4 h2l 2 hT T 62 dT l2inv C 2I T2 2k þ C ð54Þ kak2W kT 1;1 ðT Þ kuk2H kT ðT Þ ; I 2 2k T T 4 hT p T pT T 2Th where we have used the definition of PT and, again, the approximation property (22) in the last step. Summarizing (51)– (54), we obtain 2
jjju U h jjj 6 C
X h2ðlT 1Þ T T 2Th
2ðk T 1Þ
pT
2
kukH kT ðT Þ M T ðhT ; pT ; dT ; a; ~ b; cÞ
4138
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138
with (
2 bkðL1 ðT ÞÞd 1 k~ þ min ; dT aT
)
h2T þ dT k~ bk2ðL1 ðT ÞÞd p2T p2 þ kakL1 ðT Þ þ dT l2inv p2T kak2L1 ðT Þ þ ð1 þ l2inv p2T Þkak2W kT 1;1 ðT Þ T2 . hT
M T :¼ ðkckL1 ðT Þ þ
h2 dT kck2L1 ðT Þ Þ T2 pT
Applying the assumption (18) yields the assertion.
h
References [1] A.N. Brooks, T.J.R. Hughes, Streamline upwind Petrov–Galerkin formulation for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg. 32 (1982) 199–259. [2] E. Burman, A. Ern, Nonlinear diffusion and discrete maximum principle for stabilized Galerkin approximations of the convection–diffusion–reaction equation, Comput. Methods Appl. Mech. Engrg. 191 (2002) 3833–3855. [3] E. Burman, A. Ern, Stabilized Galerkin approximation of convection–diffusion–reaction equations: discrete maximum principle and convergence, Math. Comput., accepted for publication. [4] E. Burmann, P. Hansbo, Edge stabilization for Galerkin approximations of convection–diffusion problems, Comput. Methods Appl. Mech. Engrg. 193 (2004) 1437–1453. [5] 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. [6] R. Codina, O. Soto, Finite element implementation of two-equation and algebraic stress turbulence models for steady incompressible flows, Int. J. Numer. Meth. Fluids 90 (3) (1999) 309–343. [7] E.G. Do Carmo, G. Alvarez, A new stabilized finite element formulation for scalar convection–diffusion problems: the streamline and approximate upwind/Petrov–Galerlin method, Comput. Methods Appl. Mech. Engrg. 192 (2003) 3379–3396. [8] A. Ern, J.-L. Guermond, Theory and Practice of Finite Elements, Springer, 2004. [9] A.C. Galeao, E.G. Dutra do Carmo, A consistent approximate upwind Petrov–Galerkin method for convection-dominated problems, Comput. Methods Appl. Mech. Engrg. 68 (1988) 83–95. [10] A.C. Galeao, R.C. Almeida, S.M.C. Malta, A.F.D. Loula, Finite element analysis of convection dominated reaction–diffusion problems, Appl. Numer. Math. 48 (2004) 205–222. [11] B.-Y. Guo, H.-P. Ma, E. Tadmor, Spectral vanishing viscosity method for nonlinear conservation laws, SIAM J. Numer. Anal. 39 (4) (2001) 1254– 1268. [12] P. Houston, E. Su¨li, Stabilized hp-finite element approximations of partial differential equations with nonnegative characteristic form, Computing 66 (2001) 99–119. [13] T.J.R. Hughes, M. Mallet, A. Mizukami, A new finite element formulation for computational fluid dynamics: II. Beyond SUPG, Comput. Methods Appl. Mech. Engrg. 54 (1986) 341–355. [14] T.J.R. Hughes, G. Scovazzi, L.P. Franca, Multiscale and stabilized methods, in: E. Stein et al. (Eds.), Encyclopedia of Computational Mechanics, vol. 3, John Wiley and Sons, 2004. [15] C. Johnson, U. Na¨vert, Pitka¨ranta: finite element methods for linear hyperbolic problems, Comput. Methods Appl. Mech. Engrg. 45 (1984) 25–38. [16] C. Johnson, A.H. Schatz, Wahlbin: crosswind smear and pointwise errors in streamline diffusion finite element methods, Math. Comput. 49 (1987) 25–38. [17] T. Knopp, G. Lube, R. Gritzki, M. Ro¨sler, Iterative substructuring methods for incompressible non-isothermal flows and its application to indoor air flow simulation, Int. J. Numer. Meth. Fluids 40 (2002) 1527–1538. [18] T. Knopp, G. Lube, G. Rapin, Stabilized finite element methods with shock capturing for advection–diffusion problems, Comput. Methods Appl. Mech. Engrg. 191 (2002) 2997–3013. [19] G. Lube, G. Rapin, Residual-based stabilized higher-order FEM for a generalized Oseen problem, Math. Mod. Meth. Appl. Sci., accepted for publication. [20] U. Na¨vert, A Finite Element Method for Convection–Diffusion Problems, Ph.D. Thesis, Chalmers University of Technology, Go¨teborg, 1982. [21] C. Schwab, p- and hp-finite element methods. Theory and Applications to Solid and Fluid Mechanics, Oxford University Press, Oxford, 1998. [22] Y.-T. Shih, H.C. Elman, Iterative methods for stabilized discrete convection-diffusion problems, IMA J. Numer. Anal. 20 (2000) 333–358. [23] A. Szepessy, Convergence of a shock-capturing streamline diffusion finite element method to a scalar conservation law in two dimensions, Math. Comput. 53 (1989) 527–545. [24] A. Toselli, hp-finite element approximations on non-matching grids for partial differential equation with non-negative characteristic form, M2AN 37 (2003) 91–115. [25] G. Zhou, Local pointwise error estimates for the streamline diffusion method applied to nonstationary hyperbolic problems, East-West J. Numer. Math. 3 (1995) 217–235. [26] G. Zhou, R. Rannacher, Pointwise superconvergence of the streamline diffusion finite element method, Numer. Methods Part. Diff. Equat. 12 (1996) 123–145.
Residual-based stabilized higher-order FEM for advection-dominated problems Gert Lube *, Gerd Rapin Mathematics Department, NAM, University of Go¨ttingen, Gottingen D-37083, Germany Received 4 January 2005; received in revised form 7 July 2005; accepted 26 July 2005
Abstract We reconsider the numerical solution of linear(ized) advection–diffusion–reaction problems using higher-order finite elements together with stabilized Galerkin methods of streamline-diffusion type (SUPG) and with shock-capturing stabilization. The analysis improves the a priori analysis in our previous paper [T. Knopp, G. Lube, G. Rapin, Stabilized finite element methods with shock capturing for advection–diffusion problems, Comput. Methods Appl. Mech. Engrg. 191 (2002) 2997–3013]. The theoretical results are supported by some numerical experiments. Ó 2005 Elsevier B.V. All rights reserved. Keywords: Advection–diffusion; Incompressible flow; Stabilized finite elements; SUPG
1. Introduction The motivation of the present paper stems, e.g., from the finite element simulation of the non-isothermal and incompressible Navier–Stokes problem ot~ u r ðmr~ uÞ þ ð~ u rÞ~ u þ rp ¼ bh~ g; r ~ u ¼ 0; ot h þ ð~ u rÞh r ðarhÞ ¼ q_ V =cp
ð1Þ ð2Þ ð3Þ
for velocity ~ u, pressure p and temperature h in a polyhedral domain X Rd, d 6 3, with source terms bh~ g and q_ V =cp . This model describes, e.g., the air flow in buildings, etc. [17]. The momentum and continuity Eqs. (1) and (2) describe the fluid motion; the heat transfer is driven by the advection–diffusion Eq. (3). Turbulence may occur at high Rayleigh or Reynolds numbers. A standard approach is to consider the Reynolds averaged Navier–Stokes equations (RANS) together with, e.g., the k– turbulence model. Within a statistical turbulence model only averaged values are considered. An eddy viscosity ansatz for turbulent effects is modeled as an additional diffusion term with eddy viscosity mt. Then the averaged values for ~ u, p and h are determined by (1)–(3) with m and a replaced with (variable) viscosities me = m + mt and ae = a + mt/Prt. The eddy viscosity term mt is determined, e.g. in the k– turbulence model, by mt = clk2/ where the turbulent kinetic energy k and the dissipation rate of k are defined by additional (nonlinear) advection–diffusion–reaction equations. A standard algorithmic treatment of the coupled model is to semi-discretize, in an outer loop, in time (with possible step control) using an A-stable method and then, in an inner loop, to decouple and linearize the resulting system. A block *
Corresponding author. E-mail addresses: [email protected] (G. Lube), [email protected] (G. Rapin).
0045-7825/$ - see front matter Ó 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cma.2005.07.017
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138
4125
Gauss–Seidel method with fixed point or Newton-type iteration per time step leads to linearized Navier–Stokes problems (of Oseen-type) and linearized advection–diffusion–reaction problems as auxiliary problems, see [17]. A proper numerical approach to the latter model Lu :¼ r ðaruÞ þ ~ b ru þ cu ¼ f
ð4Þ
in X
is an important ingredient of the approach and will be discussed in this paper. The streamline upwind Petrov–Galerkin (SUPG) method, proposed by Brooks and Hughes [1], was the first variationally consistent, stable and accurate finite element model for advection-dominated problems. It initiated the development of stabilization techniques for advection-dominated and related problems. For an overview, see, e.g., [14]. A first relevant analysis of the SUPG method can be found in Johnson et al. [15] in case of regular solutions. The theory has been refined over the years in several directions. Here we only mention the analysis of a hp-version by Houston and Su¨li [12]. Nevertheless, for non-smooth solutions, localized oscillations of the SUPG solution may still exist in the neighborhood of steep gradients. As a remedy, discontinuity- or shock-capturing terms can be added to enhance the stability. Linear (but non-consistent) schemes for low-order elements are considered, e.g., in [16,22]. Mizukami/Hughes [13] introduced the first nonlinear discontinuity-capturing schemes DC1 and DC2. The idea was to enhance, additionally to streamline upwinding, numerical viscosity in the direction of $uh. The consistent approximate upwind (CAU) in [9] provided a further refinement. For recent developments of the CAU scheme to higher-order elements, we refer to [10]. Moreover, Codina considered in [5,6] the discontinuity-capturing/crosswind-dissipation (DC/CD) scheme with additional anisotropic viscosity. The important question of low-order (nonlinear) schemes which satisfy a discrete maximum principle is discussed in the recent papers [2,3], see also the monograph [8]. A first theoretical result for such nonlinear schemes is seemingly due to Szepessy [23]. In our previous paper [18] we considered the a priori analysis of a rather general class of shock-capturing schemes. The goal of the present paper is an extension of the stabilized higher-order FE method of the recent paper [12] to the case of shock-capturing stabilization. In particular, we address the choice of the stabilization parameters, see Section 2, extend and refine the analysis of shockcapturing schemes given in [18], see Section 3, and provide some numerical experiments, see Section 4. 2. Stabilized FEM for advection–diffusion–reaction model Following basically [12], we describe and analyze the SUPG-stabilization of the advection–diffusion–reaction model. In contrast to [12], we give a refined definition of the stabilization parameters depending on all critical parameters. 2.1. Problem statement For the advection–diffusion–reaction scheme (4), we assume a; c 2 L1 ðXÞ, ~ b 2 ðH 1 ðXÞÞ \ ðL1 ðXÞÞ , f 2 L2(X) and d
ðr ~ bÞðxÞ ¼ 0;
cðxÞ P x P 0;
aðxÞ P a0 > 0;
a.e. in X.
d
ð5Þ
For simplicity only, we analyze the homogeneous Dirichlet problem u¼0
ð6Þ
on oX.
The basic variational formulation of (4)–(6) reads: Find u 2 V :¼ H 10 ðXÞ
s.t. Aðu; vÞ ¼ lðvÞ
8v 2 V
ð7Þ
with Aðu; vÞ ¼ ðaru; rvÞX þ ð~ b ru þ cu; vÞX ; lðvÞ ¼ ðf ; vÞX .
ð8Þ ð9Þ
2.2. Finite element discretization Suppose a family of admissible triangulations Th ¼ fT g of the polyhedral domain X where h is the piecewise constant mesh function with h(x) = hT = diam (T), x 2 T. We assume that Th is shape-regular, i.e. there exists a constant Cr 5 Cr(h) such that C r hdT 6 meas ðT Þ
8T 2 [h Th .
ð10Þ
Moreover, we assume that each element T 2 Th is a smooth bijective image of a given reference element Tb , i.e., T ¼ F T ð Tb Þ for all T 2 Th . Here, Tb is the (open) unit simplex or the (open) unit hypercube in Rd.
4126
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138
For p 2 N, we denote by Pp ð Tb Þ the set of polynomials of degree at most p on Tb . To each T 2 Th an integer pT P 1 is assigned; then we define the vector p ¼ fpT : T 2 Th g. We introduce a conforming finite element space on Th by V ph ¼ X ph \ H 10 ðXÞ;
b X ph ¼ fv 2 CðXÞjvj T F T 2 P pT ð T Þ 8T 2 Th g.
ð11Þ
A key point in the subsequent analysis are the following local inverse inequalities 8w 2 X ph
krwkðL2 ðT ÞÞd 6 linv p2T h1 T kwkL2 ðT Þ
ð12Þ
on T 2 Th . The constant linv depends only on the shape-regularity parameter Cr in (10), see [21]. Interpolation estimates will be given later on. The standard Galerkin-FEM is, in general, not appropriate for the numerical solution of problem (7)–(9) in the singularly perturbed case. Therefore we consider the well-known streamline upwinding (SUPG)-stabilization: Find uh 2 V ph : As ðuh ; vÞ ¼ Ls ðvÞ with As ðu; vÞ :¼ Aðu; vÞ þ Ls ðvÞ :¼ lðvÞ þ
X
X
8v 2 V ph
ð13Þ
dT ð b Lu; ~ b rvÞT ;
ð14Þ
T 2Th
dT ðf ; ~ b rvÞT ;
ð15Þ
T 2Th
b LujT :¼ r PT ðaruÞ þ ~ b ru þ cu 2
ð16Þ
d
d
and the orthogonal projection PT : ½L ðT Þ ! ½P pT ðT Þ . This modification of the (variable) diffusion term allows the application of inverse inequalities [12]. Remark. Other residual-based stabilization variants are the Galerkin/least-squares (GLS) method and the algebraic subgrid scale (ASGS) method where the test function dT ð~ b rvÞ is replaced with dTLv and dT(L v), respectively. All these P 1=2 dT k~ b rðÞk2 . Our numerical methods combine stability and accuracy and allow additional error control of T 2Th
0;T
experience shows that the additional test terms of the GLS and ASGS methods are less important. For this reason, we restrict the subsequent analysis to the SUPG-case. 2.3. Basic quasi-optimal error estimate The following estimates of the bilinear form As, using the stabilized norm ! 1=2 X pffiffiffi p ffiffi ffi 2 2 2 jjjvjjj :¼ k arvk 2 d þ k cvk 2 þ dT k~ b rvk 2 ; L ðT Þ
ðL ðT ÞÞ
ð17Þ
L ðT Þ
T 2Th
are slight modifications of results given in [12,18] and will be important in Section 3, too. For the convenience of the reader, we include the proofs in the Appendix. First of all, we obtain the following stability and continuity results. Lemma 2.1. Let the assumptions (5) and the following condition be valid: ! 1 h2T 1 0 6 dT 6 min 4 2 ; . 2 pT linv kakL1 ðT Þ kckL1 ðT Þ
ð18Þ
Then the bilinear form As in (14) satisfies the stability estimate 1 2 As ðv; vÞ P jjjvjjj 2
8v 2 V h .
ð19Þ
Moreover, for u 2 H 10 ðXÞ with ($ Æ (a$u))kT 2 L2(T) and vh 2 V ph , we obtain As ðu; vh Þ 6 Qs ðuÞjjjvh jjj Qs ðuÞ :¼ jjjujjj þ
X T 2Th
with aT := infx2Ta(x).
(
) !1=2 2 bkðL1 ðT ÞÞd 1 k~ 2 min ; þ kukL2 ðT Þ dT aT
X T 2Th
!1=2 dT k r PT ðaruÞ þ
2 cukL2 ðT Þ
ð20Þ
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138
4127
We combine both results to prove an a priori estimate for the SUPG scheme. Lemma 2.2. Let u 2 V be a solution of (7)–(9) with ($ Æ (a$u))kT 2 L2(T) for all T 2 Th and let the assumptions (5) and (18) be valid. This implies the quasi-optimal a priori error estimate of the SUPG scheme (13)–(16): !1=2 X 2 b LÞðuÞk 2 jjju uh jjj 6 C 1 inf Qs ðu vh Þ þ C 2 dT kð L ð21Þ L ðT Þ
p
vh 2V h
T 2Th
with data-independent constants Ci. In order to determine the stabilization parameter dT we prove an error estimate using appropriate approximation properties. Let be v 2 Hk(T) with k P r > d2. The standard Lagrange interpolation I Th;p v satisfies kv I Th;p vkH s ðT Þ 6 C I
hls T kvkH k ðT Þ ; pks T
0 6 s 6 l ¼ minðpT þ 1; kÞ.
ð22Þ
But we need a global interpolation operator. Unfortunately we cannot use the global Lagrange interpolation Ih,p for nonconstant pT, since then the Lagrange interpolation is in general not continuous and therefore it is not included in V ph . Nevertheless, in order to avoid additional technical problems, we require that I h;p v 2 V ph is satisfied, which is valid if p = pT for all T 2 Th . Then the global Lagrange interpolation Ih,pv for all v 2 H 10 ðXÞ \ H r ðXÞ; r > d2 with vjT 2 Hk(T), k P r, satisfies (22) with s 2 {0, 1} for all T 2 Th . Remark (i) The assumption p = pT on Th is quite restrictive. A remedy is a non-conforming domain decomposition method where a constant polynomial degree is assumed in each subdomain. The polynomial degree can differ from subdomain to subdomain, cf. [24]. (ii) The assumption r > d2 can be relaxed to r > 12 if the Scott–Zhang quasi-interpolant is used. But the right-hand side norm in (22) has to be replaced by the Hk-norm over a patch of elements containing T . The following result is a slight modification of Theorem 9 in [12]. Corollary 2.1. Let the assumptions (5), (18) be satisfied. Assume that the weak solution u 2 H 10 ðXÞ of (7)–(9) belongs locally to H kT ðT Þ, k T > d2 and a 2 W kT 1;1 ðT Þ for all T 2 Th . Moreover, we assume that I h;p u 2 V ph ðXÞ. Then the following error estimate for the solution uh of the SUPG scheme (13)–(16) is valid: X h2ðlT 1Þ T
2
jjju uh jjj 6 C
T 2Th
þ
2
2
kukH kT ðT Þ kakL1 ðT Þ þ 2ðk T 1Þ
kakW kT 1;1 ðT Þ
pT
kckL1 ðT Þ h2T p2T
kakL1 ðT Þ
(
þ
dT k~ bk2ðL1 ðT ÞÞd
2 bkðL1 ðT ÞÞd 1 k~ þ min ; dT aT
)
! h2T . p2T
ð23Þ
2.4. Parameter design Up to now, we have not specified the stabilization parameter dT. Here we revisit the analysis already considered in [12]. Corollary 2.2. Let the assumptions of Corollary 2.1 be valid. The choice of the stabilization parameters according to ( ) hT 1 h2T dT min ; ; ð24Þ kckL1 ðT Þ p4T l2inv kakL1 ðT Þ p k~ bk 1 T
L ðT Þ
is consistent with condition (18). Moreover, it leads to X h2ðlT 1Þ opt 2 2 T M T kukH kT ðT Þ jjju uh jjj 6 C 2ðk T 1Þ T 2Th pT
ð25Þ
with
2
M opt T
:¼ kakL1 ðT Þ 1 þ
kakW kT 1;1 ðT Þ 2
kakL1 ðT Þ
þ PeT þ CT þ min
kakL1 ðT Þ aT
Pe2T ; p2T l2inv
! ð26Þ
4128
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138
and PeT :¼
hT k~ bkðL1 ðT ÞÞd pT kakL1 ðT Þ
CT :¼
;
kckL1 ðT Þ h2T p2T kakL1 ðT Þ
ð27Þ
.
Proof. The idea is to equilibrate the dT-dependent right-hand side terms of the error estimate (23) under the side conditions (18) and continuity of dT. 1 dT
We consider the case are equilibrated if dT
k~ bk2
6
ðL1 ðT ÞÞd
aT
in the min-term of the right-hand side of (23). Then the dT -dependent terms in (23)
hT . ~ pT kbkðL1 ðT ÞÞd
ð28Þ
In order to fulfil the condition (18), we set dT according to (24). This leads to ( ) ( ( ) 2 ) 2 2 bkðL1 ðT ÞÞd h2T hT k~ bkðL1 ðT ÞÞd hT k~ bkðL1 ðT ÞÞd kckL1 ðT Þ h2T 2 2 1 k~ Z T :¼ min ; ¼ min max ; ; pT linv kakL1 ðT Þ ; dT aT p2T pT p2T p2T aT kakL1 ðT Þ 2 2 2 ¼ kakL1 ðT Þ min PeT ; max PeT ; CT ; pT linv ; aT hence M opt T :¼ kakL1 ðT Þ þ 6 kakL1 ðT Þ 6 kakL1 ðT Þ
kak2W kT 1;1 ðT Þ kakL1 ðT Þ
þ
kckL1 ðT Þ h2T p2T
2
þ dT k~ bkðL1 ðT ÞÞd þ Z T
! ZT 1þ þ PeT þ CT þ kakL1 ðT Þ kak2L1 ðT Þ ! 2 kakW kT 1;1 ðT Þ kakL1 ðT Þ 2 2 2 1þ þ PeT þ CT þ min PeT ; pT linv . 2 aT kakL1 ðT Þ kak2W kT 1;1 ðT Þ
This leads to the estimate (25)–(27).
h
Remark. Let us discuss the optimality of the result. For brevity, we consider the case a(x) a. The estimate (25)–(27) gives in the advection/reaction-dominated limit maxfPeT ; CT g P p2T l2inv the estimates X
hT
T 2Th
pT k~ bkðL1 ðT ÞÞd
X T 2Th
kð~ b rÞðu uh Þk2L2 ðT Þ 6 C
X h2lT 1 kuk2H kT ðT Þ ; 2k T 1 p T 2Th
X h2lT pffiffiffi 1 2 2 k cðu uh ÞkL2 ðT Þ 6 C kukH kT ðT Þ ; 2k T kckL1 ðT Þ p T 2Th
PeT P CT ;
CT P PeT .
These estimates are simultaneously optimal in h and p with respect to the (weighted) streamline-derivative in the advectiondominated limit PeT P CT P p2T l2inv , and in the (weighted) L2-norm in the reaction-dominated limit CT P PeT P p2T l2inv , respectively, see also [12], Remark 10a. Please note that both estimates are also valid for the complete error term jjj Æ jjj2. In the diffusion-dominated limit with PeT 6 1 ð6 p2T l2inv Þ, we obtain for the diffusion-relevant part of the jjj Æ jjj-norm 2
akrðu uh ÞkL2 ðXÞ 6 C
X T 2Th
a
h2ðlT 1Þ 2 kukH kT ðT Þ ; p2ðkT 1Þ
PeT 6 1
which is simultaneously optimal in h and p with respect to the H1-seminorm, see also [12], Remark 10a. Nevertheless, there remains a transition region between the diffusion-dominated limit and the advection/reaction-dominated limit with loss of optimality with respect to p by p1=2 (at least using the present technique). We want to point out that a similar analysis has been given in Remark 10 of [12], but the analysis is seemingly not complete in the transition region. Remark. We would like to point out that the dependence of the parameters dT on the spectral order p is specified. For fixed hT, it tends to zero with increasing p. This is in (formal) agreement with the behaviour of the artificial viscosity in the spectral vanishing viscosity method, see, e.g., [11].
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138
4129
3. Shock-capturing stabilization The SUPG method stabilizes most of the unphysical oscillations of the Galerkin discretization which are caused by ? dominating advection. Nevertheless, there may remain spurious localized oscillations, often in crosswind-directions ~ b . Such perturbations are very dangerous in the numerical solution of coupled nonlinear problems, see Section 1. A wellknown remedy is the shock-capturing variant of the SUPG method. For a review of such methods we refer to [22]. Versions which add linear crosswind-stabilization terms can be constructed for piecewise linear and bilinear elements. Energy estimates and even maximum norm estimates can be found e.g. in [16,22,26]. In practice, these methods suffer from smearing in crosswind-directions since, due to the linearity of the schemes, the stabilization effect cannot be restricted to certain subregions (e.g. close to shocks). An interesting new (but nonlinear) scheme in [2,3] for piecewise linear elements allows to satisfy a discrete maximum principle. 3.1. Framework of shock-capturing The extension of shock-capturing schemes to higher order elements requires a nonlinear construction. Different proposals in the literature [5,9,13,22] are covered by the following class of schemes [18]: Find U h 2 V ph : As ðU h ; vÞ þ Asc ðU h ; U h ; vÞ ¼ Ls ðvÞ with Asc ðw; u; vÞ ¼
X
8v 2 V ph
ðsT ðwÞDsc ru; rvÞT ;
ð29Þ
ð30Þ
T 2Th
where Dsc denotes a symmetric, positive semi-definite matrix function with kDsc kL1ðXÞd d 6 1. The non-negative limiter function sT(w) has to restrict the effect of shock-capturing to subregions where the (modified) residual b Lw f is too large. The common feature of different proposals is that sT is a function of the scaled (modified) residual, namely sT ðwÞ ¼ sT ðRT ðwÞÞ;
RT ðwÞ :¼
kb Lw f kL2 ðT Þ kwkH 1 ðT Þ þ jT
ð31Þ
;
where we introduced a regularization parameter jT > 0. Standard shock-capturing schemes [5,9,13,22] with jT = 0 lead to ill-posed nonlinear problems. For a possible choice of jT, see Example 3.1. Some results for the h-version of this scheme can be found in [18]. Here we refine and extend the analysis to the hp-version of the shock-capturing scheme (29) and (30). Lemma 2.1 allows to prove the solvability of the nonlinear scheme (29)– (31). Unfortunately, a (local) uniqueness result is still open. Theorem 3.1. Let the assumptions (5) and (18) be valid. Moreover, assume that the limiter function sT(w) in (30) is continuous w.r.t. w. The SUPG-scheme with shock-capturing stabilization (29) and (30) has at least one solution U h 2 V ph . The following a priori estimate is valid for each solution 2 X pffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1=2 2 ð32Þ jjjU h jjj2 þ sT ðU h ÞDsc rU h 2 6 Cjjjf jjj L ðT Þ
T 2Th
with the dual norm jjjf jjj :¼ supv2V p ;jjjvjjj¼1 ðf ; vÞX . h
Proof. The proof relies on a variant of Brouwers fixed point theorem where Lemma 2.1 is used; for details see [18], Theorem 3.5. h 3.2. Shock-capturing with isotropic diffusion A priori error estimates of scheme (29) and (30) require more information about the shock-capturing term Asc. Let us first consider the case of isotropic diffusion with almost quadratical dependence of the function sT(w) on the scaled residual RT ðwÞ according to 2
Diso sc ¼ I;
2 iso siso T ðwÞ :¼ lT ðwÞ½RT ðwÞ ¼
liso T ðwÞ½RT ðwÞ
½kwkH 1 ðT Þ þ jT 2
ð33Þ
Lwh f kL2 ðT Þ . The class of schemes (29) and (30) with (33) covers, in particular, the consistent approximate with RT ðwh Þ ¼ k b upwind (CAU) scheme of [9] with jT = 0; it is also closely related to the original discontinuity-capturing schemes DC1 and
4130
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138
DC2 of [13]. The basic idea of these schemes was to introduce, additionally to streamline upwinding, numerical viscosity in the direction of $Uh. The following global a priori estimate shows that the convergence properties of the method are at least as good as for the SUPG-FEM. Moreover, we have stronger control of the gradient, including the crosswind-direction(s). Let again vh ¼ I h;p u 2 V ph be an appropriate interpolant of the solution u 2 V of (7)–(9). In the forthcoming analysis we introduce the quantities iso iso F iso h ¼ F h ðI h;p u; U h Þ :¼ sup F T ðI h;p u; U h Þ; T
F iso T ðvh ; U h Þ :¼
jvh jH 1 ðT Þ jU h jH 1 ðT Þ þ jT
.
ð34Þ
Assume, for simplicity, that Ih,pu can be taken as the Lagrange interpolant of the (sufficiently smooth) solution u. Then it is shown in [18], Remark 4.6 that the constant F iso h is uniformly bounded w.r.t. the mesh size h and viscosity a. See also next remark. p Theorem 3.2. Let the assumptions (5), (18), and (33) with 0 6 liso T ðwÞ 6 qdT for all w 2 V h and a sufficiently small q > 0 be 2 satisfied. Assume that the weak solution u 2 V of (7)–(9) satisfies locally ($ Æ (a$u))jT 2 L (T) and belongs to H kT ðT Þ; k T > d2 for all T 2 Th . Moreover, we assume that vh :¼ I h;p u 2 V ph . Then the following error estimate for the solution Uh of the stabilized discrete problem (29), (30) and (33) is valid: X X
2 2 2 2 siso þ K2 dT kð b L LÞðuÞkL2 ðT Þ ð35Þ jjju U h jjj þ T ðU h Þju U h jH 1 ðT Þ 6 K 1 Qs ðu I h;p uÞ T 2Th
T 2Th
with data-independent constants Ki. Proof. We set e := Uh vh. and g := u vh. Lemma 2.1 and the definition of Aiso sc imply 1 iso iso iso jjjejjj2 þ Aiso sc ðU h ; e; eÞ 6 As ðe; eÞ þ Asc ðU h ; e; eÞ ¼ As ðU h ; eÞ þ Asc ðU h ; U h ; eÞ As ðvh ; eÞ Asc ðU h ; vh ; eÞ 2 X ¼ As ðg; eÞ Aiso dT ðð b L LÞu; ~ b reÞT I þ II þ III. sc ðU h ; vh ; eÞ
ð36Þ
T 2Th
We estimate the right-hand side terms separately. Lemma 2.1 implies 1 2 2 I ¼ As ðg; eÞ 6 jjjejjj þ 2½Qs ðgÞ . 8
ð37Þ
The estimate of term II takes advantage of the structure (33) of siso T ðwÞ and of the terms defined in (34): X
II 6
T
6
!1=2 2 siso T ðU h ÞjejH 1 ðT Þ
X
!1=2 2 siso T ðU h Þjvh jH 1 ðT Þ
T
3 X iso 1 X iso 2 2 2 sT ðU h ÞjejH 1 ðT Þ þ s ðU h Þ½F iso T ðjU h jH 1 ðT Þ þ jT Þ 4 T 3 T T
3 1 X iso 2 2 6 Aiso l ðU h Þ½F iso sc ðU h ; e; eÞ þ T RT ðU h Þ 4 3 T T X 3 1 iso 2 6 Aiso dT R2T ðU h Þ sc ðU h ; e; eÞ þ q½F h 4 3 T X 3 iso 2 2 2 2 b b b 6 Asc ðU h ; e; eÞ þ q½F iso d k Lek þ k Lgk þ kð L LÞuk 2 2 2 T L ðT Þ L ðT Þ L ðT Þ . h 4 T 2Th
ð38Þ
In the last substep we used the Young inequality together with LU h f kL2 ðT Þ 6 k b LðU h vh ÞkL2 ðT Þ þ k b Lðvh uÞkL2 ðT Þ þ kð b L LÞukL2 ðT Þ . RT ðU h Þ ¼ k b For the last right-hand side term in (36) we find III 6
X 1 X 2 2 dT k~ b rekL2 ðT Þ þ 2 dT kð b L LÞukL2 ðT Þ . 8 T 2Th T 2Th
ð39Þ
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138
4131
Summarizing (36)–(39) we obtain X X 1 1 2 2 2 2 2 iso 2 iso 2 b b ½ jjjejjj þ Aiso ðU ; e; eÞ 6 2 Q ðgÞ þ q½F d k Lek þ k Lgk dT kð b L LÞukL2 ðT Þ . 2 2 h T s L ðT Þ L ðT Þ þ ðq½F h þ 2Þ h sc 4 4 T 2Th T 2Th ð40Þ Inverse inequalities (12) and assumption (18) imply X pffiffiffi 2 1 X pffiffiffi 2 2 dT k r PT ðareÞ þ cekL2 ðT Þ 6 k arvkðL2 ðT ÞÞd þ k cvkL2 ðT Þ . 2 T 2Th T 2Th 2
1 Using the triangle inequality together with the latter estimate, we find with q½F iso h 6 16 that X X 2 2 2 2 2 q½F iso dT k b LekL2 ðT Þ 6 2q½F iso dT k~ b rekL2 ðT Þ þ kr PT ðareÞ cekL2 ðT Þ h h T 2Th
T 2Th
1 2 2 2 6 2q½F iso h jjjejjj 6 jjjejjj . 8
ð41Þ
Please note that we take advantage of the assumption that vh = Ih,pu is the Lagrange interpolant of u 2 V in order to have a proper bound of F iso h . Furthermore, we obtain via triangle inequality and the definition of Qs(g) that 1 X X 2 2 2 2 2 iso 2 ~ b q½F iso d k Lgk 6 2q½F d k b rgk þ kr P ðargÞ cgk ð42Þ 2 2 2 T T T L ðT Þ L ðT Þ L ðT Þ 6 ½Qs ðgÞ . h h 8 T 2Th T 2Th We collect (40)–(42) to X 1 2 2 2 ~ ~ jjjejjj þ Aiso dT kð b L LÞukL2 ðT Þ sc ðU h ; e; eÞ 6 K 1 ½Qs ðu vh Þ þ K 2 2 T 2Th ~ i . The triangle inequality concludes the proof, using that the term Aiso ðU h ; g; gÞ can be estimated with appropriate positive K SC similarly as the term II. h Remark. In order to relax the regularity assumptions on u 2 V in Theorem 3.2, let vh ¼ I h;p u 2 V ph be chosen as a quasiinterpolant of u 2 V with the local stability estimate jvh jH 1 ðT Þ 6 CjujH 1 ðV ðT ÞÞ on a suitable neighborhood of V(T) of element T. Then F iso h is uniformly bounded w.r.t. h, see [18], Remark 4.6. The present analysis yields so far only an upper bound for the limiter function liso T ðU h Þ 6 qdT ; it gives no hint to its explicit construction. Example 3.1. For the special case of a = const. and of piecewise linear elements (i.e. p = 1), the CAU-scheme in [9] is determined by ! ~
hT j~ bk j d0 hT cau lT ðU h Þ ¼ max 0; dk dT ; dk ¼ min 1; ; 6a j~ bk j ~ brU h rU h if j$Uhj 5 0 and ~ bk ¼ 0 else. Numerical results, see, e.g. [7,9,13], show that this scheme is suitable for where ~ bk ¼ jrU j2 h
solutions with strong gradients in layers. For problems with smooth solutions, however, it suffers from smearing in crosswind-direction(s). An improved version of the CAU scheme, even for higher-order elements, has been discussed in the recent paper [10]. The analysis (with some open ends) is based on an linearized iterative scheme. The present analysis allows to set lcau T 6 qdT . Other analytical tools are necessary to close the gap between theoretical and practical results. 3.3. Shock-capturing with crosswind-diffusion The crosswind-smearing of the DC- and CAU-schemes lead Codina in [5,6] to the construction of a shock-capturing scheme with anisotropic diffusion in crosswind-direction(s). This discontinuity-capturing crosswind-dissipation (DC/CD) scheme can be seen as a special case of the general schemes (29) and (30) with crosswind-diffusion and almost linear dependence on RT ðwÞ: 8 ~ ~ > jwj : H 1 ðT Þ þ jT ~ 0; b¼0
4132
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138
Codinas scheme is introduced in [6] with jT = 0 and the limiter function 2kakL1 ðT Þ dc=cd cd lT ðwÞ ¼ lT ðwÞ :¼ l0 hT max 0; b hT RT ðwÞ
ð44Þ
h R ðwÞ
T T can be seen as a pseudo mesh Peclet number. Numerical experiments show that b is a critical with l0 ¼ 12. The term 2kak L1 ðT Þ parameter. Codina proposed b = 0.7 for p = 1 and b = 0.35 for p = 2. Here we try to mimic the proof of Theorem 3.2. The following result shows that this is possible under the following restriction on the limiter function
0 6 lcd T ðwÞ 6 qdT RT ðwÞ;
w 2 V ph
ð45Þ
with a suitable parameter q > 0. This means that both variants of the shock-capturing scheme with isotropic and with anisotropic diffusion have the same upper bounds of the additional shock-capturing term. Theorem 3.3. Let the assumptions of Theorem 3.2 be valid. Then we obtain for the shock-capturing scheme (29), (30) and (43) with assumption (45) the estimate: X X
2 2 2 2 cd 1=2 scd rðu U h ÞkL2 ðT Þ 6 K 3 Qs ðu I h;p uÞ þ K 4 dT kð b L LÞukL2 ðT Þ ð46Þ jjju U h jjj þ T ðU h Þk½Dsc T 2Th
T 2Th
with data-independent constants Ki. Proof. We modify the proof of Theorem 3.2 and introduce the quantities cd F cd h ðI h;p u; U h Þ :¼ sup F T ðI h;p u; U h Þ; T
F cd T ðI h;p u; U h Þ :¼
1=2 k½Dcd rI h;p ukL2 ðT Þ sc
jU h jH 1 ðT Þ þ jT
ð47Þ
.
The constant F cd h is uniformly bounded w.r.t. h and a. We start again from (36), but the estimate of term II is different depending on the structure of scd T ðwÞ. Proceeding similarly as in (38) and using (47), (43) and (45) we obtain
II 6
X
!1=2 2 cd 1=2 scd rI h;p ukL2 ðT Þ T ðU h Þk½Dsc
X
T
6
!1=2 2 cd 1=2 scd rekL2 ðT Þ T ðU h Þk½Dsc
T
3 X cd 1 X cd 1=2 2 2 2 sT ðU h Þk½Dcd rekL2 ðT Þ þ s ðU h Þ½F cd sc T ðjU h jH 1 ðT Þ þ jT Þ 4 T 3 T T
3 1 X cd 2 6 Acd l ðU h Þ½F cd sc ðU h ; e; eÞ þ T RT ðU h ÞðjU h jH 1 ðT Þ þ jT Þ 4 3 T T X 3 1 cd 2 6 Asc dT R2T ðU h Þ. sc ðU h ; e; eÞ þ q½F h 4 3 T Assuming a sufficiently small q, we proceed as in the proof of Theorem 3.2.
ð48Þ h
As in the previous case (cf. discussion of Theorems 3.2), 3.3 provides only an upper bound of the limiter function lcd T ðwÞ 6 qdT RT ðwÞ. As already discussed for the CAU-scheme, other analytical tools are necessary to refine the design of the limiter functions lT(Æ). p Remark. We propose to choose the parameter jT in such a way that lcd T ðwÞ remains bounded. For all w 2 V h , the inverse estimate (12) leads to
RT ðwÞ ¼ 6
k r PT ðarwÞ þ ~ b rw þ cw f kL2 ðT Þ jwjH 1 ðT Þ þ jT linv p2T hT
1 kakL1 ðT Þ þ k~ bkðL1 ðT ÞÞd þ kckL1 ðT Þ kwkL2 ðT Þ þ kf kL2 ðT Þ . jT
Choosing jT ¼ kckL1 ðT Þ kwkL2 ðT Þ þ kf kL2 ðT Þ
ð49Þ
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138
we obtain
(
dT RT ðwÞ
1 h2 ; ; 4 2 T kckL1 ðT Þ pT linv kakL1 ðT Þ L1 ðT Þ
hT 6 min pT k~ bk
)
4133
linv p2T kakL1 ðT Þ þ k~ bkðL1 ðT ÞÞd þ 1 . hT
This means that the artificial viscosity in Acd sc is at most of order hT/pT. 4. Numerical experiments Let us verify the error estimates of Sections 2 and 3, starting with a problem with a smooth solution. Example 4.1. Consider problem (4)–(6) with ~ b ¼ ðx2 ; x1 ÞT , a 2 {103,106}, and c 2 {0,103} in the domain X = (0, 1)2. The source term f is chosen in such a way that uðxÞ ¼ sinðpx1 Þ sinðpx2 Þex1 x2 is the (smooth) exact solution of (4)–(6). The numerical solution of the SUPG-scheme (13)–(16) with the parameter choice according to (24) is calculated with FEMLAB 2.3 on quasi-uniform meshes for Pp -elements with p 2 {1, 2, . . ., 6}. The additional shock-capturing stabilization is not required here. pffiffiffi 2 In Fig. 1 we present the p-convergence of the numerical error measured in the energy norm jjjvjjjE :¼ ðk arvkL2 ðXÞ þ 2 1=2 b cover all cases kvkL2 ðXÞ Þ for different values of h. The chosen values of c, a and h, p together with the variable flow field ~ of the parameter design for dT according to (24). The convergence rates, including the h-convergence for fixed p (not shown here), confirm the theoretical result of Corollary 2.2. In general, we observed for problems with a smooth solution a remarkable robustness of the results with respect to the design of the stabilization parameters. This is a consequence of the residual-based stabilization terms. We consider now the shock-capturing method (29)–(31) in case of rough solutions where typically unphysical oscillations appear around layers of the continuous solution. Interior error estimates, i.e. in subdomains X 0 X away from 0
0
10
10 h=1/4 h=1/8 h=1/16 h=1/32
–2
10
10
–4
–4
10
energy
energy
10
–6
10
–8
–6
10
–8
10
10
–10
–10
10
10
–12
10
h=1/4 h=1/8 h=1/16 h=1/32
–2
1
–12
2
3
p
4
5
10
6
1
2
3
0
5
6
0
10
10 h=1/4 h=1/8 h=1/16 h=1/32
–2
10
h=1/4 h=1/8 h=1/16 h=1/32
–2
10
–4
–4
10
10
–6
energy
energy
4
—
—
10
–8
10
–10
–8
10
10
–12
–12
10
10
–14
1
–6
10
–10
10
10
p
–14
2
3
p
—
4
5
6
10
1
2
3
p
4
5
6
—
Fig. 1. p-convergence for Example 4.1 with smooth solution and different values of a 2 {103, 106}, c 2 {0, 103}, and h 2 f14 ; 18 ; 161 ; 321 g.
4134
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138
boundary and interior layers, were given, e.g., by [20,25]. It is well-known that the nonlinear discrete problem (29)–(31) is very ill-posed for jT = 0, cf. [2,22]. Here we solve the nonlinear problem via simple iteration n nþ1 s n 2 N0 : As ðU nþ1 h ; vÞ þ Asc ðU h ; U h ; vÞ ¼ L ðvÞ
8v 2 V ph .
ð50Þ
Using the initial guess U 0h ¼ 0, the first iterate is given by the solution of the SUPG scheme (13)–(16). For a related approach for the CAU scheme, we refer to [10]. The iteration (50) with the rather strong stopping criterion kU hnþ1 U nh k1 6 104 is robust and fast for p = 1. For p P 2 we sometimes observed that the simple iteration did not converge. Then the application of Newton-type methods is recommended, see, e.g., [22]. Example 4.2. Consider problem (4) and (6) in the domain X = (0, 1)2 with a = 106, ~ bðxÞ ¼ p1ffiffi5 ð1; 2ÞT , c = 0. The source 1 2x x 2 4 ffi Þ is the exact solution. It is characterized by a characteristic term f is chosen in such a way that uðxÞ ¼ 12 ð1 tanh 1pffiffiffi 5a pffiffiffi interior layer of thickness Oð aj ln ajÞ around 2x1 x2 ¼ 14. We consider the solutions of the stabilized discrete scheme (29)–(31) without and with shock-capturing with different Ppelements for p 2 {1, 2, 3, 4} on a quasi-uniform mesh with h ¼ 641 . More precisely, we applied the DC/CD variant (44) with b = 0.7, j = 104, l0 = 1 for p 5 1 and l0 = 2 for p = 2. We stopped the simple iteration scheme (50) after 10 iterations if it had not reached the stopping criterion. In Table 1 we compare the error in different norms for the SUPG and the DC/CD scheme for Example 4.2 for p 2 {1, . . ., 4}. First of all, we observe the error reduction with increasing p. Although the SUPG scheme always shows a slightly smaller error, one observes that the errors are of the same magnitude. The larger errors of the shock-capturing scheme can be explained by the slightly increased crosswind-diffusion. We want to study the crosswind-diffusion effects more precisely. In Fig. 2 the crosswind-diffusion parameter s is presented. One can observe very well, that the additional diffusion is located around the layer. As expected, away from the layer no additional crosswind-diffusion is added. Moreover, we present in Fig. 3 contour-plots of the SUPG- and DC/ CD-solutions for p = 1 and p = 4. Indeed we observe a slight additional crosswind-effect of the shock-capturing. Finally, in Fig. 4 we plot the cross-section of the SUPG- and of the DC/CD-solutions in the crosswind-direction at x1 + 2x2 = 1. We observe the classical Gibbs phenomenon of the SUPG solution in the neighborhood of the boundary layer, i.e. there appear significant over- and undershoots of the discrete solution. The oscillations of the SUPG scheme in the characteristic layer are clearly damped with shock-capturing. It is important to observe that the strong gradient
Table 1 The error for the SUPG scheme without and with shock capturing for Example 4.2 with a = 106, c = 0, j = 104, b = 0.7, l0 = 1, h = 1/64 in different norms p
1 2 3 4
k kL2 ðXÞ
jjj Æ jjjE
k kL1 ðXÞ
SUPG
SC/CD
SUPG
SC/CD
SUPG
SC/CD
0.0559 0.0310 0.0216 0.0160
0.0587 0.0315 0.0234 0.0181
0.0589 0.0354 0.0266 0.0211
0.0615 0.0358 0.0280 0.0227
0.4421 0.3734 0.3098 0.2582
0.4985 0.3773 0.3515 0.3174
Fig. 2. The crosswind-diffusion parameter s for Example 4.2 with a = 106, c = 0, j = 104, b = 0.7, l0 = 1, h = 1/64, p = {1, 4}.
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138 SUPG
SC
1
1
0.9
0.9
0.8
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0 0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0 0.1
0.2
0.3
0.4
SUPG 1
0.9
0.9
0.8
0.8
0.7
0.7
0.6
0.6
0.5
0.5
0.4
0.4
0.3
0.3
0.2
0.2
0.1
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.5
0.6
0.7
0.8
SC
1
0 0.1
4135
0.5
0.6
6
0.7
0.8
0 0.1
0.2
0.3
0.4
4
Fig. 3. Contour plot for Example 4.2 with a = 10 , c = 0, j = 10 , b = 0.7, l0 = 1, h = 1/64, p = {1, 4}. On the left-hand side without shock capturing and on the right-hand side with shock capturing.
of the SUPG solution in the layer is still preserved by the shock-capturing scheme. Although not reflected in the present analysis of Section 3, this seems to be a result of the proper construction of the limiter function of the DC/CD-scheme with (44). Moreover, we observe that, for fixed h, the resolution of the strong gradient improves with increasing p. 5. Summary In this paper we extended the a priori analysis in [12] of a hp-version of stabilized Galerkin methods of streamline-diffusion type (SUPG) for diffusion–advection–reaction problems to the case of shock-capturing stabilization. This additional stabilization in crosswind-directions is an important ingredient within the solution of coupled nonlinear flow models, e.g., for turbulence models. In particular, we extended our previous results in [18] by studying carefully the dependence of stabilization parameters on the spectral order and on characteristic local quantities (like the mesh Peclet number). The analysis is supported by some numerical experiments. We improved our previous results in [18] for certain shock-capturing variants of the SUPG-scheme. In particular, we give better upper bounds of the limiter function which are in better agreement with schemes presented in the literature. Nevertheless, there are still open problems: Besides the uniqueness problem, a refined control of the crosswind-error (via the limiter function) would be desirable. Moreover, the convergence of iteration schemes for the arising nonlinear discrete problem is open. It seems to be useful to provide a careful and fair comparison of the present residual-based stabilization approach to the concept of edge-stabilization [2,4]. Moreover, it is necessary to extend the present concept to discontinuous Galerkin methods (eventually patchwise) in order to simplify the application of hp-adaptivity and domain decomposition techniques. In paper [19], we extend the present approach to the linearized incompressible Navier–Stokes problem, thus completing the analysis of auxiliary problems appearing within the iterative solution of coupled problems as shown in Section 1.
4136
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138 1.2
SC SUPG
1
1.2 1
0.8
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0
–0.2 0
0.05
0.1
0.15
0.2
0.25
1.2
–0.2
0
0.05
0.1
0.15
0.2
0.25
1.2
SC SUPG
1
0.8
0.6
0.6
0.4
0.4
0.2
0.2
0
0
0.05
0.1
0.15
0.2
0.25
SC SUPG
1
0.8
–0.2 0
SC SUPG
–0.2 0
0.05
0.1
0.15
0.2
0.25
6
Fig. 4. SUPG-FEM without and with SC-stabilization for Example 4.2 with a = 10 , c = 0 and h ¼ 641 ; p 2 f1; 2; 3; 4g.
Appendix Here we summarize some basic estimates for the SUPG-scheme applied to the linear advection–diffusion–reaction problems, see Section 2. Proof of Lemma 2.1. Setting u = v in (14) and integrating the advective term by parts together with (5), we obtain X pffiffiffi pffiffiffi 2 2 2 ðk arvkðL2 ðT ÞÞd þ k cvkL2 ðT Þ þ dT k~ b rvkL2 ðT Þ þ dT ðr PT ðarvÞ þ cv; ~ b rvÞT Þ. As ðv; vÞ ¼ T 2Th
We need the following auxiliary estimates, which use the local inverse estimates (12) together with the assumptions (18): pffiffiffi l p2 1=2 dT ðr PT ðarvÞ; ~ b rvÞT 6 dT inv T kakL1 ðT Þ k arvkðL2 ðT ÞÞd k~ b rvkL2 ðT Þ hT pffiffiffiffiffi 1 pffiffiffi b rvkL2 ðT Þ ; 6 pffiffiffi k arvkðL2 ðT ÞÞd dT k~ 2 pffiffiffi 1=2 b rvÞT 6 dT kckL1 ðT Þ k cvkL2 ðT Þ k~ b rvkL2 ðT Þ dT ðcv; ~ p ffiffiffiffiffi 1 pffiffiffi 6 pffiffiffi k cvkL2 ðT Þ dT k~ b rvkL2 ðT Þ . 2 We combine both estimates by using Youngs inequality pffiffiffi pffiffiffi 2 1 X ~ 1 2 2 2 2 As ðv; vÞ P jjjvjjj dT kb rvkL2 ðT Þ þ k arvkðL2 ðT ÞÞd þ k cvkL2 ðT Þ P jjjvjjj . 2 T 2Th 2
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138
4137
We now prove the continuity result: The symmetric part of As(u, vh) is bounded by jjjujjj jjjvhjjj. The remaining terms can be bounded as follows: ð~ b ru; vh Þ ¼ ðu; ~ b rvh Þ ( ) !1=2 2 X bkðL1 ðT ÞÞd 1 k~ 2 6 min ; kukL2 ðT Þ dT aT T 2Th X
dT ðr PT ðaruÞ þ cu; ~ b rvh ÞT 6
T 2Th
X
! 1=2 X pffiffiffi 2 2 k arvh kL2 ðT Þ þ dT k~ b rvh kL2 ðT Þ ; T 2Th
!1=2 dT k r PT ðaruÞ þ
cuk2L2 ðT Þ
T 2Th
X
!1=2 dT k~ b rvh k2L2 ðT Þ
.
T 2Th
Summarizing all terms and dividing by jjjvhjjj, we find the assertion.
h
Proof of Lemma 2.2. For arbitrary vh 2 V ph we set f := uh vh and g := u vh. Lemmas 2.1 and 2.2 imply 0 " #1=2 1 X X 1 2 2 Ajjjfjjj. jjjfjjj 6 As ðf; fÞ ¼ As ðg; fÞ þ dT ððL b LÞu; ~ b rfÞT 6 @Qs ðgÞ þ dT kðL b LÞukL2 ðT Þ 2 T 2Th T 2Th This estimate together with the triangle inequality conclude the proof.
h
Proof of Corollary 2.1. According to Lemma 2.2, the right-hand side of (21) is bounded by ( ) 2 X bkðL1 ðT ÞÞd 1 k~ 2 2 2 jjju U h jjj ¼ C jjju I h;p ujjj þ min ; ku I h;p ukL2 ðT Þ d a T T T 2Th þ
X
dT k r PT ðarðu I h;p uÞÞ þ cðu
2 I h;p uÞkL2 ðT Þ
T 2Th
þ
X
! 2 b dT kð L LÞukL2 ðT Þ .
ð51Þ
T 2Th
Now we estimate term by term using (18), and the interpolation estimate (22). The first one can be estimated by 2lT 2 X h2 hT kakL1 ðT Þ þ kckL1 ðT Þ T2 þ dT k~ bk2ðL1 ðT ÞÞd kuk2H kT ðT Þ . jjju I h;p ujjj2 6 3C 2I 2k T 2 p p T T T 2Th
ð52Þ
The second term in (51) is bounded in a similar way. Using the inverse inequality (12) and the definition of PT the third term of (51) can be bounded by ! T X X p4T hT2lT 2 h2l 2 2 2 2 2 T dT k r PT ðarðu I h;p uÞÞ þ cðu I h;p uÞkL2 ðT Þ 6 2C I dT linv kakL1 ðT Þ 2 2kT 2 þ kckL1 ðT Þ 2kT kuk2H kT ðT Þ . hT p T pT T 2Th T 2Th ð53Þ The last term of (51) can be estimated as follows: X X dT kð b L LÞuk2L2 ðT Þ ¼ dT kr ðaruÞ r PT ðaruÞk2L2 ðT Þ T 2Th
T 2Th
62
X
2 2 dT kr ðaruÞ r I Th;p ðaruÞkL2 ðT Þ þ kr I Th;p ðaruÞ r PT ðaruÞkL2 ðT Þ ;
T 2Th
I Th;p
is the local Lagrange interpolation. Using the approximation property (22) and the inverse inequality (12) yields where for the consistency error ! 2lT 4 4 X X 2 2 2 2 hT 2 pT T b dT kð L LÞukL2 ðT Þ 6 2 dT linv 2 kI h;p ðaruÞ PT ðaruÞkL2 ðT Þ þ C I 2kT 4 karukH kT 1 ðT Þ hT pT T 2Th T 2Th ! 2lT 4 T 2 X p4 h2l 2 hT T 62 dT l2inv C 2I T2 2k þ C ð54Þ kak2W kT 1;1 ðT Þ kuk2H kT ðT Þ ; I 2 2k T T 4 hT p T pT T 2Th where we have used the definition of PT and, again, the approximation property (22) in the last step. Summarizing (51)– (54), we obtain 2
jjju U h jjj 6 C
X h2ðlT 1Þ T T 2Th
2ðk T 1Þ
pT
2
kukH kT ðT Þ M T ðhT ; pT ; dT ; a; ~ b; cÞ
4138
G. Lube, G. Rapin / Comput. Methods Appl. Mech. Engrg. 195 (2006) 4124–4138
with (
2 bkðL1 ðT ÞÞd 1 k~ þ min ; dT aT
)
h2T þ dT k~ bk2ðL1 ðT ÞÞd p2T p2 þ kakL1 ðT Þ þ dT l2inv p2T kak2L1 ðT Þ þ ð1 þ l2inv p2T Þkak2W kT 1;1 ðT Þ T2 . hT
M T :¼ ðkckL1 ðT Þ þ
h2 dT kck2L1 ðT Þ Þ T2 pT
Applying the assumption (18) yields the assertion.
h
References [1] A.N. Brooks, T.J.R. Hughes, Streamline upwind Petrov–Galerkin formulation for convection dominated flows with particular emphasis on the incompressible Navier–Stokes equations, Comput. Methods Appl. Mech. Engrg. 32 (1982) 199–259. [2] E. Burman, A. Ern, Nonlinear diffusion and discrete maximum principle for stabilized Galerkin approximations of the convection–diffusion–reaction equation, Comput. Methods Appl. Mech. Engrg. 191 (2002) 3833–3855. [3] E. Burman, A. Ern, Stabilized Galerkin approximation of convection–diffusion–reaction equations: discrete maximum principle and convergence, Math. Comput., accepted for publication. [4] E. Burmann, P. Hansbo, Edge stabilization for Galerkin approximations of convection–diffusion problems, Comput. Methods Appl. Mech. Engrg. 193 (2004) 1437–1453. [5] 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. [6] R. Codina, O. Soto, Finite element implementation of two-equation and algebraic stress turbulence models for steady incompressible flows, Int. J. Numer. Meth. Fluids 90 (3) (1999) 309–343. [7] E.G. Do Carmo, G. Alvarez, A new stabilized finite element formulation for scalar convection–diffusion problems: the streamline and approximate upwind/Petrov–Galerlin method, Comput. Methods Appl. Mech. Engrg. 192 (2003) 3379–3396. [8] A. Ern, J.-L. Guermond, Theory and Practice of Finite Elements, Springer, 2004. [9] A.C. Galeao, E.G. Dutra do Carmo, A consistent approximate upwind Petrov–Galerkin method for convection-dominated problems, Comput. Methods Appl. Mech. Engrg. 68 (1988) 83–95. [10] A.C. Galeao, R.C. Almeida, S.M.C. Malta, A.F.D. Loula, Finite element analysis of convection dominated reaction–diffusion problems, Appl. Numer. Math. 48 (2004) 205–222. [11] B.-Y. Guo, H.-P. Ma, E. Tadmor, Spectral vanishing viscosity method for nonlinear conservation laws, SIAM J. Numer. Anal. 39 (4) (2001) 1254– 1268. [12] P. Houston, E. Su¨li, Stabilized hp-finite element approximations of partial differential equations with nonnegative characteristic form, Computing 66 (2001) 99–119. [13] T.J.R. Hughes, M. Mallet, A. Mizukami, A new finite element formulation for computational fluid dynamics: II. Beyond SUPG, Comput. Methods Appl. Mech. Engrg. 54 (1986) 341–355. [14] T.J.R. Hughes, G. Scovazzi, L.P. Franca, Multiscale and stabilized methods, in: E. Stein et al. (Eds.), Encyclopedia of Computational Mechanics, vol. 3, John Wiley and Sons, 2004. [15] C. Johnson, U. Na¨vert, Pitka¨ranta: finite element methods for linear hyperbolic problems, Comput. Methods Appl. Mech. Engrg. 45 (1984) 25–38. [16] C. Johnson, A.H. Schatz, Wahlbin: crosswind smear and pointwise errors in streamline diffusion finite element methods, Math. Comput. 49 (1987) 25–38. [17] T. Knopp, G. Lube, R. Gritzki, M. Ro¨sler, Iterative substructuring methods for incompressible non-isothermal flows and its application to indoor air flow simulation, Int. J. Numer. Meth. Fluids 40 (2002) 1527–1538. [18] T. Knopp, G. Lube, G. Rapin, Stabilized finite element methods with shock capturing for advection–diffusion problems, Comput. Methods Appl. Mech. Engrg. 191 (2002) 2997–3013. [19] G. Lube, G. Rapin, Residual-based stabilized higher-order FEM for a generalized Oseen problem, Math. Mod. Meth. Appl. Sci., accepted for publication. [20] U. Na¨vert, A Finite Element Method for Convection–Diffusion Problems, Ph.D. Thesis, Chalmers University of Technology, Go¨teborg, 1982. [21] C. Schwab, p- and hp-finite element methods. Theory and Applications to Solid and Fluid Mechanics, Oxford University Press, Oxford, 1998. [22] Y.-T. Shih, H.C. Elman, Iterative methods for stabilized discrete convection-diffusion problems, IMA J. Numer. Anal. 20 (2000) 333–358. [23] A. Szepessy, Convergence of a shock-capturing streamline diffusion finite element method to a scalar conservation law in two dimensions, Math. Comput. 53 (1989) 527–545. [24] A. Toselli, hp-finite element approximations on non-matching grids for partial differential equation with non-negative characteristic form, M2AN 37 (2003) 91–115. [25] G. Zhou, Local pointwise error estimates for the streamline diffusion method applied to nonstationary hyperbolic problems, East-West J. Numer. Math. 3 (1995) 217–235. [26] G. Zhou, R. Rannacher, Pointwise superconvergence of the streamline diffusion finite element method, Numer. Methods Part. Diff. Equat. 12 (1996) 123–145.