Discontinuous Galerkin methods for a contact problem with Tresca friction arising in linear elasticity

Accepted Manuscript Discontinuous Galerkin methods for a contact problem with Tresca friction arising in linear elasticity

Kamana Porwal

PII: DOI: Reference:

S0168-9274(16)30213-6 http://dx.doi.org/10.1016/j.apnum.2016.10.012 APNUM 3118

To appear in:

Applied Numerical Mathematics

Received date: Revised date: Accepted date:

22 October 2015 9 September 2016 24 October 2016

Please cite this article in press as: K. Porwal, Discontinuous Galerkin methods for a contact problem with Tresca friction arising in linear elasticity, Appl. Numer. Math. (2016), http://dx.doi.org/10.1016/j.apnum.2016.10.012

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Discontinuous Galerkin methods for a contact problem with Tresca friction arising in linear elasticity Kamana Porwal∗ Department of Mathematics and Center for Computation and Technology, Louisiana State University, Baton Rouge, LA 70803

Abstract In this article, we propose and analyze discontinuous Galerkin (DG) methods for a contact problem with Tresca friction for the linearized elastic material. We derive a residual based a posteriori error estimator for the proposed class of DG methods. The reliability and the efficiency of a posteriori error estimator is shown. We further investigate a priori error estimates under the minimal regularity assumption on the exact solution. An important property shared by a class of DG methods, allow us to carry out the analysis in a unified framework. Numerical experiments are reported to illustrate theoretical results. Keywords: finite element, discontinuous Galerkin, a posteriori error estimate, frictional contact problem, variational inequalities, linear elasticity 2000 MSC: 65N30, 65N15, 65N12

1. Introduction In the last few decades, different methods have been studied for the numerical approximation of variational inequalities (VIs), which form an important class of nonlinear PDEs with the applications in mechanics and physical sciences. The elliptic variational inequalities can be studied in two categories, VIs of the first kind and VIs of the second kind. In VIs of the first kind, the variational formulation is posed over a closed and convex set, while in the VIs of second kind, the nonlinearity in the variational formulation is due to the presence of a nondifferentiable term in the bilinear form. In recent years, a lot of research has been done on the numerical analysis of VIs of the first kind. We refer to [18, 19, 32, 31, 57, 4, 44, 51, 54] for a priori analysis and [5, 7, 12, 25, 34, 35, 36, 40, 41, 52, 46, 53, 58] for a posteriori analysis of VIs of the first kind. There are only handful works on the numerical approximation of VIs of the second kind. A classical nonconforming finite element method for a fourth order variational inequality of the second kind ∗ Corresponding

author Email address: [email protected] (Kamana Porwal)

Preprint submitted to Applied Numerical Mathematics

October 25, 2016

has been studied in [39]. The article [37] is devoted to the analysis of continuous interior penalty method for a fourth order variational inequality of the second kind. In [1], a dual-dual variational inequality formulation is studied for frictional contact problems. For a posteriori analysis of VIs of the second kind, we refer to articles [8, 11, 9, 10, 29, 38, 56, 48, 22]. A recovery based a posteriori error estimator of conforming finite element method have been developed in [8, 11, 10], while [56] studied a recovery as well as residual based a posteriori error estimator. A reliable a posterirori error estimator of conforming finite element method for variational inequalities of the second kind is derived in [38] and a goal oriented a posteriori error analysis for frictional contact problems is established in [48]. The articles [29, 22] studied adaptive hp finite element method over quadrilateral and simplical elements, respectively. The analysis in [22] is presented in unified framework compassing variational inequalities of the first and second kind, therein a priori analysis of hp finite element methods is also discussed. In [9], a residual based a posteriori error estimator of conforming finite element method for a frictional contact problem arising in linear elasticity has been derived. In this article, we study discontinuous Galerkin (DG) methods for a frictional contact problem in linear elasticity, which is a variational inequality of the second kind. Discontinuous Galerkin methods are well known for their high order accuracy, flexibility for hp-adaptivity and high suitability for parallel computing. DG method was first introduced by Reed and Hill [49] as a technique to solve neutron transport equation. Since then, DG methods have been studied for several linear and nonlinear PDEs. The detailed analysis of the DG methods can be find in [3]. In this article, we propose and analyze the DG methods for a frictional contact problem in linear elasticity. We derive a reliable and efficient residual based a posteriori error estimator. In addition to that, we also establish a priori estimates assuming the minimal regularity of the exact solution. A common property shared by many DG methods, allows us to study them in a unified framework. Let Ω ⊂ R2 be a bounded polygonal domain with boundary ∂Ω = Γ partitioned into three parts ΓD , ΓN and ΓC which are open and mutually disjoint sets with meas(ΓD ) > 0. The linearized strain tensor ε and stress tensor σ are defined, respectively, as ε(u) =

 1 ∇u + ∇uT , 2

σ(u) = λ (trε(u))I + 2μ ε(u),

(1.1) (1.2)

where u : Ω → R2 is the displacement vector, λ > 0 and μ > 0 are the Lame’s coefficients and I denotes 2 × 2 identity matrix. Let n be the unit outward normal vector and t be the tangential vector to ∂Ω. For a vector valued function v, its normal and tangential components on the boundary are defined by vn = v · n and vt = v − vn n, respectively . Similarly, for a matrix valued function q, its normal and tangential components are defined 2

by qn = (qn) · n and qt = qn − qn n, respectively. Define the space V by V : = {v ∈ [H 1 (Ω)]2 : v = 0 on ΓD , vn |ΓC = 0},

(1.3)

where H 1 (Ω) denotes the standard Sobolev space. Given f1 ∈ [L2 (Ω)]2 , f2 ∈ [L2 (ΓN )]2 and g0 ∈ L2 (ΓC ) with g0 > 0, the variational formulation for the contact problem with Tresca friction is to find u ∈ V such that a(u, v − u) + j(v) − j(u) ≥ l(v − u)

∀ v ∈ V,

where the bilinear form a(·, ·), the linear functional l(·) and the functional j(·) are defined by  σ(w) : ε(v) dx ∀ v, w ∈ V, a(w, v) = Ω  l(v) = f1 · v dx + f2 · v ds ∀ v ∈ V, Ω ΓN  g0 |vt | ds ∀ v ∈ V. j(v) =

(1.4)

(1.5) (1.6) (1.7)

ΓC

It is well known that the problem (1.4) admits a unique solution from the theory of elliptic variational inequalities [4, 30, 32, 44, 45, 47]. Define the space Λ as Λ = {μ ∈ [L∞ (ΓC )]2 : |μ| ≤ 1 a.e. on ΓC }. The following characterization of the problem (1.4) holds [9]: Lemma 1.1. There exists λt ∈ Λ such that a(u, v) + g(λt , v) = l(v) v ∈ V λt · ut = |ut | a.e. on ΓC

(1.8) (1.9)



where g(λt , v) =

ΓC

g0 λt · vt ds.

(1.10)

The model problem (1.4) satisfies the following strong formulation whenever the solution u is sufficiently regular [9]: σ(u) = λ (trε(u))I + 2μ ε(u) −divσ(u) = f1

in

Ω,

u=0

on

ΓD ,

σ(u)n = f2

on

ΓN , 3

in

Ω,

un = 0,

|σ(u)t | ≤ g0 ,

|σ(u)t | < g0 =⇒ ut = 0, |σ(u)t | = g0 =⇒ ut = −kσ(u)t for some k ≥ 0

on

ΓC .

