Reliable and efficient a posteriori error estimates of DG methods for a frictional contact problem Fei Wang1 and Weimin Han2 4 Abstract. A posteriori error estimators are studied for discontinuous Galerkin 1 0 methods for solving a frictional contact problem, which is a representative elliptic 2 variational inequality of the second kind. The estimators are derived by relating n a the error of the variational inequality to that of a linear problem. Reliability and J efficiency of the estimators are shown. 0 2 Keywords. Elliptic variational inequality, discontinuous Galerkin method, a pos- ] teriori error estimators, reliability, efficiency A N AMS Classification. 65N15, 65N30, 49J40 . h t a m 1 Introduction [ 1 For more than three decades, adaptive finite element method (AFEM) has been an active v 1 research field in scientific computing. As an efficient numerical approach, it has been widely 1 used for solving a variety of differential equations. Each loop of AFEM consists of four steps, 0 5 . 1 Solve → Estimate → Mark → Refine. 0 4 That is, in each loop, we first solve the problem on an mesh, then use a posteriori error 1 v: estimators to mark those elements to be refined, and finally, refine the marked elements and Xi geta newmesh. Wecancontinue thisprocess until theerror satisfies certainsmallness criterion. r The adaptive finite element method can achieve high accuracy with lower memory usage and a less computation time. A posteriori error estimators are computable quantities that indicate the contribution of error on each element to the global error. They are used in adaptive algorithms to indicate which elements need to be refined or coarsened. To capture the true error as precisely as possible, they should have two properties: reliability and efficiency ([1, 4]). Hence, obtaining 1Department of Mathematics, Pennsylvania State University, University Park, PA 16802, USA. School of Mathematics and Statistics, Huazhong University of Science and Technology, Wuhan 430074, China. Email: [email protected] 2Department of Mathematics & Programin Applied Mathematical and Computational Sciences, University of Iowa, Iowa City, Iowa 52242,USA. 1 reliable and efficient error estimators is the key for successful adaptive algorithms. A variety of different a posteriori error estimators have been proposed and analyzed. Many error estimators can be classified as residual type or recovery type ([1, 4]). Various residual quantities are used to capture lost information going from u to u , such as residual of the equation, residual from h derivative discontinuity and so on. Another type of error estimators is gradient recovery, i.e., ||G(∇u )−∇u || is used to approximate ||∇u−∇u ||, where a recovery operator G is applied h h h to the numerical solution u to rebuild the gradient of the true solution u. A posteriori error h analysis have been well established for standard finite element methods for solving linear partial differential equations, and we refer the reader to [1, 4, 28]. Due to the inequality feature, it is more difficult to develop a posteriori error estimators for variational inequalities (VIs). However, numerous articles can be found on a posteriori error analysis of finite element methods for the obstacle problem, which is an elliptic variational inequality (EVI) of the first kind, e.g., [5, 15, 22, 24, 27, 32]. In [11], Braess demonstrated that a posteriori error estimators for finite element solutions of the obstacle problem can be derived by applying a posteriori error estimates for an associated linear elliptic problem. For VIs of the second kind, in [7, 8, 9, 10], the authors studied a posteriori error estimation and established a framework through the duality theory, but the efficiency was not completely proved. In [29], the ideas in [11] were extended to give a posteriori error analysis for VIs of the second kind. Moreover, a proof was provided for the efficiency of the error estimators. In recent years, thanks to the flexibility in constructing feasible local shape function spaces and the advantage to capture non-smooth or oscillatory solutions effectively, discontinuous Galerkin (DG) methods have been widely used for solving various types of partial differential equations. When applying h-adaptive algorithm with standard finite element methods, one needs to choose the mesh refinement rule carefully to maintain mesh conformity and shape regularity. In particular, hanging nodes are not allowed without special treatment. For discon- tinuous Galerkinmethods, theapproximate functions areallowed to bediscontinuous acrossthe element boundaries, so general meshes with hanging nodes and elements of different shapes are accepted. Advantages of DG methods include the flexibility of mesh-refinements and construc- tion of local shape function spaces (hp-adaptivity), and the increase of locality in discretization, which is of particular interest for parallel computing. A historical account of DG methods’ de- velopment can be found in [16]. In [2, 3], Arnold et al. established a unified error analysis of nine DG methods for elliptic problems and several articles provided a posteriori error analysis of DG methods for elliptic problems (e.g. [6, 12, 14, 21, 23, 25]). Carstensen et al. presented a unified approach to a posteriori error analysis for DG methods in [13]. In [30], the authors extended ideas of the unified framework about DG methods for elliptic problems presented in [3] to solve the obstacle problem and a simplified frictional contact problem, and obtained a priori error estimates, which reach optimal order for linear elements. In [31], reliable a posteri- ori error estimators of the residual type were derived for DG methods for solving the obstacle problem, and efficiency of the estimators is theoretically explored and numerically confirmed. 2 A posteriori error analysis of DG methods for the obstacle problem was also studied in [20]. In this paper, we study a posteriori error estimates of DG methods for solving a frictional contact problem. The paper is organized as follows: in Section 2 we introduce a frictional contact problem and the DG schemes for solving it. Then we derive a reliable residual type a posteriori error estimators for the DG methods of a frictional contact problem in Section 3. Finally, we prove efficiency of the proposed error estimators in Section 4. 2 A frictional contact problem and DG formulations 2.1 A frictional contact problem . Let Ω ⊂ Rd (d = 2,3) be an open bounded domain with Lipschitz boundary Γ that is divided into two mutually disjoint parts, i.e., Γ = Γ ∪Γ . Here Γ is a relatively closed subset 1 2 1 of Γ, and Γ = Γ\Γ . Given f ∈ L2(Ω) and a constant g > 0, the frictional contact problem is: 2 1 find u ∈ V = H1 (Ω) := {v ∈ H1(Ω) : v = 0 a.e. on Γ } such that Γ1 1 a(u,v−u)+j(v)−j(u) ≥ (f,v −u) ∀v ∈ V, (2.1) where (·,·) denotes the L2 inner product in the domain Ω and a(u,v) = ∇u·∇vdx+ uvdx, ZΩ ZΩ j(v) = g|v|ds. ZΓ2 The frictional contact problem is an example of elliptic variational inequalities of the second kind and has a unique solution u ∈ V ([18, 19]). Moreover, there exists a unique Lagrange multiplier λ ∈ L∞(Γ ) such that 2 a(u,v)+ gλvds = (f,v) ∀v ∈ V, (2.2) ZΓ2 |λ| ≤ 1, λu = |u| a.e. on Γ . (2.3) 2 From (2.2) and (2.3), we know that the solution u of (2.1) is the weak solution of the following boundary value problem −△u+u = f in Ω, u = 0 on Γ , 1 ∇u·n = −gλ on Γ , 2 3 where n is the unit outward normal vector. For any v ∈ V, set ℓ(v) = f vdx− gλvds. ZΩ ZΓ2 Then we have by (2.2) a(u,v) = ℓ(v) ∀v ∈ V. (2.4) Similar with the argument in [29], given a triangulation T of Ω, for a Lipschitz subdomain h ω ⊂ Ω, define a (v,w) := (∇v ·∇w +vw)dx ω,h KX∈ThZω∩K and kvk := a (v,v)1/2. 1,ω,h ω,h Then define |λ| := sup gλvds : v ∈ H1(ω), kvk = 1 , (2.5) ∗,γ,h h 1,ω,h (cid:26)Zγ (cid:27) where γ ⊂ ∂ω ∩Γ is a measurable subset and H1(ω) = {v ∈ L2(ω) : v| ∈ H1(K ∩ω)}. If 2 h K∩ω ω = Ω and γ = Γ , the subscript ω and γ are omitted. We have 2 |λ| = kwk , (2.6) ∗,γ,h 1,ω,h where w ∈ H1(ω) is the solution of the following auxiliary equation h a (w,v) = g λ vds ∀v ∈ H1(ω). (2.7) ω,h h Zγ The formula (2.6) can be proved by an argument similar to that found in [29]. 2.2 Discontinuous Galerkin formulations First, we introduce some notations. Let {T } be a family of triangulations of Ω such that h the minimal angle condition is satisfied. For a triangulation T , let E be the set of all edges, h h Ei ⊂ E the set of all interior edges, Eb := E \Ei the set of all boundary edges, E0 ⊂ E the set h h h h h h h of all edges not lying on Γ , E1 := E0\Ei, E2 := E \E0, and define E(K) as the set of sides of 2 h h h h h h K. Let h = diam(K) for K ∈ T , h = length(e) for e ∈ E , and N denote the set of nodes K h e h h of T . For any element K ∈ T , define the patch set ω := ∪{T ∈ T , T ∩K 6= Ø}, and for h h K h any edge e shared by two elements K+ and K−, define ω := K+ ∪ K−. For a scalar-valued e function v and a vector-valued function q, let vi = v| , qi = q| , and ni = n| be the ∂Ki ∂Ki ∂Ki 4 unit normal vector external to ∂Ki with i = ±. Define the average {·} and the jump [·] on an interior edge e ∈ Ei as follows: h 1 {v} = (v+ +v−), [v] = v+n+ +v−n−, 2 1 {q} = (q+ +q−), [q] = q+ ·n+ +q− ·n−. 2 For a boundary edge e ∈ Eb, we let h [v] = vn, {q} = q, where n is the outward unit normal. Let us define the following linear finite element spaces V = {v ∈ L2(Ω) : v | ∈ P (K) ∀K ∈ T }, h h h K 1 h W = {w ∈ [L2(Ω)]2 : w | ∈ [P (K)]2 ∀K ∈ T }. h h h K 1 h We denote by ∇ the broken gradient whose restriction on each element K ∈ T is equal to ∇. h h Define some seminorms and norms by the following relations: kvk2 = v2dx, |v|2 = k∇vk2 , kvk2 = v2ds, K 1,K K e ZK Ze kvk2 = kvk2 , |v|2 = |v|2 , kvk2 = kvk2 +|v|2 . 0,h K 1,h 1,K 1,h 0,h 1,h KX∈Th KX∈Th Throughout this paper, “. ···” stands for “≤ C···”, where C denotes a generic positive constant dependent on the minimal angle condition but not on the element sizes, which may take different values at different occurrences. Now, let us introduce the Discontinuous Galerkin methods for solving the variational in- equality (2.1). Here, we take the local DG method (LDG) as an example to show how to derive a posteriori error estimators of DG methods for solving the frictional contact problem (2.1). The derivation and analysis for the LDG method in this paper can be extended to other DG methods studied in [30]. The LDG method ([17]) for solving the frictional contact problem is to find u ∈ V such that h h B (u ,v −u )+j(v )−j(u ) ≥ (f,v −u ) ∀v ∈ V , (2.8) h h h h h h h h h h where B (u,v) = (∇ u·∇ v +uv)dx− [u]·{∇ v}ds− {∇ u}·[v]ds h h h h h ZΩ ZEh0 ZEh0 − β ·[u][∇ v]ds− [∇ u]β ·[v]ds h h Ei Ei Z h Z h +(r ([u])+l(β ·[u]),r ([v])+l(β ·[v]))+αj(u,v). (2.9) 0 0 0 5 Here β ∈ [L2(Ei)]2 is a vector-valued function which is constant on each edge of Ei, and h h αj(u,v) = η[u]·[v]ds is the penalty term with the penalty weighting function η : E0 → R 0 E0 h h given by η h−1 on each e ∈ E0, η being a positive number on e. For any w ∈ W , the lifting eRe h e h h operators r : [L2(E0)]2 → W and l : L2(Ei) → W are defined by 0 h h h h r (q)·w dx = − q ·{w }ds, l(v)·w dx = − v[w ]ds ∀w ∈ W . (2.10) 0 h h h h h h ZΩ ZEh0 ZΩ ZEhi The bilinear form B is continuous and elliptic with respect to certain DG-norm, and there- h fore, in particular, the discrete problem has a unique solution u ∈ V (see [3, 30]). Similar to h h the continuous problem, there exists a unique Lagrange multiplier λ ∈ L∞(Γ ) such that ([19]) h 2 B (u ,v )+ gλ v ds = (f,v ) ∀v ∈ V , (2.11) h h h h h h h h ZΓ2 |λ | ≤ 1, λ u = |u | a.e. on Γ . (2.12) h h h h 2 For any v ∈ V , let h h ℓ (v ) = (f,v )− gλ v ds. h h h h h ZΓ2 Then (2.11) becomes B (u ,v ) = ℓ (v ) ∀ v ∈ V . (2.13) h h h h h h h For any v ∈ V, we know that [u] = 0 and [v] = 0 on e ∈ E0. Then we have from (2.2) that h B (u,v) = a(u,v) = ℓ(v) ∀v ∈ V (2.14) h Obviously, u is also the finite element approximation of the solution z ∈ V of the linear h problem: B (z,v) = ℓ (v) ∀v ∈ V, (2.15) h h which is the weak formulation of the boundary value problem −∆z +z = f in Ω, (2.16) z = 0 on Γ , 1 ∂z = −gλ on Γ . h 2 ∂n 2.3 A bridge between u − u and u − z h h Next, we relate the error e := u −u to u −z, namely, h h 1/2 kek +|λ−λ | . ku −zk + h−1k[u ]k2 . (2.17) 1,h h ∗,h h 1,h e h e eX∈Eh0 6 Then we will use this relation to derive a posteriori error estimators for DG solutions of the frictional contact problem by utilizing a posteriori error estimators of the related linear ellip- tic problem (2.16). Note that a similar approach can be applied to other elliptic variational inequalities of the second kind. To derive the inequality (2.17), we first define a continuous piecewise linear function in V ∩ H1 (Ω), whose value is close to the numerical solution. For any given v ∈ V , written h Γ1 h h v = 3 α(j)φ(j), where φ(j), 1 ≤ j ≤ 3, are the linear basis functions corresponding h K∈Th j=1 K K K to the three vertices of K, we construct a function χ ∈ V ∩H1 (Ω) as follows: At every interior P P h Γ1 node and the nodes on Γ of the conforming mesh T , the value of χ is set to be the average of 2 h the values of v computed from all the elements sharing that node, and χ = 0 at the boundary h nodes on Γ . For each ν ∈ N , let ω = {K ∈ T : ν ∈ K} and denote its cardinality by |ω |, 1 h ν h ν which is bounded by a constant depending only on the minimal angle condition of the mesh. To each node ν, the associated basis function φ(ν) is given by suppφ(ν) = K, φ(ν)| = φ(j) for x(j) = ν. K K K K[∈ων Then we define χ ∈ V ∩H1 (Ω) by h Γ1 1 χ = β(ν)φ(ν), where β(ν) = α(j) if ν ∈ N and ν 6∈ Γ . (2.18) |ω | K h 1 νX∈Nh ν xX(j)=ν K For nonconforming meshes, let N0 be the set of all hanging nodes. Then we construct χ from h v same as conforming mesh case on all the nodes ν ∈ N \N0. For an upper bound of the h h h error v −χ, we quote a result from [21] (which is Theorem 2.2 there for conforming meshes; h the same result also holds for nonconforming meshes, which is Theorem 2.3 in [21]). Lemma 2.1 Let T be a conforming triangulation. Then for any v ∈ V , we can construct a h h h continuous function χ ∈ V ∩H1 (Ω) from v , such that h Γ1 h kv −χk2 ≤ C h1−2ik[v ]k2, i = 0,1, (2.19) h i,K e h e KX∈Th eX∈Eh0 where the constant C is independent of mesh size and v but which may depend on the lower h bound of the minimal angle of the elements in T . h Now, let us derive the inequality (2.17). From (2.14) and (2.15), for all v ∈ V, we have B (u −u,v) = B (u −z,v)+B (z −u,v) = B (u −z,v)+ g(λ−λ )vds. h h h h h h h h ZΓ2 7 By the definition (2.9) and noticing [v] = 0 on each e ∈ E0, the above equation becomes h a(e,v)− [e]·{∇ v}ds− β ·[e][∇ v]ds h h E0 Ei Z h Z h e = a(u −z,v)− [u −z]·{∇ v}ds h h h E0 Z h −e β ·[u −z][∇ v]ds+ g(λ−λ )vds, h h h ZEhi ZΓ2 where a(u,v) = (∇ u·∇ v +uv)dx. h h ZΩ Then, e a(e,v) = a(u −z,v)− [u−z]·{∇ v}ds− β ·[u−z][∇ v]ds+ g(λ−λ )vds. h h h h ZEh0 ZEhi ZΓ2 Noete that [ue−z] = 0 on each e ∈ E0. We have h a(e,v) = a(u −z,v)+ g(λ−λ )vds. (2.20) h h ZΓ2 Let χ ∈ V ∩H1 (Ω) be the function constructed from u , satisfying (2.19) for v = u . Taking h Γ1 e e h h h v := χ−u = χ−u +u −u in (2.20) and using Cauchy-Schwarz inequality, we have h h kek2 ≤ ku −zk (kχ−u k +kek )+kek kχ−u k + g(λ−λ )(χ−u)ds 1,h h 1,h h 1,h 1,h 1,h h 1,h h ZΓ2 = kek (ku −zk +kχ−u k )+ku −zk kχ−u k 1,h h 1,h h 1,h h 1,h h 1,h + g(λ−λ )(χ−u)ds h ZΓ2 1 1 ≤ kek2 + (ku −zk +kχ−u k )2 +ku −zk kχ−u k 2 1,h 2 h 1,h h 1,h h 1,h h 1,h + g(λ−λ )(χ−u)ds. (2.21) h ZΓ2 Note that by (2.3) and (2.12), we have g(λ−λ )(u −u)ds = g λ u ds− g λ uds− g λ u ds+ g λ uds h h h h h h ZΓ2 ZΓ2 ZΓ2 ZΓ2 ZΓ2 ≤ g |u |ds− g |u|ds− g |u |ds+ g |u|ds = 0. h h ZΓ2 ZΓ2 ZΓ2 ZΓ2 and g(λ−λ )(χ−u )ds ≤ |λ−λ | kχ−u k h h h ∗,h h 1,h ZΓ2 1 ≤ ǫ|λ−λ |2 + kχ−u k2 . h ∗,h 4ǫ h 1,h 8 Hence, kek2 . ku −zk2 +kχ−u k2 +ǫ|λ−λ |2 . (2.22) 1,h h 1,h h 1,h h ∗,h Recalling (2.6), we have |λ−λ | = ku−zk ≤ kek +ku −zk . h ∗,h 1,h 1,h h 1,h Then, we obtain the following result kek +|λ−λ | . ku −zk +kχ−u k . 1,h h ∗,h h 1,h h 1,h Using (2.19) to bound kχ−u k , the above inequality can be rewritten as h 1,h 1/2 kek +|λ−λ | . ku −zk + h−1k[u ]k2 . (2.23) 1,h h ∗,h h 1,h e h e eX∈Eh0 The relation (2.23) serves as a starting point for derivation of reliable and efficient error estimators of DG methods for a frictional contact problem. In this paper, we focus on the derivation and analysis of residual type error estimators derived from the inequality (2.23). A similar approach can also be applied to recovery type error estimators. 3 Reliable residual-type estimators Now we follow the ideas in [29] to obtain a posteriori error estimators of DG methods for solving the frictional contact problem. The detailed derivation and analysis of a posteriori error estimators is given for the LDG method [17]. For other DG methods discussed in [30], similar results could be obtained by similar arguments. To bound the first term ku −zk , we recall one result in [13]. Note that the a posteriori h 1,h erroranalysisin[13]wasonlyforthePoissonproblemwithhomogenousDirichletboundarycon- dition, but it is easy to extend the result to general elliptic problems with Neumann boundary conditions. For the second-order elliptic problem ∂u −∆u+u = f in Ω, u = 0 on Γ , = g on Γ , 1 2 ∂n rewrite it as the first order system ∂u p = ∇u, −∇·p+u = f in Ω, u = 0 on Γ , = g on Γ . (3.1) 1 2 ∂n 9 Then the DG formulation for this problem is p ·τ dx = − u ∇ ·τ dx+ uˆ n ·τ ds ∀τ ∈ W , (3.2) h h h h h h K h h h ZΩ ZΩ KX∈ThZ∂K (p ·∇ v +u v )dx = f v dx+ pˆ ·n v ds ∀v ∈ V , (3.3) h h h h h h h K h h h ZΩ ZΩ KX∈ThZ∂K where uˆ and pˆ are numerical fluxes. Different choices of the numerical fluxes lead to different h h DG methods. The following theorem (see [13]) holds for the LDG method and other methods discussed in [3]. Theorem 3.1 Assume u ∈ H1 (Ω) and p ∈ W := [L2(Ω)]2 are the solution of the problem Γ1 (3.1), and u ∈ V and p ∈ W are the solution of the problem (3.2)–(3.3). Then, h h h h kp−p k ≤ C(η +ζ ), h ∗ ∗ where η2 := h2 kdivp −u +fk2 + h k[p ]k2 + h kp ·n−gk2, ∗ K h h K e h e e h e KX∈Th eX∈Ehi eX∈Eh2 ζ2 := h−1k[u ]k2 ∗ e h e eX∈Eh0 and C is a mesh-sizeindependentconstantwhich depends only on the domain Ω and the minimal angle condition. From the relation between p and u ([3, 13]), we deduce the following result. h h Corollary 3.2 With the same notation as in Theorem 3.1, we have k∇u−∇ u k ≤ C(η +ζ ), h h ∗ where η2 := h2 k∆u −u +fk2 + h k[∇ u ]k2 + h k∇ u ·n−gk2. K h h K e h h e e h h e KX∈Th eX∈Ehi eX∈Eh2 Proof. By [26, Lemma 7.2], kr ([v ])k2 ≤ C h−1k[v ]k2, kl(β ·[v ])k2 ≤ C h−1k[v ]k2, ∀v ∈ V . 0 h e h e h e h e h h eX∈Eh0 eX∈Ehi 10