Note that, the boundary condition on ΓC represent Tresca’s friction law. Under the assumption that u is sufficiently smooth, using a integration by parts in the equation (1.8) together with (1.10) and strong formulation, it can be easily obtained that σ(u)t = −g0 λt . The rest of the article is organized as follows: In next section, we present some useful preliminary results and introduce an enriching map with its important properties. In Section 3, we present discontinuous Galekin methods for the continuous problem (1.4), followed by that in Section 4, we discuss examples of several DG methods for (1.4). A posteriori analysis of DG methods for the frictional contact problem (1.4) has been established in Section 5. Further, in Section 6, we derive a priori bounds assuming the minimal regularity on the exact solution u. Numerical experiments are reported in the Section 7 to illustrate theoretical results. Finally, we present the conclusions of this article in Section 8. 2. Preliminaries 2.1. Notations We introduce the following notations, which will be used throughout this article: Th := a regular simplicial triangulations of Ω,

T := a triangle in Th , h := max{hT : T ∈ Th },

hT := diameter of T, Vhi := set of all vertices in Th that are in Ω, Vhb := set of all vertices in Th that are on ∂Ω, ¯D, VhD := set of all vertices in Th that are on Γ VhN := set of all vertices in Th that are on ΓN , ¯C , VhC := set of all vertices in Th that are on Γ VT := set of three vertices of T, Ve := set of two vertices of edge e, Ehi := set of all interior edges of Th , Ehb := set of all boundary edges on Γ EhD := set of all edges on ΓD ,

Eh := Ehi ∪ EhD

4

Vh := V i ∪ V b ,

EhN := {e ∈ Ehb : e ⊂ ΓN },

EhC := {e ∈ Ehb : e ⊂ ΓC },

Tp := set of all triangles sharing the vertex p, Ep := set of all edges connected to the vertex p, he := length of an edge e. To define the jump and the mean of discontinuous functions across inter element boundaries, we define the broken Sobolev space as H 1 (Ω, Th ) = {v ∈ L2 (Ω) : vT = v|T ∈ H 1 (T )

∀ T ∈ Th }.

For any d = (di ) ∈ Rn and z = (zi ) ∈ Rn , define the product d ⊗ z ∈ Rn×n as (d ⊗ z)ij = di zj . Let e ∈ Ehi be the interior edge shared by two neighboring triangles T+ and T− i.e. e = T¯+ ∩ T¯− . Let n+ be the unit normal of e pointing from T+ to T− , and n− = −n+ (cf. Fig. 2.1). For a scalar valued function v ∈ H 1 (Ω, Th ), a vector valued function w ∈ [H 1 (Ω, Th )]2 and a matrix valued function q ∈ [H 1 (Ω, Th )]2×2 , jumps [[·]] and averages {{·}} across the edge e are defined as follows [[v]] = v− n− + v+ n+ , and {{v}} =

1 (v− + v+ ), 2

[[w]] = w− ⊗ n− + w+ ⊗ n+ , and {{w}} =

[[q]] = q− n− + q+ n+ , and {{q}} =

1 (w− + w+ ), 2

1 (q + q+ ), 2 −

where  v ± = v T , ± w± = w|T± , q± = q|T± . We also define jump and mean on the boundary Γ. For any edge e ∈ Ehb , let T ∈ Th be a triangle such that e = ∂T ∩ ∂Ω. Let ne be the unit normal of e that points outside T . For any v ∈ H 1 (T ), we set on e ∈ Ehb [[v]] = vne and {{v}} = v, 5

P+

B @ @

@ @

@

T+

@

@

@ @

e

@  te @ . @ T− @ @ R ne @

@

@ A

@ P−

Figure 2.1: The edge e = ∂T− ∩ ∂T+ is shared by two neighboring triangles T− and T+ with initial node A, end node B and unit normal ne . The orientation of ne = n+ = −n− equals the outer normal of T+ , and hence, points into T− .

for w ∈ [H 1 (T )]2 , we define [[w]] = w ⊗ ne , and {{w}} = w, and for a matrix value function q ∈ [H 1 (T )]2×2 , we set [[q]] = qne , and {{q}} = q. The discontinuous finite element space Vh is defined by ˜ h : = {vh ∈ [L2 (Ω)]2 : vh |T ∈ [P1 (T )]2 V

∀ T ∈ Th },

˜ h : (v|e · n)(p) = 0 ∀e ∈ E C ∀p ∈ Ve }, Vh : = {v ∈ V h where P1 (T ) denotes the polynomial space of degree less than or equal to one defined on T . We will also require a conforming finite element space Vc := Vh ∩ V, which is chosen to be the Lagrange linear finite element space. Note that, normal and tangent vectors to edges are taken to be limits of normal and tangent vectors from the interior of e ⊂ ΓC if p ∈ e ⊂ ΓC is a corner of Γ. Moreover, in the definition of Vh , we do not put the constraint on the normal component of discrete functions on the elements T ∈ Th if T¯ ∩ ΓC is just a single point. Henceforth, C denotes a positive generic constant which takes different values at different appearances. We recall the following Clement type approximation result [17, Section 4.8] required in the further analysis in the article. Lemma 2.1. Let v ∈ V. Then there exists vh ∈ Vc and a positive constant C independent of h such that on any T ∈ Th ,

v − vh H l (T ) ≤ Ch1−l T v H 1 (TT ) , Here TT = {T  ∈ Th : T¯ ∩ T¯ = ∅}. 6

l = 0, 1.

(2.1)

We also need the following trace and inverse inequalities [17, Section 1.6, Section 4.5], [26, Theorem 3.2.6]. Lemma 2.2. Let v ∈ [H 1 (T )]2 for T ∈ Th and e be an edge of T . Then, it holds that

v L2 (e) ≤ C(h−1/2

v L2 (T ) + h1/2 e e ∇v L2 (T ) ).

(2.2)

Lemma 2.3. There exists a positive C independent of h such that for vh ∈ Vh ,

vh L∞ (e) ≤ Ch−1/2

vh L2 (e) , e

vh L2 (e) ≤

∇vh L2 (T ) ≤

Ch−1/2

vh L2 (T ) e −1 ChT vh L2 (T )

(2.3) ∀ T ∈ Th , ∀ T ∈ Th ,

(2.4) (2.5)

where e is an edge of T . 2.2. Enriching map In the analysis of discontinuous Galerkin methods, enriching functions which map nonconforming function to a conforming function plays an important role [13, 14, 15, 16, 43, 34]. Here, we construct an enriching function Eh : Vh → Vc by averaging as follows: If p ∈ Vhi ∪ VhN ∪ VhD , define

⎧ ⎪ ⎪ ⎪ ⎪ ⎨ Eh vh (p) =

1 |Tp |

⎪ ⎪ ⎪ ⎪ ⎩ 0

T ∈Tp

vh |T (p)

for p ∈ Vhi ∪ VhN

for p ∈ VhD ,

where |Tp | denotes the cardinality of Tp . If p ∈ VhC , define

(Eh vh · n)(p) = 0, (Eh vh · t)(p) =

1  (vh |e · t)(p), |EpC | C e∈Ep

where EpC = Ep ∩ EhC and ν(resp. τ ) is the unit normal (resp. tangent) vector to e ⊂ ΓC at p. Clearly Eh vh ∈ Vc for every vh ∈ Vh . We will require the following approximation properties of the smoothing map Eh : Lemma 2.4. It hold that  T ∈Th

and

 T ∈Th

2 h−2 T Eh v − v L2 (T ) ≤ C

 1

[[v]] 2L2 (e) ∀ v ∈ Vh , he

∇(Eh v − v) 2L2 (T ) ≤ C

 1

[[v]] 2L2 (e) ∀ v ∈ Vh . he

e∈Eh

e∈Eh

7

Proof. Let T ∈ Th be arbitrary. It follows from the scaling that  2 (Eh v − v|T )2 (p). h−2 T Eh v − v L2 (T ) ≤ C p∈VT

For p ∈ Vhi ∪ VhN , using the definition of Eh , we find (Eh v − v|T )(p) =

 1   v|T  (p) − v|T (p) . |Tp | 

(2.6)

T ∈Tp

Further, there are Ti ∈ Tp for i = 1, 2..., n, n ≥ 1, such that T1 = T and Tn = T  . Therefore, v|T  (p) − v|T (p) = v|T  (p) − v|Tn−1 (p) + v|Tn−1 (p) − v|Tn−2 (p) + v|Tn−2 (p) + ... + v|T3 (p) − v|T2 (p) + v|T2 (p) − v|T (p). By using the inverse inequality (2.3), we have |v|T  (p) − v|T (p)|2 ≤ C





[[v]] 2L∞ (e) ≤ C

i e∈Ep ∩Eh



i e∈Ep ∩Eh

therefore 2

|(Eh v − v|T )(p)| ≤ C



 e

i e∈Ep ∩Eh

e

1 [[v]]2 ds, he

1 [[v]]2 ds. he

For p ∈ VhD , we get (Eh v − v|T )(p) = −(v|T )(p), Using the same arguments as above and inverse inequality (2.3), we find   1 |(Eh v − v|T )(p)|2 ≤ C [[v]]2 ds. h e e D e∈Ep ∩Eh

For p ∈ VhC , we have  2  2 |Eh v − v|T |2 (p) = (Eh v − v) · n (p) + (Eh v − v) · t (p). Then, it follows from the definition of Eh and inverse inequality (2.3), |(Eh v − v|T )(p)|2 ≤ C

 i e∈Ep ∩Eh

 e

1 [[v]]2 ds. he

This proves the first inequality of this lemma. The second inequality of the Lemma 2.4, follows using the first inequality together with the inverse inequality (2.5).

8

3. Discrete problem For any v ∈ Vh , define ∇h (v), divh (v), εh (v) and σh (v) as follows: on any T ∈ Th ∇h (v) = ∇v, divh (v) = div(v), εh (v) = ε(v), σh (v) = 2μ εh (v) + λ tr(εh (v))I. Discontinuous Galerkin methods for (1.4) is to find uh ∈ Vh such that Ah (uh , vh − uh ) + j(vh ) − j(uh ) ≥ l(vh − uh )

∀ v h ∈ Vh ,

(3.1)

where the bilinear form Ah (v, w) is given by Ah (v, w) = ah (v, w) + bh (v, w)

(3.2)

with ah (v, w) =

 

2μ ε(v) : ε(w) + λ tr(ε(v))tr(ε(w)) dx,

T ∈Th



= Ω

T

2μ εh (v) : εh (w) + λ tr(εh (v))tr(εh (w)) dx,

(3.3)

and the bilinear form bh (v, w) consists of all the consistency and stability terms. Here, we do not specify a particular form for bh , rather, we assume that the bilinear from bh satisfies the following estimate: |bh (vh , wh )| ≤ C

 1

1/2 [[vh ]]2 ds

∇wh L2 (Ω) e he

∀ vh ∈ Vh and wh ∈ Vc .

(3.4)

e∈Eh

We refer to Section 4 for DG methods which satisfy the property (3.4). We further assume that the discrete problem (3.1) admits a unique solution. Before introducing a discrete Lagrange multiplier, we define Zh := {vh ∈ [L2 (ΓC )]2 : vh |e ∈ [P1 (e)]2 ∀e ∈ EhC }. Let γh : Vh → Zh be the trace map defined by γh (vh ) = vh |Γc . As in the case of the continuous problem, we have the following characterization of the discrete problem: Lemma 3.1. There exists a unique Lagrange multiplier λht ∈ Λ such that the solution uh of (3.1) can be characterized by the following equations Ah (uh , vh ) + g(λht , vht ) = l(vh ) λht · uht = |uht | 9

∀ v h ∈ Vh ,

(3.5)

a.e on

(3.6)

ΓC .

Proof. First we prove (3.5) and (3.6) using (3.1). By taking vh = 0 and vh = 2uh in (3.1), we find that Ah (uh , uh ) + j(uh ) = l(uh ).

(3.7)

Ah (uh , vh ) + j(vh ) ≥ l(vh ).

(3.8)

A use of (3.1) and (3.7) yields, Replacing vh by −vh in last equation, we get Ah (uh , vh ) − j(vh ) ≤ l(vh ).

(3.9)

Combining (3.8) and (3.9) , we obtain |Ah (uh , vh ) − l(vh )| ≤ j(vh ).

(3.10)

0 0 ⊥ ⊕ (Vht ) , where Let Wht be the restriction of tangential trace operator to Vh . We have that Vh = Vht 0 0 ⊥ 0 Vht = ker(Wht ) ∩ Vh and (Vht ) denotes the orthogonal compliment of Vht in Vh . Since Wht restricted to 0 ⊥ 0 ⊥ ) is an isomorphism, by setting Hht = Wht ((Vht ) ),the functional L : Hht → R defined by (Vht

˜h ) − l(˜ vh ) L(vh ) = Ah (uh , v vh ) = vh . Further, using is a continuous linear functional, where v ˜h ∈ Vh is a function such that Wht (˜ equation (3.10), we have that  |L(vh )| ≤ g0 |vht | ds. ΓC

2

Since Hht ⊂ [L1 (ΓC )] , by Hahn Banach theorem this mapping can be extended to [L1 (ΓC )]2 . Then, by the Riesz representation theorem we ensures the existence of λht ∈ [L∞ (ΓC )]2 , |λht | ≤ 1 a.e. in ΓC such that (3.5) holds. Further, using (3.7) and the fact that |λht | ≤ 1, we get λht · uht = |uht |. Conversely, we assume (3.5)-(3.6) hold and prove (3.1). To this end, using (3.5) and (3.6), we get   g0 λht · vht ds − g0 |uht | ds = (f, vh − uh ). Ah (uh , vh − uh ) + ΓC

Thus the inequality (3.1) follows, since j(vh ) ≥

 ΓC

ΓC

g0 λht · vht ds.

4. Examples of DG methods In this section, we study DG methods for which requirements of Section 3 are fulfilled and we verify that (3.4) holds. In [36] several DG methods have been discussed for a frictionless contact problem, for the sake of completeness, we are again presenting them here. Let r0 and re denote the global and local lifting operators, respectively [3, 55]. The bilinear form ah (., .) will be the same as in (3.3) for all DG methods. Here, we present the choice of bh for various DG methods:

1. SIPG method [2, 59, 42, 55]:  bh (v, w) = −

Eh

[[v]] : (2μ{{εh (w)}} + λtr({{εh (w)}}I)) ds 10

(4.1)

 −

 Eh

[[w]] : (2μ{{εh (v)}} + λtr({{εh (v)}}I)) ds +

Eh

η [[v]] : [[w]] ds he

for v, w ∈ Vh and η ≥ η0 > 0.

2. NIPG method [50, 55]:  [[v]] : (2μ{{εh (w)}} + λtr({{εh (w)}}I)) ds bh (v, w) = Eh   − [[w]] : (2μ{{εh (v)}} + λtr({{εh (v)}}I)) ds + Eh

(4.2)

Eh

η [[v]] : [[w]] ds he

for v, w ∈ Vh and η > 0.

3. Brezzi et al. [20, 21, 55]:  [[v]] : (2μ{{εh (w)}} + λtr({{εh (w)}}I)) ds bh (v, w) = − E  h − [[w]] : (2μ{{εh (v)}} + λtr({{εh (v)}}I)) ds E  h + r0 ([[w]]) : (2μr0 ([[v]]) + λtr(r0 ([[v]])I) dx Ω  + η(2μ re ([[v]]) : re ([[w]]) + λtr(re ([[v]])tr(re ([[w]]) dx e∈Eh

(4.3)

Ω

for v, w ∈ Vh and η > 0.

4. Bassi et al. [6, 55]:  [[v]] : (2μ{{εh (w)}} + λtr({{εh (w)}}I)) ds bh (v, w) = − E  h − [[w]] : (2μ{{εh (v)}} + λtr({{εh (v)}}I)) ds Eh  + η(2μ re ([[v]]) : re ([[w]]) + λtr(re ([[v]])tr(re ([[w]]) dx e∈Eh

(4.4)

Ω

for v, w ∈ Vh and η > 3.

5. LDG method [27, 24, 23, 55]:  [[v]] : (2μ{{εh (w)}} + λtr({{εh (w)}}I)) ds bh (v, w) = − E  h − [[w]] : (2μ{{εh (v)}} + λtr({{εh (v)}}I)) ds Eh

11

(4.5)





+ Ω

r0 ([[w]]) : (2μr0 ([[v]]) + λtr(r0 ([[v]])I) dx +

Eh

η [[v]] : [[w]] ds he

for v, w ∈ Vh and η > 0.

The norm . h on Vh is defined by

v 2h :=

 T ∈Th

ε(v) 2L2 (T ) +

 1

[[v]] 2L2 (e) . he

e∈Eh

Notice that, by using Korn’s type inequality [16], for any v ∈ Vh we get,  T ∈Th

∇v 2L2 (T ) ≤

 T ∈Th

ε(v) 2L2 (T ) +

 1

[[v]] 2L2 (e) . h e i

e∈Eh

The bilinear form Ah (., .) is bounded and coercive for all the DG methods discussed in this section [55] Ah (v, w) ≤ C2 v h w h Ah (v, v) ≥ C1 v 2h

∀ v, w ∈ Vh ,

∀ v ∈ Vh ,

(4.6) (4.7)

for some positive constants C1 and C2 . Hence, the discrete problem (3.1) admits a unique solution [32, 45] and the bound (3.4) on bh (., .) follows from its definition. 5. A posteriori analysis In this section, we study a posteriori error analysis of the error u − uh h . We introduce the following error estimators: η12 :=

 T ∈Th

η22 :=



h2T f1 2L2 (T ) ,

he [[σh (uh )]] 2L2 (e) ,

i e∈Eh

η32 :=



e∈Eh

η42 :=



2 h−1 e [[uh ]] L2 (e) ,

he σht (uh ) + g0 λht 2L2 (e) ,

C e∈Eh

η52 :=



he σh (uh )n − f2 2L2 (e) .

N e∈Eh

The total residual based estimator ηh is defined by ηh2 = η12 + η22 + η32 + η42 + η52 . 12

(5.1)

5.1. Reliability of a posteriori error estimator The following theorem establishes a relation between the error and the error estimator. Theorem 5.1. Let u and uh be the solution of (1.4) and (3.1), respectively. There exists positive constant C, depending only on the shape regularity of Th , μ and λ such that

u − uh 2h ≤ Cηh . Proof. We have,

u − uh h ≤ u − Eh uh h + Eh uh − uh h , where Eh is the enriching map defined in Section 2.2. Let φ = u − Eh uh and φh ∈ Vc be an approximation of φ as in Lemma 2.1. Denoting |v|2h = ah (v, v) and using the coercivity of the bilinear form a(·, ·), Lemma 1.1, equation (3.2) and Lemma 3.1, we get C|u − Eh uh |2h ≤ a(u − Eh uh , φ) = l(φ) − g(λt , φ) − a(Eh uh , φ) = l(φ − φh ) + l(φh ) − g(λt , φ) + ah (uh − Eh uh , φ) − ah (uh , φ) = l(φ − φh ) − ah (uh , φ − φh ) + l(φh ) − ah (uh , φh ) + ah (uh − Eh uh , φ) − g(λt , φ) = l(φ − φh ) − ah (uh , φ − φh ) + bh (uh , φh ) + l(φh ) − Ah (uh , φh ) + ah (uh − Eh uh , φ) − g(λt , φ) = l(φ − φh ) − ah (uh , φ − φh ) + bh (uh , φh ) + g(λht , φh ) + ah (uh − Eh uh , φ) − g(λt , φ) = l(φ − φh ) − ah (uh , φ − φh ) + bh (uh , φh ) − g(λht , φ − φh ) + ah (uh − Eh uh , φ) + g(λht , φ) − g(λt , φ). Now we evaluate the terms of the right hand side in the following manner. First, by integrating by parts, we find   ah (uh , φ − φh ) = σh (uh ) : εh (φ − φh ) dx = =

T ∈Th

T

T ∈Th

∂T

 



0 e∈Eh

+

e

σh (uh )n · (φ − φh ) ds

[[(φ − φh )]] : {σh (uh )}} ds +



N ∪E C e∈Eh h

 i e∈Eh

 e

e

{{(φ − φh )}} · [[σh (uh )]] ds

(φ − φh ) · σh (uh )n ds.

Notice that, on any e ⊂ ΓC , φn = 0 = φhn . Therefore, l(φ − φh ) − ah (uh , φ − φh )+bh (uh , φh ) − g(λht , φ − φh )    = f1 · (φ − φh ) dx − {{(φ − φh )}} · [[σh (uh )]] ds T ∈Th

T

i e∈Eh

13

e



+ bh (uh , φh ) −

0 e∈Eh

 

+

N e∈Eh

 

+

C e∈Eh

e

e

e

[[(φ − φh )]] : {σh (uh )}} ds

(φ − φh ) · (f2 − σh (uh )n) ds

(φt − φht ) · (−g0 λht − σht (uh ) ds

≤ ηh φ h , since, bh (uh , φh ) +

 0 e∈Eh

e

[[φ]] : {σh (uh )}} ds ≤

 1

1/2 [[uh ]]2 ds |φ|h . e he e∈Eh

Using bounded-ness of the bilinear form Ah (·, ·) and Lemma 2.4 , we find ah (uh − Eh uh , φ) ≤ uh − Eh uh h |φ|h ≤

2 1/2  1   |φ|h . [[uh ]] he 0,e e∈Eh

By the characterization (1.1), (3.6) and |λht | ≤ 1 , we obtain    g0 λht · ut ds − g0 λht · (Eh uh )t ds − g(λht , φ) − g(λt , φ) = ΓC

ΓC

 ΓC

g0 |ut | ds +

ΓC

g0 λt · (Eh uh )t ds

≤ g(λt , Eh uh ) − g(λht , Eh uh )  ≤ g0 |(Eh uh )t | ds − g(λht , Eh uh ) ΓC   g0 |(Eh uh )t − uht | ds + g0 |uht | ds − g(λht , Eh uh ) ≤ ΓC ΓC   g0 |(Eh uh )t − uht | ds + g0 λht · uht ds − g(λht , Eh uh ) ≤ ΓC ΓC  g0 |(Eh uh )t − uht | ds + g(λht , uh − Eh uh ) ≤ ΓC

≤2



g0 0,e Eh uh − uh 0,e .

C e∈Eh

Thus proof follows by using Lemma 2.4. 5.2. Efficiency of a posteriori error estimator We will use the following lemma to obtain efficiency estimates: Lemma 5.2. Let vh ∈ Vh . It hold that,  h2T f1 2L2 (T ) ≤ C( u − vh 2h + Osc(f1 )2 ), 

(5.2)

T ∈Th

he [[σh (vh )]] 2L2 (e) ≤ C( u − vh 2h + Osc(f1 )2 ),

i e∈Eh

14

(5.3)



he σht (vh ) + g0 λt 2L2 (e) ≤ C( u − vh 2h + Osc(f1 )2 + Osc1 (λt )2 + Osc1 (g0 )2 ),

(5.4)

C e∈Eh



he σh (vh )n − f2 2L2 (e) ≤ C( u − vh 2h + Osc(f1 )2 + Osc1 (f2 )2 ),

(5.5)

N e∈Eh

where, Osc(f1 )2 :=



h2T f1 − f1 20,T ,

T ∈Th 2

Osc1 (f2 ) :=



he f2 − f2 20,e ,

N e∈Eh

Osc1 (g0 )2 :=



he g0 − g 0 20,e ,

C e∈Eh

Osc1 (λt )2 :=



he λt − λt 20,e ,

C e∈Eh

with v denoting the L2 projection of any function v to the space of piece-wise constant functions. Proof. The estimates (5.2), (5.3) and (5.5) can be proved using the same argument as in [36]. Proof of (5.4) : Let e ⊂ ΓC and T ∈ Th be a triangle such that e ⊂ ∂T . Let be ∈ P2 (T ) denotes the bubble function which takes the unit value at the midpoint of edge e and vanishes on ∂T \e. Define φe = ξbe on T where ξ ∈ P0 (T ) is such that ξ = σht (vh ) + g 0 λt on e. Let φ be the extension of φe by zero to Ω. A use of integration by parts and Lemma 1.1 yields,   C σht (vh ) + g0 λt 2L2 (e) ≤ (σht (vh ) + g0 λt ) · φ ds + (g0 λt − g0 λt ) · φ ds e e ≤ ah (vh , φ) − a(u, φ) + l(φ) + ((g0 − g0 )λt + g0 (λt − λt )) · φ ds e

−1/2 ≤ (h−1

(g0 − g0 )λt 0,e T ∇(u − vh ) 0,T + f1 0,T + he

+ h−1/2

g0 (λt − λt ) 0,e ) φ 0,T . e Therefore, Ch1/2 e σht (vh ) + g0 λt L2 (e) ≤ ( ∇(u − vh ) 0,T + hT f1 0,T 1/2 + h1/2 e (g0 − g0 )λt 0,e + he g0 (λt − λt ) 0,e ) φ 0,T .

Now (5.4) follows by using (5.2) and inverse inequality (2.4). We now deduce the following efficiency estimates from Lemma 5.2: Theorem 5.3. Let u and uh be the solution of (1.4) and (3.1), respectively. Then, it holds that ηh2 =

5  i=1

ηi2 ≤ C( u − uh 2h +



he λt − λht 2L2 (e) + Osc(f1 )2 + Osc1 (f2 )2 + Osc1 (λt )2 + Osc1 (g0 )2 ).

C e∈Eh

15

Proof. Bound on η1 , η2 and η5 follow by taking vh = uh in (5.2), (5.3) and (5.5), respectively. From the definition of η3 , it is easy to observe that η3 ≤ u − uh h . To get the bound on η4 , observe that

σht (uh ) + g0 λht L2 (e) ≤ σht (uh ) + g0 λt L2 (e) + g0 (λht − λt ) L2 (e) . Taking vh = uh in equation (5.4) yields,   he σht (uh ) + g0 λt 2L2 (e) + he λt − λht 2L2 (e) ) η42 ≤ C( C e∈Eh

≤ C( u − uh 2h +



C e∈Eh

he λt − λht 2L2 (e) + Osc(f1 )2 + Osc1 (f2 )2 + Osc1 (λt )2 + Osc1 (g0 )2 ).

C e∈Eh

6. Medius analysis In this section, we establish the convergence analysis with the minimal regularity assumption on the exact solution u of (1.4) i.e. u ∈ [H 1+s (Ω)]2 with s ≥ 0 using the techniques from a priori and a posteriori analysis [33]. For any vh ∈ Vh , define η(vh ) as η(vh )2 =

 T ∈Th

+



h2T f1 + div(σh (vh )) 2L2 (T ) +



he [[σh (vh )]] 2L2 (e) +

i e∈Eh

he σht (vh ) + g0 λt 2L2 (e) +

C e∈Eh



 e∈Eh

2 h−1 e [[vh ]] L2 (e)

he σh (vh )n − f2 2L2 (e) .

N e∈Eh

Theorem 6.1. Let u and uh be the solution of (1.4) and (3.1), respectively. Then, 1

1

2 2

u − uh h ≤ C inf ( u − vh h + g0 0,Γ

u − vh 0,Γ ) + Osc(f1 ) + Osc1 (f2 ) + Osc1 (g0 ) + Osc1 (λt ). C C

vh ∈Vh

To prove this theorem, we will require the following lemma: Lemma 6.2. Let vh ∈ Vh be arbitrary and φ = ψ − Eh ψ for some ψ ∈ Vh . Then it holds that Ah (vh , φ) − l(φ) + g(λt , φ) ≤ η(vh ) ψ h . Proof. From (3.2), we have Ah (vh , φ) − l(φ) + g(λt , φ) = ah (vh , φ) + bh (vh , φ) − l(φ) + g(λt , φ) Integration by parts yields,   σh (vh ) : εh (φ) dx ah (vh , φ) = T ∈Th

=−

T

 

T ∈Th

T

div(σh (vh )) · φ dx +

  T ∈Th

16

∂T

σh (vh )n · φ ds

(6.1)

=−

  T

T ∈Th

+

div(σh (vh )) · φ dx +

0 e∈Eh



 N ∪E C e∈Eh h



e

e

[[φ]] : {σh (vh )}} ds +

 i e∈Eh

e

{{φ}} · [σh (vh )]] ds

φ · σh (vh )n ds.

Since on Γc , φn = 0, we get σh (vh )n · φ = [(σht (vh ) + σhn (vh ))n] · φ = σht (vh ) · φ. Therefore, Ah (vh , φ) − l(φ) + g(λt , φ) = −

  T ∈Th



+

0 e∈Eh

e

N e∈Eh

(f1 + div(σh (vh ))) · φ dx + bh (vh , φ)

[[φ]] : {σh (vh )}} ds +

 

+

T

e

 i e∈Eh

φ · (σh (vh ) − f2 )n ds +

e

{{φ}} · [σh (vh )]] ds

  C e∈Eh

e

φ · (σht (vh ) + g0 λt ) ds.

Now the desired estimate can be obtain by using the same argument as in Lemma 5.1 together with the observation that for all the DG methods introduced in Section 4, the following estimate holds  1

1/2  bh (vh , φ) + [[φ]] : {σh (vh )}} ≤ [[vh ]]2 ds

φ h . e e he 0 e∈Eh

e∈Eh

Now, we are ready to provide the proof of the Theorem 6.1. Proof. Let vh ∈ Vh be arbitrary. Setting ψ = vh − uh and using the triangle inequality to get,

u − uh h ≤ u − vh h + vh − uh h . Note that, using Lemma 1.1, we have a(u, Eh ψ) + g(λt , Eh ψ) = l(Eh ψ) By coercivity of Ah (·, ·) and equations (3.1), (6.2), we find C1 vh − uh 2h ≤ Ah (vh − uh , vh − uh ) = Ah (vh , ψ) − Ah (uh , vh − uh ) ≤ Ah (vh , ψ) − l(ψ) + j(vh ) − j(uh ) = Ah (vh , ψ − Eh ψ) − l(ψ − Eh ψ) − l(Eh ψ) + Ah (vh , Eh ψ) + j(vh ) − j(uh ) = Ah (vh , ψ − Eh ψ) − l(ψ − Eh ψ) + g(λt , ψ − Eh ψ) + Ah (vh , Eh ψ) − a(u, Eh ψ) − g(λt , ψ) + j(vh ) − j(uh ).

17

(6.2)

By using Lemma 6.2 and Lemma 5.2, we obtain Ah (vh , ψ − Eh ψ) − l(ψ − Eh ψ) + g(λt , ψ − Eh ψ) ≤ η(vh ) ψ h

≤ u − vh h + Osc(f1 ) + Osc1 (f2 ) + Osc1 (g0 ) + Osc1 (λt ) ψ h . Using (3.2), definitions of a(·, ·) , ah (·, ·), estimate (3.4) and Lemma 2.4, we get Ah (vh , Eh ψ) − a(u, Eh ψ) = ah (vh , Eh ψ) + bh (vh , Eh ψ) − a(u, Eh ψ)

1/2  1   [[vh − u]]2 ds

∇Eh ψ L2 (Ω) (σh (vh ) − σ(u)) : ε(Eh ψ) dx + ≤ T e he T ∈Th

e∈Eh

≤ u − vh h Eh ψ h ≤ C u − vh h ψ h . Now, a use of |λt | ≤ 1 and Lemma 3.1 yields,     j(vh ) − j(uh ) − g(λt , ψ) = g0 |vht | ds − g0 |uht | ds − g0 λt · vht ds + g0 λt · uht ds ΓC ΓC ΓC ΓC   g0 |vht | ds − g0 λt vht ds ≤ ΓC ΓC   g0 (|vht | − |ut |) ds − g0 λt · (vht − ut ) ds = ΓC ΓC  g0 |vht − ut | ds ≤ 2 g0 0,ΓC u − vh 0,ΓC . ≤2 ΓC

Thus,

vh − uh 2h ≤ C( u − vh h +Osc(f1 )+Osc1 (f2 )+Osc1 (g0 )+Osc1 (λt )) vh − uh h +2 g0 0,ΓC u−vh 0,ΓC . Hence



u − uh 2h ≤ C u − vh 2h + ( u − vh h + Osc(f1 ) + Osc1 (f2 ) + Osc1 (g0 ) + Osc1 (λt )) vh − uh h

+ g0 0,ΓC u − vh 0,ΓC .

Finally, a use of triangle inequality and Young’s inequality completes the proof. 3 s

1

Remark 6.3. If u ∈ [H (1+s) (Ω)]2 , s ≥ 0, then we obtain that u − uh h converges with order hmin { 4 , 2 + 4 } due to the suboptimal convergence of the term inf u − vh 0,ΓC . If u ∈ [H 1 (Ω)]2 , then we can show that inf u − vh h → 0 and

vh ∈Vh

vh ∈Vh

inf u − vh 0,ΓC → 0 as h → 0 using the density argument, see for example

vh ∈Vh

[26, Theorem 3.2.3]. Hence, u − uh h → 0 as h → 0. 7. Numerical results In this section, we present numerical examples to illustrate the theoretical results discussed in this article. We use the Uzawa algorithm [32, 38] to compute the discrete solution uh . The following algorithm has been used for the adaptive mesh refinement 18

SOLVE → ESTIMATE → MARK → REFINE In the step SOLVE, we solve the discrete problem to find the discrete solution uh . Next, in the step ESTIMATE, we compute the error estimator ηh element-wise and use it for marking the elements in the step MARK with D¨ orfler marking strategy [28] with parameter θ = 0.3. In the last step REFINE, we refine the marked elements by using the newest vertex bisection algorithm to obtain new mesh and repeat the algorithm. Here, we present three test examples. The first test example is constructed such that the exact solution u is known in the closed form. Therein, we can compute the exact error and compare it with a posteriori error estimator. The data of the Example 2 and Example 3 are motivated from [9]. Therein, the exact solution u is unknown. We have discussed the numerical results for the SIPG and NIPG finite element methods. The penalty parameter η is chosen to be 10μ for the SIPG method and 10 for the NIPG method. In the examples, E and ν denote the Young’s modulus and the Poisson ratio, respectively and then the Lam´ e parameters μ and λ are given by μ=

E 2(1 + ν)

and

λ=

Eν . (1 + ν)(1 − 2ν)

Example 1. We consider Ω = (0, 1)2 and g0 = 1.10. The top of the unit square is kept fixed and a Neumann force is acting on the left and right sides of the square, which can be computed by the exact solution u = (y(y − 1)ey , y(y − 1)). The bottom of the unit square is the contact boundary. In this example, we take μ = 1, λ = 1. The convergence history of SIPG and NIPG methods on the uniform mesh is shown in Table 7.1, figures 7.1 and 7.2, respectively. Figures 7.3 and 7.4 illustrate the behavior of the error estimator and the energy norm error for SIPG and NIPG methods, respectively on adaptive meshes. Comparing figures 7.1, 7.2 to √ figures 7.3, 7.4, we observe that the error and the estimator converges optimally (1/ N , N = degrees of freedom) on the adaptive mesh while convergence on the uniform mesh is sub-optimal. Also, from figures

1/2 7.1 and 7.2, we observe that the term e∈EC he λt − λht L2 (e) appearing in reliability estimate converges with higher rate than the error u − uh h . Figures 7.3 and 7.4 ensure the reliability of the error estimator. Example 2. Here, the domain Ω := (0, 4) × (0, 4). The Neumann boundary condition is applied on the left side and on the top of the square, namely {0} × (0, 4) ∪ ((0, 4) × {4}). The bottom of the square is in contact with a rigid foundation, which is the contact boundary ΓC , i.e. ΓC = ((0, 4) × {0}). The right side 19

h 1/2 1/4 1/8 1/16 1/32 1/64

∇h (u − uh )L2 (Ω) 1.5403 6.3494 3.7353 2.3296 1.5301 1.0386

order of conv.

h

– 1.2785 0.7654 0.6811 0.6065 0.5591

1/2 1/4 1/8 1/16 1/32 1/64

×10−0 ×10−1 ×10−1 ×10−1 ×10−1 ×10−1

∇h (u − uh )L2 (Ω) 1.4830 6.1448 3.6907 2.3153 1.5224 1.0332

×10−0 ×10−1 ×10−1 ×10−1 ×10−1 ×10−1

order of conv. – 1.2711 0.7355 0.6727 0.6049 0.5593

Table 7.1: Errors and orders of convergence for SIPG and NIPG methods on uniform mesh for Example 1

10

1

True Error Estimator Lambda Error 10

0

1

Errors and Estimator

0.29

10

-1

0.28 1

10

-2

0.49 1

10

-3

10 1

10 2

10 3

10 4

10 5

10 6

Degrees of Freedom

Figure 7.1: u − uh h , ηh and h1/2 λt − λht L2 (E C ) on uniform mesh for SIPG for Example 1 h

10

1

True Error Estimator Lambda Error 10

0

1

Errors and Estimator

0.29

10

-1

0.28 1

10

-2

0.49 1

10

-3

10 1

10 2

10 3

10 4

10 5

10 6

Degrees of Freedom

Figure 7.2: u − uh h , ηh and h1/2 λt − λht L2 (E C ) on uniform mesh for NIPG for Example 1 h

of the square is kept fixed, i.e. the zero Dirichlet boundary condition is applied on {4} × (0, 4). We consider the following data E = 2000, f1 = 0,

ν = 0.4,

g0 = 150

f2 = (500(5 − y), 500), 20

True Error Estimator Optimal Rate 0

Error and Estimator

10

−1

10

−2

10

2

10

3

4

10

5

10

10

Degrees of Freedom

Figure 7.3: Error and Estimator on adaptive mesh for SIPG for Example 1

True Error Estimator Optimal Rate 0

Error and Estimator

10

−1

10

−2

10

2

10

3

4

10

10

5

10

Degrees of Freedom

Figure 7.4: Error and Estimator on adaptive mesh for NIPG for Example 1

Let N denotes degrees of freedom at a refinement level. Table 7.2 shows the behavior of log(ηh ) with respect to log(N ) for SIPG and NIPG methods, demonstrating the convergence of order N −α , α ≈ 0.5 of ηh . The decay of the error estimator with the increasing number of degrees of freedom for SIPG and NIPG methods have been shown in Figures 7.5 and 7.6, respectively. We observe that the error estimator converges optimally for both the methods, exhibiting some pre-asymptotic regime. Figures 7.7 and 7.8 show the adaptive mesh refinement at certain level for the SIPG and NIPG methods, respectively. The higher mesh refinement can be observed near the intersection of the boundaries and near the contact edge, which is due to the deformation in the body under the action of traction. Thus the error estimator captures well the singular behavior of the discrete solution. Example 3. In this example, the domain Ω := (0, 8) × (0, 4) and we have the following data:

21

log(N ) 1.0095 1.0412 1.0734 1.1048 1.1363 1.1679 1.1999 1.2318

×10+1 ×10+1 ×10+1 ×10+1 ×10+1 ×10+1 ×10+1 ×10+1

log(ηh )

order of conv.

6.5665 6.4299 6.2620 6.1059 5.9513 5.8193 5.6580 5.4839

0.4309 0.5214 0.4971 0.4908 0.4177 0.5041 0.5458

log(N )

log(ηh )

order of conv.

6.6142 6.4822 6.3340 6.1651 6.0136 5.8725 5.7289 5.5657

0.4314 0.4750 0.5484 0.4903 0.4522 0.4588 0.5231

×10+1 ×10+1 ×10+1 ×10+1 ×10+1 ×10+1 ×10+1 ×10+1

1.0036 1.0342 1.0654 1.0962 1.1271 1.1583 1.1896 1.2208

Table 7.2: log(N ), log(ηh ) and orders of convergence (from refinement level 22 to level 30) for SIPG and NIPG methods for Example 2

4

10

Estimator Optimal Rate

3

Estimator

10

2

10

1

10

2

10

3

4

10

10

5

10

6

10

Degrees of Freedom

Figure 7.5: Estimator for SIPG for Example 2

4

10

Estimator Optimal Rate

3

Estimator

10

2

10

1

10

2

10

3

4

10

10

5

10

Degrees of Freedom

Figure 7.6: Estimator for NIPG for Example 2

ΓC = (0, 8) × {0}, ΓD = (0, 8) × {4}, ΓN = {0} × (0, 4) ∪ {8} × (0, 4)

22

6

10

4

3.5

3

2.5

2

1.5

1

0.5

0

0

0.5

1

1.5

2

2.5

3

3.5

4

Figure 7.7: Adaptive mesh for SIPG at level 25 for Example 2

4

3.5

3

2.5

2

1.5

1

0.5

0

0

0.5

1

1.5

2

2.5

3

3.5

4

Figure 7.8: Adaptive mesh for NIPG at level 25 for Example 2

E = 1000, ν = 0.3, g0 = 150 f1 = 0 and f2 = (1000, 0). The decay of log(ηh ) with respect to log(N ) for SIPG and NIPG methods is shown in Table 7.3. Figures 7.9 and 7.10 show convergence behavior of the error estimator ηh with the increasing number of degrees of freedom for SIPG and NIPG methods, respectively. From these figures and Table 7.3, it is clear that 23

the error estimator converges optimally, except both SIPG and NIPG methods have some pre-asymptotic regime. The adaptive mesh refinement at the level 25 for SIPG and NIPG methods have been shown in Figures 7.11 and 7.12. The local mesh refinement is the effect of traction and the change in the behavior of the solution (from sliding to sticking). log(N ) 1.0064 1.0404 1.0748 1.1099 1.1434 1.1782 1.2132 1.2489

×10+1 ×10+1 ×10+1 ×10+1 ×10+1 ×10+1 ×10+1 ×10+1

log(ηh )

order of conv.

6.3422 6.2225 6.0794 5.8797 5.6820 5.5435 5.3899 5.2010

0.3521 0.4160 0.5690 0.5902 0.3980 0.4389 0.5291

log(N ) 9.9440 1.0281 1.0628 1.0982 1.1324 1.1665 1.2018 1.2382

log(ηh )

order of conv.

6.4116 6.2667 6.1419 5.9628 5.7519 5.5742 5.4403 5.2741

0.4300 0.3597 0.5059 0.6167 0.5211 0.3794 0.4566

×10+0 ×10+1 ×10+1 ×10+1 ×10+1 ×10+1 ×10+1 ×10+1

Table 7.3: log(N ), log(ηh ) and orders of convergence (from refinement level 22 to level 29) for SIPG and NIPG methods for Example 3

4

10

Estimator Optimal Rate 3

Estimator

10

2

10

1

10

2

10

3

4

10

10

5

10

6

10

Degrees of Freedom

Figure 7.9: Estimator for SIPG for Example 3

8. Conclusions In this paper, we have proposed discontinuous Galerkin (DG) methods for a frictional contact problem in linear elasticity. Therein, we derived a residual based a posteriori error estimator for the proposed class of DG methods. The reliability and the efficiency of the error estimator is discussed. We have introduced a smoothing function mapping DG finite element space to conforming finite element space which plays a crucial role in the analysis. We have also investigated a priori convergence analysis with minimal regularity assumption on the exact solution. Numerical results illustrate the performance of the error estimator. Though results of this article are discussed in two dimensions, they can also be extended in three dimensions using similar arguments. The results of this article also hold for conforming finite element methods. 24

4

10

Estimator Optimal Rate

3

Estimator

10

2

10

1

10

2

3

10

4

10

5

10

10

Degrees of Freedom

Figure 7.10: Estimator for NIPG for Example 3

4

3.5

3

2.5

2

1.5

1

0.5

0

0

1

2

3

4

5

6

7

8

Figure 7.11: Adaptive mesh for SIPG at level 26 for Example 3

9. Acknowledgment The author would like to thank anonymous reviewers for their helpful and constructive comments that lead to the improvement of this article. References [1] M. Andres, Dual-dual formulations for frictional contact problems in mechanics, PhD thesis, Leibniz Universit¨ at Hannover, 2011.

25

4

3.5

3

2.5

2

1.5

1

0.5

0

0

1

2

3

4

5

6

7

8

Figure 7.12: Adaptive mesh for NIPG at level 24 for Example 3

[2] D.N. Arnold, An interior penalty finite element method with discontinuous elements, SIAM J. Numer. Anal. 19 (1982) 742–760. [3] D.N. Arnold, F. Brezzi, B. Cockburn, L.D. Marini, Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM J. Numer. Anal. 39 (2002) 1749–1779. [4] K. Atkinson, W. Han, Theoretical Numerical Analysis. A functional analysis framework, Third edition, Springer, 2009. [5] S. Bartels, C. Carstensen, Averaging techniques yield relaible a posteriori finite element error control for obstacle problems, Numer. Math. 99 (2004) 225–249. [6] F. Bassi, S. Rebay, G. Mariotti, S. Pedinotti, M. Savini, A high order accurate discontinuous finite element method for inviscid and viscous turbomachinery flows, Proceedings of 2nd European Conference on Turbomachinery, Fluid Dynamics and Thermodynamics, R. Decuypere and G. Dibelius, eds., Technologisch Instituut, Antwerpen, Belgium, 1997, pp. 99–178. [7] H. Blum, F.T. Suttmeier, An adaptive finite element discretization for a simplified Signorini problem, Calcolo 37 (2000) 65–77. [8] V. Bostan, W. Han, Recovery-based error estimation and adaptive solution of elliptic variational inequalities of the second kind, Commun. Math. Sci. 2 (2004) 1–18.

26

[9] V. Bostan, W. Han, A posteriori error analysis for finite element solutions of a frictional contact problem, Comput. Methods Appl. Mech. Engrg. 195 (2006) 1252–1274. [10] V. Bostan, W. Han, Adaptive finite element solution of variational inequalities with application in contact problem in: D.y. gao, h.d. sherali (eds.), Advances in Applied Mathematics and Global Optimization, Springer, New York, 2009, pp. 25–106. [11] V. Bostan, W. Han, B. Reddy, A posteriori error estimation and adaptive solution of elliptic variational inequalities of the second kind, Appl. Numer. Math. 52 (2004) 13–38. [12] D. Braess, A posteriori error estimators for obstacle problems-another look, Numer. Math. 101 (2005) 415–421. [13] S. Brenner, Two-level additive schwarz preconditioners for nonconforming finite element methods, Math. Comp. 65 (1996) 897–921. [14] S. Brenner, Convergence of nonconforming multigrid methods without full elliptic regularity, Math. Comp. 68 (1999) 25–53. [15] S. Brenner, Ponicar´e-friedrichs inequalities for piecewise H 1 functions, SIAM J. Numer. Anal 41 (2003) 306–324. [16] S. Brenner, Korn’s inequalities for piecewise H 1 vector fields, Math. Comp. 73 (2004) 1067–1087. [17] S. Brenner, L. Scott, The Mathematical Theory of Finite Element Methods, Third edition, SpringerVerlag, New York, 2008. [18] S. Brenner, L.Y. Sung, Y. Zhang, Finite element methods for the displacement obstacle problem of clamped plates, Math. Comp. 81 (2012) 1247–1262. [19] F. Brezzi, W.W. Hager, P.A. Raviart, Error estimates for the finite element solution of variational inequalities, Part I. Primal theory, Numer. Math. 28 (1977) 431–443. [20] F. Brezzi, G. Manzini, D. Marini, P. Pietra, A. Russo, Discontinuous finite elements for diffusion problems, in Atti Convegno in onore di F. Brioschi (Milan, 1997), Istituto Lombardo, Accademia di Scienze e Lettere, Milan, Italy, 1999, pp. 197–217. [21] F. Brezzi, G. Manzini, D. Marini, P. Pietra, A. Russo, Discontinuous Galerkin approximations for elliptic problems, Numer. Meth. PDE 16 (2000) 365–378.

27

[22] M. B¨ urg, A. Schr¨oder, A posteriori error control of hp-finite elements for variational inequalities of the first and second kind, Computers and Mathematics with Applications 70 (2015) 2783–2802. [23] R. Bustinza, F.J. Sayas, Error estimates for an LDG method applied to a Signorini type problems, J. Sci. Comput. 52 (2012) 322–339. [24] P. Castillo, B. Cockburn, I. Perugia, D. Sch¨otzau, An a priori error analysis of the local discontinuous Galerkin method for elliptic problems, SIAM J. Numer. Anal. 38 (2000) 1676–1706. [25] Z. Chen, R. Nochetto, Residual type a posteriori error estimates for elliptic obstacle problems, Numer. Math. 84 (2000) 527–548. [26] P. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978. [27] B. Cockburn, C.W. Shu, The local discontinuous Galerkin method for time-dependent convectiondiffusion systems, SIAM J. Numer. Anal. 35 (1998) 2440–2463. [28] W. D¨orlfer, A convergent adaptive algorithm for Poisson’s equation, SIAM J. Numer. Anal. 33 (1996) 1106–1124. [29] P. D¨orsek, J.M. Melenk, Adaptive hp-fem for the contact problem with Tresca friction in linear elasticity: The primal–dual formulation and a posteriori error estimation, Appl. Numer. Math. 60 (2010) 689–704. [30] G. Duvaut, J. Lions, Inequalities in Mechanics and Physics, Springer, Berlin, 1976. [31] R.S. Falk, Error estimates for the approximation of a class of variational inequalities, Math. Comp. 28 (1974) 963–971. [32] R. Glowinski, Numerical Methods for Nonlinear Variational Problems, Springer-Verlag, New York, 1984. [33] T. Gudi, A new error analysis for discontinuous finite element methods for linear elliptic problems, Math. Comp. 79 (2010) 2169–2189. [34] T. Gudi, K. Porwal, A posteriori error control of discontinuous Galerkin methods for elliptic obstacle problems, Math. Comp. 83 (2014) 579–602. [35] T. Gudi, K. Porwal, A remark on the a posteriori error analysis of discontinuous Galerkin methods for obstacle problem, Comput. Meth. Appl. Math. 14 (2014) 71–87. [36] T. Gudi, K. Porwal, An a posteriori error estimator for a class of discontinuous Galerkin methods for Signorini problem, J. Comp. Appl. Math. 292 (2016) 257–278. 28

[37] T. Gudi, K. Porwal, A C0 interior penalty method for a fourth-order variational inequality of the second kind, Numer. Meth. PDE 32 (2016) 36–59. [38] D. Hage, N. Klein, F. Suttmeier, Adaptive finite elements for a certain class of variational inequalities of the second kind, Calcolo 48 (2011) 293–305. [39] W. Han, L. Wang, Nonconforming finite element analysis for a plate contact problem, SIAM J. Numer. Anal. 40 (2002) 1683–1697. [40] P. Hild, S. Nicaise, A posteriori error estimates of residual type for Signorini’s problem, Numer. Math. 101 (2005) 523–549. [41] P. Hild, S. Nicaise, Residual a posteriori error estimators for contact problems in elasticity, ESAIM:M2AN 41 (2007) 897–923. [42] J. J. Douglas, T. Dupont, Interior penalty procedures for elliptic and parabolic Galerkin methods, Lecture Notes in Phys., 58 : Springer-Verlag, Berlin, 1976. [43] O. Karakashian, F. Pascal, A posteriori error estimates for a discontinuous Galerkin approximation of second-order elliptic problems, SIAM J. Numer. Anal. 41 (2003) 2374–2399. [44] N. Kikuchi, J.T. Oden, Contact problems in elasticity, SIAM, Philadelphia, 1988. [45] D. Kinderlehrer, G. Stampacchia, An Introduction to Variational Inequalities and Their Applications, SIAM, Philadelphia, 2000. [46] R. Nochetto, T.V. Petersdorff, C.S. Zhang, A posteriori error analysis for a class of integral equations and variational inequalities, Numer. Math. 116 (2010) 519–552. [47] P. Panagiotopoulos, Inequality Problems in Mechanics and Applications, Birkh¨auser, Boston, 1985. [48] A. Rademacher, A. Schr¨ oder, Dual weighted residual error control for frictional contact problems, Comput. Meth. Appl. Math. 15 (2015) 391–413. [49] W.H. Reed, T.R. Hill, Triangular mesh methods for the neutron transport equation, Technical Report LA-UR-73-479, Los Alamos Scientific Laboratory, 1973. [50] B. Rivi`ere, M. Wheeler, V. Girault, A priori error estimates for finite element methods based on discontinuous approximation spaces for elliptic problems, SIAM J. Numer. Anal. 39 (2001) 902–931.

29

[51] A. Schr¨oder, Mixed finite element methods of higher-order for model contact problems, SIAM J. Numer. Anal. 49 (2011) 2323–2339. [52] F.T. Suttmeier, General approach for a posteriori error estimates for finite element solutions of variational inequalities, Comput. Mech. 27 (2001) 317–323. [53] A. Veeser, Efficient and reliable a posteriori error estimators for elliptic obstacle problems, SIAM J. Numer. Anal. 39 (2001) 146–167. [54] F. Wang, W. Han, X.Cheng, Discontinuous Galerkin methods for solving elliptic variational inequalities, SIAM J. Numer. Anal. 48 (2010) 708–733. [55] F. Wang, W. Han, X.Cheng, Discontinuous Galerkin methods for solving the Signorini problem, IMA J. Numer. Anal. 41 (2011) 1754–1772. [56] F. Wang, W. Han, X.Cheng, Another view for aposteriori error estimates for variational inequalities of the second kind, Appl. Numer. Math. 78 (2013) 225–233. [57] L. Wang, On the quadratic finite element approximation to the obstacle problem, Numer. Math. 92 (2002) 771–778. [58] A. Weiss, B.I. Wohlmuth, A posteriori error estimator and error control for contact problems, Math. Comp. 78 (2009) 1237–1267. [59] M. Wheeler, An elliptic collocation-finite-element method with interior penalties, SIAM J. Numer. Anal. 15 (1978) 152–161.

30