Symmetric and non-symmetric variants of Nitsche’s method for

contact problems in elasticity: theory and numerical experiments

∗ † ‡
Franz Chouly Patrick Hild Yves Renard

September 5, 2013

Abstract
A general Nitsche method, which encompasses symmetric and non-symmetric variants, is
proposed for frictionless unilateral contact problems in elasticity. The optimal convergence of
the method is established both for two and three-dimensional problems and Lagrange affine
and quadratic finite element methods. Two and three-dimensional numerical experiments
illustrate the theory.

Key words: unilateral contact, finite elements, Nitsche’s method.

AMS Subject Classification: 65N12, 65N30, 74M15.

1 Introduction
In solid mechanics, the numerical implementation of contact and impact problems generally uses
the Finite Element Method (FEM) (see [28, 18, 21, 19, 31, 41, 39]). The non-linear contact
conditions involved on a part of the boundary lead naturally to a variational inequality (see, e.g.,
[15]).
We consider in this paper a special FEM inspired from Nitsche’s method [34] (see also [3] for an
early extension to the Discontinuous Galerkin framework) in which the coercivity condition for
the domain bilinear form Aθγ (see (9)) is the same as in the standard formulation (see, e.g., [4]).
This method aims at treating the boundary or interface conditions in a weak sense, thanks to a
consistent penalty term. So it differs from standard penalization techniques which are typically
non-consistent [28]. Moreover, unlike mixed methods (see, e.g., [21]), no additional unknown
(Lagrange multiplier) is needed. Nitsche’s method has been widely applied during these last years
to problems involving linear conditions on the boundary of a domain or in the interface between
sub-domains: see, e.g,. [37] for the Dirichlet problem or [4] for domain decomposition with non-
matching meshes. More recently, in [20, 22] it has been adapted for bilateral (persistent) contact,
which still involves linear boundary conditions on the contact zone (note that an algorithm for
unilateral contact which involves Nitsche’s method in its original form is given and implemented
in [20]). An extension to large strain bilateral contact has been performed in [42].
∗
Laboratoire de Mathématiques de Besançon - UMR CNRS 6623, Université de Franche Comté, 16 route de
Gray, 25030 Besançon Cedex, France. email: [Link]@[Link]
†
Institut de Mathématiques de Toulouse - UMR CNRS 5219, Université Paul Sabatier, 118 route de Narbonne,
31062 Toulouse Cedex 9, France. email: [Link]@[Link]
‡
Université de Lyon, CNRS, INSA-Lyon, ICJ UMR5208, LaMCoS UMR5259, F-69621, Villeurbanne, France.
email: [Link]@[Link]

1
In a previous work a symmetric Nitsche-based method has been proposed for (non-linear) uni-
lateral contact conditions (see [9]). We established its optimal convergence in the H 1 (Ω)-norm
1 3
of order O(h 2 +ν ) for a solution of regularity H 2 +ν (Ω) (0 < ν ≤ 21 ). Furthermore, no additional
assumption on the contact zone is needed for the proof, such as an increased regularity of the
contact stress or a finite number of transitions between contact and non-contact. Besides, the
standard FEM for contact consists in a direct conforming approximation of the variational in-
equality, with the elastic displacement as the only unknown. For this standard FEM and also
for all the other approaches such as mixed/hybrid methods (e.g., [23, 6, 30]), stabilized mixed
methods (e.g., [25]), penalty methods (e.g., [10]), no such proof of optimal convergence has been
3
established to the best of our knowledge in the case the solution u is in H 2 +ν (Ω) (0 < ν ≤ 21 ).
We refer e.g. to [27, 26] for recent reviews on a priori error estimates for contact problems.
In the present paper we introduce an extension of the symmetric method of [9], which allows
nonsymmetric variants, depending upon a new parameter called θ. The symmetric case is recov-
ered when θ = 1. In this new method the advantage of the symmetric formulation consisting in
the positivity of the Nitsche penalty term is generally lost. Nevertheless this extension presents
some other advantages, mostly from the numerical viewpoint. In particular, one of its variants
(θ = 0) involves a reduced quantity of terms, which makes it easier to implement and to extend
to contact problems involving non-linear elasticity. Also, for θ = −1, the well-posedness of the
discrete formulation and the optimal convergence are preserved irrespectively of the value of the
Nitsche parameter. For all the values of θ, the optimal convergence rates can still be obtained.
We prove this optimal convergence property and illustrate it on several numerical experiments.
Note that the behavior of the generalized Newton algorithm when applied to this method has
already been studied in [36], in the frictionless and frictional cases and for two values of θ (0 and
1). This study illustrated the better numerical performances of the nonsymmetric variant θ = 0,
which requires less Newton iterations to converge, for a wider range of the Nitsche parameter.
Values of θ different from 0 and 1 have not been studied previously.
Our paper is outlined as follows. In Section 2, we recall the continuous (strong and weak)
formulations for unilateral contact problems and introduce our Nitsche-based FEM. In Section
3, we carry out the numerical analysis of this method: we prove its consistency, the existence
and uniqueness of solutions, and at last its optimal convergence. Numerical experiments in 2D
and 3D are described in Section 4, which illustrate the convergence properties for different values
of the numerical parameters, and in particular θ. In Section 5 conclusions are drawn and some
perspectives are given.
Let us introduce some useful notations. In what follows, bold letters like u, v, indicate vector or
tensor valued quantities, while the capital ones (e.g., V, K . . .) represent functional sets involving
vector fields. As usual, we denote by (H s (.))d , s ∈ R, d = 1, 2, 3, the Sobolev spaces in one, two or
three space dimensions (see [1]). The usual norm of (H s (D))d is denoted by k · ks,D and we keep
the same notation when d = 1 or d > 1. The letter C stands for a generic constant, independent
of the discretization parameters.

2 Setting
2.1 The unilateral contact problem
We consider an elastic body whose reference configuration is represented by the domain Ω in
Rd with d = 2 or d = 3. Small strain assumptions are made, as well as plane strain when
d = 2. The boundary ∂Ω of Ω is polygonal or polyhedral and we suppose that ∂Ω consists

2
in three nonoverlapping parts ΓD , ΓN and the contact boundary ΓC , with meas(ΓD ) > 0 and
meas(ΓC ) > 0. The contact boundary is supposed to be a straight line segment when d = 2 or
a polygon when d = 3 to simplify. The unit outward normal vector on ∂Ω is denoted n. In its
initial stage, the body is in contact on ΓC with a rigid foundation (the extension to two elastic
bodies in contact can be easily made, at least for small strain models) and we suppose that the
unknown final contact zone after deformation will be included into ΓC . The body is clamped on
ΓD for the sake of simplicity. It is subjected to volume forces f ∈ (L2 (Ω))d and to surface loads
g ∈ (L2 (ΓN ))d .
The unilateral contact problem in linear elasticity consists in finding the displacement field u :
Ω → Rd verifying the equations and conditions (1)–(2):
div σ(u) + f = 0 in Ω,
σ(u) = A ε(u) in Ω,
(1)
u=0 on ΓD ,
σ(u)n = g on ΓN ,
where σ = (σij ), 1 ≤ i, j ≤ d, stands for the stress tensor field and div denotes the divergence
T
operator of tensor valued functions. The notation ε(v) = (∇v+∇v )/2 represents the linearized
strain tensor field and A is the fourth order symmetric elasticity tensor having the usual uniform
ellipticity and boundedness property. For any displacement field v and for any density of surface
forces σ(v)n defined on ∂Ω we adopt the following notation
v = vn n + vt and σ(v)n = σn (v)n + σt (v),
where vt (resp. σt (v)) are the tangential components of v (resp. σ(v)n). The conditions
describing unilateral contact without friction on ΓC are:

un ≤ 0, (i)
σn (u) ≤ 0, (ii)
(2)
σn (u) un = 0, (iii)
σt (u) = 0. (iv)
We introduce the Hilbert space V and the convex cone K of admissible displacements which
satisfy the noninterpenetration on the contact zone ΓC :
n
d o
V := v ∈ H 1 (Ω) : v = 0 on ΓD , K := {v ∈ V : vn = v · n ≤ 0 on ΓC } .

We define also
Z Z Z
a(u, v) := σ(u) : ε(v) dΩ, L(v) := f · v dΩ + g · v dΓ,
Ω Ω ΓN

for any u and v in V. From the previous assumptions, we deduce that a(·, ·) is bilinear, symmetric,
V-elliptic and continuous on V × V. We see also that L(·) is a continuous linear form on V. The
weak formulation of Problem (1)-(2), as a variational inequality (see [15, 21, 28]), reads as:

Find u ∈ K such that:
(3)
a(u, v − u) ≥ L(v − u), ∀ v ∈ K.
Stampacchia’s Theorem ensures that Problem (3) admits a unique solution.

3
2.2 A general Nitsche-based finite element method
Let Vh ⊂ V be a family of finite dimensional vector spaces (see [11, 14, 7]) indexed by h
coming from a family T h of triangulations of the domain Ω (h = maxT ∈T h hT where hT is the
diameter of T ). The family of triangulations is supposed regular (i.e., there exists σ > 0 such that
∀T ∈ T h , hT /ρT ≤ σ where ρT denotes the radius of the inscribed ball in T ) and conformal to the
subdivision of the boundary into ΓD , ΓN and ΓC (i.e., a face of an element T ∈ T h is not allowed
to have simultaneous non-empty intersection with more than one part of the subdivision). To fix
ideas, we choose a standard Lagrange finite element method of degree k with k = 1 or k = 2, i.e.:
n o
Vh := vh ∈ (C 0 (Ω))d : vh|T ∈ (Pk (T ))d , ∀T ∈ T h , vh = 0 on ΓD . (4)

However, the analysis would be similar for any C 0 -conforming finite element method.
Let us introduce the notation [·]+ for the positive part of a scalar quantity a ∈ R:

a if a > 0,
[a]+ :=
0 otherwise,
The positive part satisfying
a ≤ [a]+ , a[a]+ = [a]2+ , ∀a ∈ R, (5)
we can deduce the following monotonicity property:
([a]+ − [b]+ )(a − b) = a[a]+ + b[b]+ − b[a]+ − a[b]+
≥ [a]2+ + [b]2+ − 2[a]+ [b]+ (6)
2
= ([a]+ − [b]+ ) ≥ 0.
Note that conditions (5) and (6) can be straightforwardly extended to real valued functions.
The derivation of a Nitsche-based method comes from a classical reformulation (see for instance
[2]) of the contact conditions (2) (i)-(iii), for a given γ > 0:
1
σn (u) = − [un − γσn (u)]+ . (7)
γ
Remark 2.1. Note that condition (7) is still equivalent to (2) (i)-(iii) on ΓC when γ is a positive
function defined on ΓC instead of a positive constant.
Let now θ ∈ R be a fixed parameter, and let u be the solution of the unilateral contact problem
in its strong form (1)–(2). We assume that u is sufficiently regular so that all the following
calculations make sense. From the Green formula, equations (1) and (2)-(iv), we get for every
v ∈ V: Z
a(u, v) − σn (u) vn dΓ = L(v).
ΓC
With the splitting vn = (vn − θγσn (v)) + θγσn (v), we obtain
Z Z
a(u, v) − θγ σn (u) σn (v) dΓ − σn (u) (vn − θγσn (v)) dΓ = L(v).
ΓC ΓC

Finally, using condition (7), we still have for every v ∈ V:

Z Z
1
a(u, v) − θγ σn (u) σn (v) dΓ + [un − γσn (u)]+ (vn − θγσn (v)) dΓ = L(v). (8)
ΓC ΓC γ

4
Formula (8) is the starting point of our Nitsche-based formulation. We remark that it may have no
sense at the continuous level if u lacks of regularity (the only assumption u ∈ V is not sufficient to
justify the above calculations). Nevertheless, and as in the stabilized Lagrange multiplier method
[25], we consider in what follows that γ is a positive piecewise constant function on the contact
interface ΓC defined for any x ∈ ΓC lying on the relative interior of ΓC ∩ T for a (closed) element
T having a nonempty intersection of dimension d − 1 with ΓC by

γ(x) = γ0 hT ,

where γ0 is a positive given constant (the value of γ on element intersections has no influence).
This allows to define a discrete counterpart of (8). Let us introduce for this purpose the discrete
linear operator
Vh → L2 (ΓC )
Pγ : h ,
v 7→ vn − γ σn (vh )
h

and also the bilinear form:

Z
Aθγ (uh , vh ) := a(uh , vh ) − θγ σn (uh )σn (vh ) dΓ. (9)
ΓC

Remark 2.2. With the previous notations, we see that problem (3) could be formally written as
follows:

 Find a sufficiently regular u ∈ V such that:

Z
1
 Aθγ (u, v) +
 [Pγ (u)]+ Pθγ (v) dΓ = L(v), for all sufficiently regular v ∈ V,
ΓC γ

with Pθγ (v) := vn − θγ σn (v). Note in particular that Pγ (u) and Pθγ (u) are well-defined and
3
belong to L2 (ΓC ) provided that u ∈ (H 2 +ν (Ω))d (ν > 0).

Our generalized Nitsche-based method then reads:


h h
 Find u ∈ V such that:

(10)
Z
h h 1
A
 θγ
 (u , v ) + [Pγ (uh )]+ Pθγ (vh ) dΓ = L(vh ), ∀ vh ∈ Vh .
ΓC γ

Remark 2.3. The introduction of an additional numerical parameter θ allows us to introduce

some new interesting variants acting on the symmetry / skew-symmetry / non-symmetry of the
discrete formulation. Moreover, a unified analysis of all these variants can be performed. Note
that some values of θ may be of special interest:

• for θ = 1 we recover the symmetric method proposed and analyzed in [9],

• for θ = 0 we recover a very simple non-symmetric version close to the classical penalty
method,

• for θ = −1 we obtain a skew-symmetric version which has the remarkable property to be well-
posed and convergent irrespectively of the value of γ0 > 0 (see Theorem 3.4 and Theorem
3.8).

5
Remark 2.4. As in the symmetric case [9], we can rewrite the formulation (10) in a mixed form
as follows:
Find (uh , λh ) ∈ Vh × L2− (ΓC ) such that:




 Z Z
h h h h
θγ(λh − σn (uh )) σn (vh ) dΓ = L(vh ), ∀ vh ∈ Vh ,

a(u , v ) − λ vn dΓ +

ΓC ΓC

 Z Z
(µ − λh )uhn dΓ + γ(µ − λh )(λh − σn (uh )) dΓ ≥ 0, ∀ µ ∈ L2− (ΓC ),




ΓC ΓC

where L2− (ΓC ) := {µ ∈ L2 (ΓC ) | µ ≤ 0 a.e. on ΓC } and with λh = − γ1 [Pγ (uh )]+ . Note that
formally, in the case θ 6= 1, this mixed form is somehow different from the stabilized method
applied to unilateral contact (see [25]). In particular the stabilization term in the first equation
vanishes when θ = 0.

3 Analysis of the Nitsche-based method

In this section, we carry out the mathematical analysis of the method (10). A difference between
Nitsche’s method and classical penalty methods [28, 10] is the property of consistency, which
we first show in §3.1. The proof of the well-posedness of the (non-linear) discrete problem (10)
is carried out in §3.2. The error analysis is finally detailed in §3.3. We show that the method
converges in an optimal way when the mesh size h vanishes.

3.1 Consistency
Like Nitsche’s method for second order elliptic problems with Dirichlet boundary conditions or
domain decomposition [4], our Nitsche-based formulation (10) for unilateral contact is consistent:
Lemma 3.1. The Nitsche-based method for contact is consistent: suppose that the solution u of
3
(1)–(2) lies in (H 2 +ν (Ω))d with ν > 0 and d = 2, 3. Then u is also solution of
Z
h 1
Aθγ (u, v ) + [Pγ (u)]+ Pθγ (vh ) dΓ = L(vh ), ∀ vh ∈ Vh .
ΓC γ
3
Proof: Let u be the solution of (1)–(2) and set vh ∈ Vh . Since u ∈ (H 2 +ν (Ω))d and ν > 0, we
have σn (u) ∈ H ν (ΓC ) ⊂ L2 (ΓC ). As a result, Aθγ (u, vh ) makes sense and Pγ (u) ∈ L2 (ΓC ). On
the one hand, we use the definition of Pγ , of Aθγ (·, ·) and the reformulation (7) of the contact
conditions to obtain:
Z
h 1
Aθγ (u, v ) + [Pγ (u)]+ Pθγ (vh ) dΓ
ΓCγ
Z Z
1
= a(u, vh ) − θγ σn (u)σn (vh ) dΓ + (−γσn (u))(vnh − θγ σn (vh )) dΓ
γ
ZΓC ΓC

= a(u, vh ) − σn (u)vnh dΓ.

ΓC

On the other hand, with equations (1)–(2) and an integration by parts, we get:
Z
a(u, vh ) − σn (u)vnh dΓ = L(vh ),
ΓC

which ends the proof. 

6
3.2 Well-posedness
To prove well-posedness of our Nitsche-based formulation, we first need the following classical
property.
Lemma 3.2. There exists C > 0, independent of the parameter γ0 and of the mesh size h, such
that:
1
kγ 2 σn (vh )k20,ΓC ≤ Cγ0 kvh k21,Ω , (11)
for all vh ∈ Vh .
Proof: It follows from the definition of σn (vh ) and the boundedness of A that:
1
kγ 2 σn (vh )k20,ΓC ≤ γ0 hkσn (vh )k20,ΓC ≤ Cγ0 hk∇vh k20,ΓC .
Then estimation (11) is obtained using a scaling argument: see [38, Lemma 2.1, p.24] for a
detailed proof in the general case (for an arbitrary degree k and any dimension d). 
Remark 3.3. Note that contrary to the penalty case or to the symmetric case of Nitsche’s method
(i.e., θ = 1), the integral term on ΓC :
Z
1
[Pγ (uh )]+ Pθγ (vh )
ΓC γ

is not necessarily non-negative when vh = uh and θ 6= 1.

We then show that Problem (10) is well-posed using an argument by Brezis for M-type and
pseudo-monotone operators [8] (see also [32] and [29]).
Theorem 3.4. Suppose that one of the two following assumptions hold:
1. θ 6= −1 and γ0 > 0 is sufficiently small,
2. θ = −1 and γ0 > 0.
Then Problem (10) admits one unique solution uh in Vh .
Proof: Using the Riesz representation theorem, we define a (non-linear) operator Bh : Vh → Vh ,
by means of the formula:
Z
1
(Bh vh , wh )1,Ω := Aθγ (vh , wh ) + [Pγ (vh )]+ Pθγ (wh ) dΓ,
ΓC γ
for all vh , wh ∈ Vh , and where (·, ·)1,Ω stands for the scalar product in (H 1 (Ω))d . Note that
Problem (10) is well-posed if and only if Bh is a one-to-one operator.
Let vh , wh ∈ Vh . We have:
1
(Bh vh − Bh wh , vh − wh )1,Ω = a(vh − wh , vh − wh ) − θkγ 2 σn (vh − wh )k20,ΓC
Z
1
+ ([Pγ (vh )]+ − [Pγ (wh )]+ )(vnh − wnh − θγσn (vh − wh )) dΓ
ΓC γ
1
= a(vh − wh , vh − wh ) − θkγ 2 σn (vh − wh )k20,ΓC
Z
1
+ ([Pγ (vh )]+ − [Pγ (wh )]+ )Pγ (vh − wh ) dΓ
ΓC γ
Z
1
+ (1 − θ) ([Pγ (vh )]+ − [Pγ (wh )]+ )γσn (vh − wh ) dΓ.
ΓC γ
(12)

7
Using (6) in (12) and Cauchy-Schwarz inequality, we get

(Bh vh − Bh wh , vh − wh )1,Ω
1
≥ a(vh − wh , vh − wh ) − θkγ 2 σn (vh − wh )k20,ΓC
1 (13)
+ kγ − 2 ([Pγ (vh )]+ − [Pγ (wh )]+ )k20,ΓC
1 1
− |1 − θ| kγ − 2 ([Pγ (vh )]+ − [Pγ (wh )]+ )k0,ΓC kγ 2 σn (vh − wh )k0,ΓC .

If θ = 1, we use the coercivity of a(·, ·) and the property (11) in the previous expression (13).
Therefore:
1
(Bh vh − Bh wh , vh − wh )1,Ω ≥ a(vh − wh , vh − wh ) − kγ 2 σn (vh − wh )k20,ΓC
(14)
≥ Ckvh − wh k21,Ω ,

when γ0 is chosen sufficiently small. This corresponds to the symmetric version already studied
in [9].
We now suppose that θ 6= 1. Let β > 0. Applying Young inequality in (13) yields:

(Bh vh − Bh wh , vh − wh )1,Ω
1
≥ a(vh − wh , vh − wh ) − θkγ 2 σn (vh − wh )k20,ΓC
1
+ kγ − 2 ([Pγ (vh )]+ − [Pγ (wh )]+ )k20,ΓC
|1 − θ| − 1 |1 − θ|β 1
− kγ 2 ([Pγ (vh )]+ − [Pγ (wh )]+ )k20,ΓC − kγ 2 σn (vh − wh )k20,ΓC (15)
2β 2
 
h h h h |1 − θ|β 1
= a(v − w , v − w ) − θ + kγ 2 σn (vh − wh )k20,ΓC
2
 
|1 − θ| 1
+ 1− kγ − 2 ([Pγ (vh )]+ − [Pγ (wh )]+ )k20,ΓC .
2β
Choosing β = |1 − θ|/2 in (15), we get:
1 1
(Bh vh − Bh wh , vh − wh )1,Ω ≥ a(vh − wh , vh − wh ) − (1 + θ)2 kγ 2 σn (vh − wh )k20,ΓC
4 (16)
h
≥ Ckv − wh k21,Ω ,

when θ 6= −1 and γ0 sufficiently small, or when θ = −1. Note that in the latter case we do not
need the smallness assumption on γ0 .
Next, let us show that Bh is also hemicontinuous. Since Vh is a vector space, it is sufficient to
show that
[0, 1] 3 t 7→ ϕ(t) := (Bh (vh − twh ), wh )1,Ω ∈ R
is a continuous real function, for all vh , wh ∈ Vh . Let s, t ∈ [0, 1], we have:

|ϕ(t) − ϕ(s)|
= |(Bh (vh − twh ) − Bh (vh − swh ), wh )1,Ω |
Z
1  
= Aθγ ((s − t)wh , wh ) + [Pγ (vh − twh )]+ − [Pγ (vh − swh )]+ Pθγ (wh ) dΓ
γ
Z ΓC
1
≤ |s − t| Aθγ (wh , wh ) + [Pγ (vh − twh )]+ − [Pγ (vh − swh )]+ |Pθγ (wh )| dΓ.
ΓC γ

8
With help of the bound |[a]+ − [b]+ | ≤ |a − b| , for all a, b ∈ R, and using the linearity of Pγ , we
deduce that:
Z
1
[Pγ (vh − twh )]+ − [Pγ (vh − swh )]+ |Pθγ (wh )| dΓ
ΓC γ
Z
1
≤ Pγ (vh − twh ) − Pγ (vh − swh ) |Pθγ (wh )| dΓ
γ
ZΓC
1
= |(s − t)Pγ (wh )||Pθγ (wh )| dΓ.
ΓC γ

It results that:
 Z 
1
|ϕ(t) − ϕ(s)| ≤ |s − t| Aθγ (wh , wh ) + |Pγ (wh )||Pθγ (wh )| dΓ ,
ΓC γ

which means that ϕ is Lipschitz, so that Bh is hemicontinuous. Since properties (14), (16) also
hold, we finally apply the Corollary 15 (p.126) of [8] to conclude that Bh is a one-to-one operator.
This ends the proof. 
Remark 3.5. When γ0 is large and θ 6= −1 we can neither conclude to uniqueness, nor to
existence of a solution. In Appendix B, we show some simple explicit examples of nonexistence
and nonuniqueness of solutions.

3.3 A priori error analysis

Our Nitsche-based method (10) converges in a optimal way as the mesh parameter h vanishes.
This is proved in the following theorems. First, we establish an abstract error estimate.
3
Theorem 3.6. Suppose that the solution u of Problem (3) belongs to (H 2 +ν (Ω))d with ν > 0
and d = 2 or d = 3.
1. Let θ ∈ R. Suppose that the parameter γ0 > 0 is sufficiently small. Then, the solution uh of
Problem (10) satisfies the following abstract error estimate:
1 1
ku − uh k1,Ω + kγ 2 (σn (u) + [Pγ (uh )]+ )k0,ΓC
γ
 1 1
 (17)
≤ C inf ku − vh k1,Ω + kγ − 2 (un − vnh )k0,ΓC + kγ 2 σn (u − vh )k0,ΓC ,
vh ∈Vh

where C is a positive constant, independent of h, u and γ0 .

2. Set θ = −1. Then for all values of γ0 > 0, the solution uh of Problem (10) satisfies the
abstract error estimate (17) where C is a positive constant, dependent of γ0 but still independent
of h and u.
Proof: Let vh ∈ Vh . We first use the V-ellipticity and the continuity of a(·, ·), as well as Young’s
inequality, to obtain:

αku − uh k21,Ω ≤ a(u − uh , u − uh )

= a(u − uh , (u − vh ) + (vh − uh ))
≤ Cku − uh k1,Ω ku − vh k1,Ω + a(u − uh , vh − uh )
α C2
≤ ku − uh k21,Ω + ku − vh k21,Ω + a(u, vh − uh ) − a(uh , vh − uh ),
2 2α

9
where α > 0 is the ellipticity constant of a(., .). We can transform the term a(u, vh − uh ) −
a(uh , vh − uh ) since u solves (3), uh solves (10) and using Lemma 3.1, which yields:
C2
Z
α h 2 h 2
ku − u k1,Ω ≤ ku − v k1,Ω − θ γ σn (uh − u)σn (vh − uh ) dΓ
2 2α ΓC
Z
1 
+ [Pγ (u )]+ − [Pγ (u)]+ Pθγ (vh − uh ) dΓ.
h
(18)
ΓC γ

The first integral term in (18) is bounded as follows:

Z
− θ γ σn (uh − u)σn (vh − uh ) dΓ
Z ΓC Z
h h h h
= θ γ σn (v − u )σn (v − u ) dΓ − θ γ σn (vh − u)σn (vh − uh ) dΓ
ΓC ΓC
1 1 1
h h
≤ θkγ σn (v − u
2 )k20,ΓC
+ |θ|kγ σn (v − u)k0,ΓC kγ 2 σn (vh − uh )k0,ΓC
2
h

β1 θ 2 1
 
h 2 1 1
≤ kγ 2 σn (v − u)k0,ΓC + + θ kγ 2 σn (vh − uh )k20,ΓC , (19)
2 2β1
with β1 > 0. The second integral term in (18) is considered next:
Z
1 
[Pγ (uh )]+ − [Pγ (u)]+ Pθγ (vh − uh ) dΓ
γ
ZΓC  Z
1  1 
= [Pγ (uh )]+ − [Pγ (u)]+ Pγ (vh − u) dΓ + [Pγ (uh )]+ − [Pγ (u)]+ Pγ (u − uh ) dΓ
ΓC γ ΓC γ
Z
1  
+ [Pγ (uh )]+ − [Pγ (u)]+ γ(1 − θ)σn (vh − uh ) dΓ. (20)
ΓC γ

Using the bound ([a]+ − [b]+ )(b − a) ≤ −([a]+ − [b]+ )2 (see (6)) and two times Cauchy-Schwarz
and Young’s inequalities in (20) we obtain
Z
1 
[Pγ (uh )]+ − [Pγ (u)]+ Pθγ (vh − uh ) dΓ
ΓC γ
1 1 1 β2 1 1 1
≤ kγ 2 (σn (u) + [Pγ (uh )]+ )k20,ΓC + kγ − 2 Pγ (vh − u)k20,ΓC − kγ 2 (σn (u) + [Pγ (uh )]+ )k20,ΓC
2β2 γ 2 γ
|1 − θ| 1 1 |1 − θ|β 3 1
+ kγ 2 (σn (u) + [Pγ (uh )]+ )k20,ΓC + kγ 2 σn (vh − uh )k20,ΓC , (21)
2β3 γ 2
with β2 > 0 and β3 > 0. Noting that
1 1 1
kγ − 2 Pγ (vh − u)k20,ΓC ≤ 2kγ − 2 (un − vnh )k20,ΓC + 2kγ 2 σn (u − vh )k20,ΓC

and putting together estimates (19) and (21) in (18) gives

C2 β1 θ 2
 
α 1 1
h 2
ku − u k1,Ω ≤ h 2
ku − v k1,Ω + + β2 kγ 2 σn (u − vh )k20,ΓC + β2 kγ − 2 (un − vnh )k20,ΓC
2 2α 2
 
1 |1 − θ| 1 1
+ −1 + + kγ 2 (σn (u) + [Pγ (uh )]+ )k20,ΓC
2β2 2β3 γ
 
1 |1 − θ|β3 1
+ +θ+ kγ 2 σn (vh − uh )k20,ΓC . (22)
2β1 2

10
The norm term on the discrete normal constraints in (22) is bounded as follows by using (11):
1 1 1
kγ 2 σn (vh − uh )k0,ΓC ≤ Cγ02 kvh − uh k1,Ω ≤ Cγ02 (kvh − uk1,Ω + kuh − uk1,Ω ). (23)

For a fixed θ ∈ R we then choose β2 and β3 large enough such that

1 |1 − θ| 1
−1 + + <− .
2β2 2β3 2
Choosing then γ0 small enough in (23) and putting the estimate in (22) establishes the first
statement of theorem.
We now consider separately the case θ = −1 for which (22) becomes:

C2
 
α β1 1 1
h 2
ku − u k1,Ω ≤ h 2
ku − v k1,Ω + + β2 kγ 2 σn (u − vh )k20,ΓC + β2 kγ − 2 (un − vnh )k20,ΓC
2 2α 2
 
1 1 1 1
+ −1 + + kγ 2 (σn (u) + [Pγ (uh )]+ )k20,ΓC
2β2 β3 γ
 
1 1
+ − 1 + β3 kγ 2 σn (vh − uh )k20,ΓC . (24)
2β1
Let be given η > 0. Set β1 = 1/(2η), β2 = 1 + (1/η) and β3 = 1 + η. Therefore (24) becomes:

C2
 
α h 2 h 2 5 1 1 + η −1
ku − u k1,Ω ≤ ku − v k1,Ω + + 1 kγ 2 σn (u − vh )k20,ΓC + kγ 2 (un − vnh )k20,ΓC
2 2α 4η η
η 1 1 1
− kγ 2 (σn (u) + [Pγ (uh )]+ )k20,ΓC + 2ηkγ 2 σn (vh − uh )k20,ΓC .
2(1 + η) γ

Let γ0 > 0. Set η = α/(16C 2 γ0 ) where C is the constant in (23) to conclude the proof of the
theorem. 
Remark 3.7. In the case θ = −1, note that the convergence result holds for any value of γ0 > 0.
However, in that case, the abstract estimate takes the form
 1
h C 2 1 1
ku − u k1,Ω + kγ 2 (σn (u) + [Pγ (uh )]+ )k0,ΓC
γ0 + C γ
1
 1 1

−
≤ C(1 + γ0 ) 2 inf ku − v k1,Ω + kγ 2 (un − vnh )k0,ΓC + kγ 2 σn (u − vh )k0,ΓC ,
h
vh ∈Vh

where C is a positive constant, independent of h, u and γ0 . As a result, the estimation is

deteriorated when γ0 increases, in particular the estimation of the contact stress on ΓC .
The optimal convergence of the method is stated below.
3
Theorem 3.8. Suppose that the solution u to Problem (3) belongs to (H 2 +ν (Ω))d with 0 < ν ≤
k− 21 (k = 1, 2 is the degree of the finite element method, given in (4)) and d = 2, 3. When θ 6= −1,
suppose in addition that the parameter γ0 is sufficiently small. The solution uh of Problem (10)
satisfies the following error estimate:
1 1 1
ku − uh k1,Ω + kγ 2 (σn (u) + [Pγ (uh )]+ )k0,ΓC ≤ Ch 2 +ν kuk 3 +ν,Ω , (25)
γ 2

where C is a positive constant, dependent of γ0 but independent of h and u.

11
Proof: We need to bound the right terms in estimate (17) and we choose vh = I h u where I h
stands for the Lagrange interpolation operator mapping onto Vh . The estimation of the Lagrange
interpolation error in the H 1 -norm on a domain Ω is classical (see, e.g., [7, 13, 14]):
1
ku − I h uk1,Ω ≤ Ch 2 +ν kuk 3 +ν,Ω , (26)
2

for − 21 < ν ≤ k − 12 .
1
The estimation of the term kγ − 2 (un − (I h u)n )k0,ΓC can be done in a very similar manner to [25].
Indeed, let E be an edge (resp. a face) of a triangle (resp. tetrahedron) T ∈ T h on ΓC :
1 −1 1
kγ − 2 (un − (I h u)n )k0,E ≤ ChT 2 h1+ν
T kun k1+ν,E = Ch
2
+ν
kun k1+ν,E ,
(see [13] for instance). By summation on all the edges (resp. faces) and the trace theorem, it
results:
1 1 1
kγ − 2 (un − (I h u)n )k0,ΓC ≤ Ch 2 +ν kun k1+ν,ΓC ≤ Ch 2 +ν kuk 3 +ν,Ω . (27)
2

From Lemma A.1 in appendix, (see also [16, 17]), the following estimate also holds:
1 1
kγ 2 σn (u − I h u)k0,ΓC ≤ Ch 2 +ν kuk 3 +ν,Ω . (28)
2

We conclude by inserting the three estimates (26)–(28) into (17). 

Remark 3.9. Note that for the theoretical estimation of Theorem 3.8, in the case θ 6= −1, and
−1 1
for γ0 small, the constant C in the estimate (25) behaves in O(γ0 2 ) (due to the term kγ − 2 (un −
vnh )k0,ΓC ), so that taking γ0 too small theoretically deteriorates the convergence. In the case
θ = −1, the same observation can be made, with the additional fact that the convergence is also
deteriorated in O(γ0 ) if γ0 is too large (see Remark 3.7). However, numerical experiments of
Section 4 reveal that in practice the method is quite robust relatively to the value of γ0 , provided
it is sufficiently small when θ 6= −1.
Actually we are not able to obtain estimates for the displacements in the L2 -norm (ku − uh k0,Ω
and also ku − uh k0,ΓC ) by using the Aubin-Nitsche argument as it is achieved in the linear
case (see [17]). Note also that the L2 (Ω)-norm estimates for contact problems are not easy
to prove and there are to our knowledge only few estimates (see [12]). Nevertheless, we can
easily obtain the following error estimate on the weighted L2 (ΓC )-norm on the normal constraint
1
kγ 2 σn (u − uh )k0,ΓC (note that σn (uh ) 6= − γ1 [Pγ (uh )]+ on ΓC contrary to the continuous case).
3
Corollary 3.10. Suppose that the solution u to Problem (3) belongs to (H 2 +ν (Ω))d with 0 <
ν ≤ k − 21 (k = 1, 2 is the degree of the finite element method, given in (4)) and d = 2, 3.
When θ 6= −1, suppose in addition that the parameter γ0 is sufficiently small. The solution uh of
Problem (10) satisfies the following error estimate:
1 1
kγ 2 σn (u − uh )k0,ΓC ≤ Ch 2 +ν kuk 3 +ν,Ω ,
2

where C is a positive constant, independent of h and u.

Proof: We use (28), (11), (26) and (25) to establish the bound:
1 1 1
kγ 2 σn (u − uh )k0,ΓC ≤ kγ 2 σn (u − I h u)k0,ΓC + kγ 2 σn (I h u − uh )k0,ΓC
1 1
≤ Ch 2 +ν kuk 3 +ν,Ω + Cγ02 (kI h u − uk1,Ω + ku − uh k1,Ω )
2
1
+ν
≤ Ch 2 kuk 3 +ν,Ω .
2

12
Remark 3.11. Although contact problems are known to be limited in regularity (the regularity
H 5/2 (Ω) can not generally be passed beyond for such inequality problems), the use of quadratic
finite element methods can be of interest in particular for the most regular solutions lying in
H s (Ω), 2 < s < 5/2 (as it is considered in e.g., [5, 24, 35, 40]).

4 Numerical experiments
In this section, the numerical results of two and three-dimensional Hertz’s contact problems of a
disk/sphere with a plane rigid foundation are presented. This slightly exceeds the scope defined
in Section 2 since a non-zero initial gap between the elastic solid and the rigid foundation is
considered in the computations. Moreover, the tests are performed with P1 and isoparametric P2
Lagrange finite elements on meshes which are approximations of the real domain.
The finite element library Getfem++1 is used. The discrete contact problem is solved by using
a generalized Newton method, which means that Problem (10) is derived with respect to uh to
obtain the tangent system. The term “generalized Newton’s method” comes from the fact that
the positive part [x]+ is non-differentiable at x = 0. However, no special treatment is considered.
If a point of non-differentiability is encountered, the tangent system corresponding to one of the
two alternatives (x < 0 or x > 0) is chosen arbitrarily. Integrals of the non-linear term on ΓC
are computed with standard quadrature formulas. Note that the situation where the solution is
non-differentiable at an integration point is very rare and corresponds to what is called a “grazing
contact” (both un = 0 and σn (u) = 0). Further details on generalized Newton’s method applied
to contact problems can be found for instance in [36] and the references therein.

4.1 Two-dimensional numerical tests

The numerical situation is represented in Fig. 1. A disc of radius 20 cm is considered with a
contact boundary ΓC which is restricted to the lower part (y < 20 cm) of the boundary. A homo-
geneous Neumann condition is applied on the remaining part of the boundary. Since no Dirichlet
condition is considered, the problem is not fully coercive. To overcome the non-definiteness com-
ing from the free rigid motions, the horizontal displacement is prescribed to be zero on the two
points of coordinates (0 cm, 10 cm) and (0 cm, 30 cm) which blocks the horizontal translation and
the rigid rotation. Homogeneous isotropic linear elasticity in plane strain approximation is con-
sidered with a Young modulus fixed at E = 25 MPa and a Poisson ratio P = 0.25. A vertical
density of volume forces of 20 MN/m3 is applied.
1
see [Link]

13
Figure 1: Example of two-dimensional mesh and reference solution (with color plot of the Von-
Mises stress).

The error curves when using a P1 Lagrange finite element method are shown in Figs. 2, 3 and
4. The solution for mesh sizes h = 0.5 cm, 1 cm, 3 cm, 4.5 cm and h = 10 cm are compared with a
reference solution on a very fine mesh (h = 0.15 cm) using quadratic isoparametric finite elements.
Moreover, the reference solution is computed with a different discretization of the contact problem
(Lagrange multipliers and Alart-Curnier augmented Lagrangian, see [36]). On all the figures, the
graph on the left represents the relative H 1 (Ω)-norm of the error between the computed solution
and the reference one and the graph on the right represents the following relative L2 (ΓC )-norm:
1
kγ 2 (γ −1 [Pγ (uh )]+ − γ −1 [Pγ (uhref )]+ )k0,ΓC
,
kγ −1 [Pγ (uhref )]+ k0,ΓC
where uh is the discrete solution and uhref the reference solution. The slopes shown in the
figures give an approximation of the convergence rate. They correspond to the slopes of the
regression lines on a logarithmic scale for the different experiments. Note that − γ1 [Pγ (uh )]+ is
an approximation of the contact stress with a convergence of order 1 (see Theorem 3.8).
Let us recall that the formulation is symmetric if and only if θ = 1. The results corresponding
to this symmetric case are shown in Fig. 2. Optimal convergence is obtained for both H 1 (Ω)
and weighted L2 (ΓC )-norms of the error, but only for the smallest value of the parameter γ0
(γ0 = 1/(100E)). For higher values of γ0 , at least the convergence of the contact stress is non-
optimal. This corroborates the theoretical result of Theorem 3.8 for which the optimal rate of
convergence is obtained for a sufficiently small γ0 . We also noted that for the largest value of
γ0 (i.e., γ0 = 100/E) the convergence of Newton’s method is not always achieved (in this case,
Newton’s method is stopped after 100 iterations). In some experiments, the residual of Newton’s
method diminishes to a value which is greatly higher than the one we considered for the tests
achieving the convergence. Our interpretation is that for large values of γ0 , when coercivity is
lost, there might be no solution to the discrete problem.
Now, when θ = 0, the error curves are plotted in Fig. 3. For this version, which is in a
sense simpler to implement than the other ones (since the number of additional terms is lower),

14
Newton’s method always converges even for the largest value of γ0 . For this latter value of γ0
the convergence remains sub-optimal. However, for the intermediate value of γ0 (γ0 = 1/E) the
optimal convergence is reached.
Concerning the version with θ = −1, which corresponds to an unconditionally coercive problem,
one can see in Fig. 4 that optimal convergence is reached for all values of γ0 . Moreover, the
smallest error on the contact stress corresponds to the intermediate value of γ0 .
Globally, it is remarkable that the error curves for the smallest value of γ0 are rather the same for
the three values of θ. A strategy to guarantee an optimal convergence is of course to consider a
sufficiently small γ0 . However, the price to pay is an ill-conditioned discrete problem. The study
presented in [36] for the versions θ = 1 and θ = 0 shows that Newton’s method has important
difficulties to converge when γ0 is small. When symmetry is not required, a better strategy seems
to consider the version with θ = −1 or an intermediate value of θ = 0 which ensure both a optimal
convergence rate and few iterations of Newton’s method to converge.
We have no interpretation of the slight super-convergence noted on most of the error curves.
Theoretically, the convergence rates should be close to 1, since, if we assume that there are only
two transition points between effective contact and non-contact, the solution to the continuous
problem should be in (H s (Ω))2 for all s < 5/2 (see [33]).

γ0=1/100E (slope=1.1952) γ0=1/100E (slope=1.4331)

γ0=1/E (slope=0.8438) γ0=1/E (slope=0.60812)
L2(ΓC) contact stress relative error

100%
γ0=100/E (slope=1.3522) γ0=100/E (slope=0.52918)
H1(Ω) relative error

100%

1%

1%
0.01%

1 10 1 10
h h

Figure 2: Error curves in the 2D case with P1 elements and θ = 1. Left: relative H 1 (Ω)-norm on
the displacements. Right: relative weighted L2 (ΓC )-norm on the contact pressure.

15
 γ0=1/100E (slope=1.2026) γ0=1/100E (slope=1.4334)
γ0=1/E (slope=1.2286) γ0=1/E (slope=1.3469)

L2(ΓC) contact stress relative error

100%
γ0=100/E (slope=0.8851) γ0=100/E (slope=0.8151)
H1(Ω) relative error

100%

1%

1%
0.01%

1 10 1 10
h h

Figure 3: Error curves in the 2D case with P1 elements and θ = 0. Left: relative H 1 (Ω)-norm on
the displacements. Right: relative weighted L2 (ΓC )-norm on the contact pressure.

γ0=1/100E (slope=1.2101) γ0=1/100E (slope=1.4338)

γ0=1/E (slope=1.3189) γ0=1/E (slope=1.3214)
L2(ΓC) contact stress relative error

100%
γ0=100/E (slope=1.6323) γ0=100/E (slope=1.4662)
H1(Ω) relative error

100%

1%

1%
0.01%

1 10 1 10
h h

Figure 4: Error curves in the 2D case with P1 elements and θ = −1. Left: relative H 1 (Ω)-norm
on the displacements. Right: relative weighted L2 (ΓC )-norm on the contact pressure.
The same numerical experiment has been extended to the P2 Lagrange isoparametric finite ele-
ment method. The corresponding error curves are presented in Figs. 5, 6 and 7. The results are
quite similar compared with the P1 Lagrange method. The convergence is even poorer for large
values of γ0 and θ = 1. The error levels are smaller compared with the P1 method. However,
the convergence rates are only slightly better. This comes probably from the regularity of the
solution to the continuous problem.

16
 γ0=1/100E (slope=1.6153) γ0=1/100E (slope=1.4477)
γ0=1/E (slope=0.13572) γ0=1/E (slope=0.43804)

L2(ΓC) contact stress relative error

100%
γ0=100/E (slope=0.99543) γ0=100/E (slope=0.58009)
H1(Ω) relative error

100%

1%

1%
0.01%

1 10 1 10
h h

Figure 5: Error curves in the 2D case with P2 elements and θ = 1. Left: relative H 1 (Ω)-norm on
the displacements. Right: relative weighted L2 (ΓC )-norm on the contact pressure.

γ0=1/100E (slope=1.6285) γ0=1/100E (slope=1.4198)

γ0=1/E (slope=1.6252) γ0=1/E (slope=1.5008)
L2(ΓC) contact stress relative error

100%
γ0=100/E (slope=1.4303) γ0=100/E (slope=1.2389)
H1(Ω) relative error

100%

1%

1%
0.01%

1 10 1 10
h h

Figure 6: Error curves in the 2D case with P2 elements and θ = 0. Left: relative H 1 (Ω)-norm on
the displacements. Right: relative weighted L2 (ΓC )-norm on the contact pressure.

17
 γ0=1/100E (slope=1.6376) γ0=1/100E (slope=1.4161)
γ0=1/E (slope=1.7495) γ0=1/E (slope=1.5539)

L2(ΓC) contact stress relative error

100%
γ0=100/E (slope=1.9389) γ0=100/E (slope=1.612)
H1(Ω) relative error

100%

1%

1%
0.01%

1 10 1 10
h h

Figure 7: Error curves in the 2D case with P2 elements and θ = −1. Left: relative H 1 (Ω)-norm
on the displacements. Right: relative weighted L2 (ΓC )-norm on the contact pressure.

4.2 Three-dimensional numerical tests

The three-dimensional tests are similar to the two-dimensional ones. We consider a sphere of
radius 20 cm with mesh sizes h = 3.6 cm, 6 cm, 11 cm, 23 cm and a P1 Lagrange finite element
method (an example of a mesh and a reference solution are presented in Fig. 8). Homogeneous
isotropic linear elasticity is considered with still a Young modulus E = 25 MPa and a Poisson ratio
P = 0.25. A vertical density of volume forces of 20 MN/m3 is also still considered. Similarly to the
two-dimensional case, the horizontal rigid motions and the rotations are blocked by prescribing
the displacement on specific chosen points. The reference solution is still computed with quadratic
isoparametric finite elements on a fine mesh (h = 1 cm) using Lagrange multipliers. The error
curves are presented in Fig. 9, 10 and 11. Very similar conclusions can be drawn compared with
the two-dimensional case. In particular, the method for θ = −1 allows us to obtain an optimal
convergence rate for any value of γ0 . However, note that the convergence for large values of γ0
for θ = 1 and θ = 0 is even worse compared to the two-dimensional case.

18
Figure 8: Three-dimensional mesh example and reference solution (sectional view, with color plot
of the Von-Mises stress).

γ0=1/100E (slope=1.615) γ0=1/100E (slope=1.1246)

γ0=1/E (slope=−0.21414) γ0=1/E (slope=0.80868)
L2(ΓC) contact stress relative error

γ0=100/E (slope=−0.46543) γ0=100/E (slope=0.51169)

H1(Ω) relative error

100%
100%

1%

1%
1 10 1 10
h h

Figure 9: Error curves in the 3D case with P1 elements and θ = 1. Left: relative H 1 (Ω)-norm on
the displacements. Right: relative weighted L2 (ΓC )-norm on the contact pressure.

19
 γ0=1/100E (slope=1.5876) γ0=1/100E (slope=1.097)
γ0=1/E (slope=1.0013) γ0=1/E (slope=1.9446)

L2(ΓC) contact stress relative error

γ0=100/E (slope=0.40117) γ0=100/E (slope=1.2549)
H1(Ω) relative error

100%
100%

1%

1%
1 10 1 10
h h

Figure 10: Error curves in the 3D case with P1 elements and θ = 0. Left: relative H 1 (Ω)-norm
on the displacements. Right: relative weighted L2 (ΓC )-norm on the contact pressure.

γ0=1/100E (slope=1.5604) γ0=1/100E (slope=1.0709)

γ0=1/E (slope=1.438) γ0=1/E (slope=1.8505)
L2(ΓC) contact stress relative error

γ0=100/E (slope=1.4222) γ0=100/E (slope=1.4133)

H1(Ω) relative error

100%
100%

1%

1%
1 10 1 10
h h

Figure 11: Error curves in the 3D case with P1 elements and θ = −1. Left: relative H 1 (Ω)-norm
on the displacements. Right: relative weighted L2 (ΓC )-norm on the contact pressure.

5 Conclusion and perspectives

In comparison with the other existing methods for unilateral contact, Nitsche’s method has the
great advantage that no other unknown than the displacement field is introduced. Indeed, most
of the methods are based on a mixed formulation which introduces Lagrange multipliers: see e.g.
[39] for a recent review and [36] for a recent numerical comparison of existing methods. From this
point of view, Nitsche’s method is close to the well-known penalty method in its simplicity, with
the advantage of remaining consistent. Moreover, it admits some variants that have interesting
properties from the theoretical and/or computational point of view: the variant θ = 1 preserves
symmetry and has a positive penalty term, the variant θ = 0 involves less terms and the variant
θ = −1 preserves well-posedness and the convergence irrespectively of the parameter γ0 . For all
the variants, i.e., all the values of θ, the discrete problem is well-posed and optimal convergence
is achieved if γ0 is sufficiently small when θ 6= −1. Numerical experiments point out the good
convergence properties of Nitsche’s formulation when solved with a generalized Newton’s method.

20
In particular, the variant θ = −1 presents very good convergence properties, for both small and
large values of γ0 .
Among future directions of work is extension of the method for friction problems, in particular
Tresca’s friction and then Coulomb’s friction.

Appendix A: interpolation error estimate for the gradient on the

boundary
3
Lemma A.1. Let u ∈ H 2 +ν (Ω) where Ω ⊂ Rd and 0 < ν ≤ k − 12 (k = 1, 2 is the degree of the
finite element method, given in (4)). Set Γ := ∂Ω and consider an element T ∈ T h such that
Γ ∩ T is a face of T . Then there exists a positive constant C independent of T and h such that:

k∇(u − I h u)k0,Γ∩T ≤ ChνT |u| 3 +ν,T ,

2

3
where hT is the diameter of T and |u| 3 +ν,T is the usual H 2 +ν (T ) seminorm of u.
2

Proof: We will use a classical scaling argument (see for instance [14, Theorem 1.103]). Let us
consider the reference element T̂ (which is independent of T and hT ) and the jacobian matrix JT
of the linear geometric transformation from T̂ to T . Then, due to the regularity of the family of
meshes T h , we have
|T | hT h
| det(JT )| = , kJT k ≤ , kJT−1 k ≤ T̂ ,
|T̂ | ρT̂ ρT
where kJT k := sup (kJT x̂k/kx̂k) is the matrix norm associated to the usual euclidean norm in Rd
x̂6=0
and |T |, |T̂ | stand for the areas of T , T̂ . Using this and the regularity of the mesh, we deduce
from a basic calculus
d−3
ˆ − Î h û)k ,
k∇(u − I h u)k0,Γ∩T ≤ ChT 2 k∇(û 0,Γ̂

where Γ̂ is the corresponding face on the reference element, ∇, ˆ Î h are the gradient and the
Lagrange interpolation operator in the reference coordinates, respectively, and û(x̂) = u(JT (x̂)).
Now, we consider the map
3
F : H 2 +ν (T̂ ) → (L2 (Γ̂))d ,
ˆ − Î h û).
û 7→ ∇(û
From standard trace theorems (see [1]), we deduce that F is continuous. Using the property
F (p̂) = 0 for all p̂ ∈ Pk (T̂ ), we can write
ˆ − Î h û)k
k∇(û 0,Γ̂ = kF (û + p̂)k0,Γ̂ ∀p̂ ∈ Pk (T̂ ),
≤ kF k 3 kû + p̂k 3 +ν,T̂ ∀p̂ ∈ Pk (T̂ ),
L (H 2 +ν (T̂ ),(L2 (Γ̂))d ) 2

≤ C inf kû + p̂k 3 +ν,T̂

p̂∈Pk (T̂ ) 2

≤ C|û| 3 +ν,T̂ .
2

The last estimate is the application of the extension to fractional order spaces of Deny-Lions
lemma. Such an extension can be found for instance in [13] (Theorem 6.1). Now, proceeding to

21
the reverse change of variable, and still using the regularity of the mesh and the expression of
|û| 3 +ν,T̂ given also for instance in [13] it comes
2

d−3
k∇(u − I h u)k0,Γ∩T ≤ ChT 2 |û| 3 +ν,T̂
2
d−3
− d−3−2ν
≤ ChT 2
hT 2 |u| 3 +ν,T
2

= ChνT |u| 3 +ν,T .

2


This result can be straightforwardly extended to the vectorial case. The global interpolation
estimate on the whole ΓC can be obtained by summation on all the faces of elements lying on
ΓC .

Appendix B: a simple example of nonexistence and nonuniqueness

of solutions
We consider the triangle Ω of vertexes A = (0, 0), B = (`, 0) and C = (0, `). We define ΓD =
[B, C], ΓN = [A, C], ΓC = [A, B] and {X1 , X2 } denotes the canonical orthonormal basis (see
Figure 12). We suppose that the volume forces f are absent and that the surface forces denoted
g = g1 X1 + g2 X2 are such that g1 and g2 are constant on ΓN .

C
X2

g
X1
ΓD
ΓN
Ω
t
A ΓC B
n

Figure 12: Setting of the problem.

We suppose that Ω is discretized with a single finite element of degree one. Consequently, the
finite element space becomes:
n o
Vh := vh = (v1h , v2h ) ∈ (P1 (Ω))2 , vh = 0 on ΓD .

Clearly, Vh is of dimension two. For vh ∈ Vh (resp. uh ∈ Vh ), we denote by (VT , VN ) (resp.

(UT , UN )) the value of vh (A) corresponding to the tangential and the normal displacements at
point A respectively (in our example, we have VT = −v1h (A) and VN = −v2h (A)). Then, for any
vh ∈ Vh , one obtains  
1 2VT VT + VN
ε(vh ) = .
2` VT + VN 2VN

22
We denote by E and ν the Young modulus and the Poisson ratio or equivalently λ = (EP )/((1 −
2P )(1 + P )) and µ = E/(2(1 + P )) are the corresponding Lamé coefficients. So we get
 
h 1 (λ + 2µ)VT + λVN µ(VT + VN )
σ(v ) = .
` µ(VT + VN ) (λ + 2µ)VN + λVT

We next consider a very simple case:

√
E = 1, P = 0 (i.e., λ = 0, µ = 1/2), θ = 1, ` = 1, γ = 1 (i.e., γ0 = 1/ 2).

Therefore
3 1
A1 (uh , vh ) = (UT VT + UN VN ) + (UT VN + UN VT ) − UN VN .
4 4
Denoting by φN the (normal component) basis function of Vh (φN (A) = 1, φN (B) = 0), we get

P1 (vh ) = VN φN − VN and [P1 (uh )]+ = (1 − φN )(−UN )+ ,

so Z
1 1
[P1 (uh )]+ P1 (vh ) dΓ = − (−UN )+ VN .
ΓC γ 3
Besides
1
L(vh ) = − (g1 VT + g2 VN ).
2
The discrete problem (10) consists then of finding (UT , UN ) ∈ R2 such that:
3 1


 UT + UN = −g1 ,
2 2
1 1 2
UT − UN − (−UN )+ = −g2 .

 (29)
2 2 3
Clearly, a solution of (29) satisfies either UN ≥ 0 or UN < 0. We now show that for certain values
of g1 and g2 , system (29) admits an infinity of solutions or no solution.
• Suppose that g1 = 3g2 , then there exists an infinity of solutions (UT , UN ) = (−2g1 /3 − x/3, x)
for any x ≤ 0.
• Suppose that g1 > 3g2 , then (29) admits no solution.

References
[1] R. A. Adams, Sobolev spaces, Academic Press, New York-London, 1975. Pure and Applied
Mathematics, Vol. 65.

[2] P. Alart and A. Curnier, A generalized Newton method for contact problems with fric-
tion, J. de Mécanique Théorique et Appliquée, 7 (1988), pp. 67–82.

[3] D. Arnold, An interior penalty finite element method with discontinuous elements, SIAM
J. Numer. Anal., 19 (1982), pp. 742–760.

[4] R. Becker, P. Hansbo, and R. Stenberg, A finite element method for domain decompo-
sition with non-matching grids, M2AN Math. Model. Numer. Anal., 37 (2003), pp. 209–225.

23
 [5] Z. Belhachmi and F. Ben Belgacem, Quadratic finite element approximation of the
Signorini problem, Math. Comp., 72 (2003), pp. 83–104.

[6] F. Ben Belgacem and Y. Renard, Hybrid finite element methods for the Signorini prob-
lem, Math. Comp., 72 (2003), pp. 1117–1145.

[7] S.-C. Brenner and L.-R. Scott, The Mathematical Theory of Finite Element Methods,
vol. 15 of Texts in Applied Mathematics, Springer-Verlag, New York, 2007.

[8] H. Brezis, Équations et inéquations non linéaires dans les espaces vectoriels en dualité,
Ann. Inst. Fourier (Grenoble), 18 (1968), pp. 115–175.

[9] F. Chouly and P. Hild, A Nitsche-based method for unilateral contact problems: numerical
analysis, SIAM J. Numer. Anal., 51 (2013), pp. 1295–1307.

[10] F. Chouly and P. Hild, On convergence of the penalty method for unilateral contact
problems, Appl. Numer. Math., 65 (2013), pp. 27–40.

[11] P.-G. Ciarlet, Handbook of Numerical Analysis (eds. P.G. Ciarlet and J.L. Lions), vol. II,
North Holland, 1991, ch. 1. “The finite element method for elliptic problems”, pp. 17–352.

[12] P. Coorevits, P. Hild, K. Lhalouani, and T. Sassi, Mixed finite element methods for
unilateral problems: convergence analysis and numerical studies, Math. Comp., 71 (2002),
pp. 1–25.

[13] T. Dupont and R. Scott, Polynomial approximation of functions in Sobolev spaces, Math.
Comp., 34 (1980), pp. 441–463.

[14] A. Ern and J.-L. Guermond, Theory and practice of finite elements, vol. 159 of Applied
Mathematical Sciences, Springer-Verlag, New York, 2004.

[15] G. Fichera, Problemi elastostatici con vincoli unilaterali: Il problema di Signorini con
ambigue condizioni al contorno, Atti Accad. Naz. Lincei Mem. Cl. Sci. Fis. Mat. Natur. Sez.
I (8), 7 (1963/1964), pp. 91–140.

[16] A. Fritz, S. Hüeber, and B. I. Wohlmuth, A comparison of mortar and Nitsche tech-
niques for linear elasticity, Universität Stuttgart, IANS Preprint, (2003).

[17] , A comparison of mortar and Nitsche techniques for linear elasticity, Calcolo, 41 (2004),
pp. 115–137.

[18] R. Glowinski and P. Le Tallec, Augmented Lagrangian and operator-splitting meth-

ods in nonlinear mechanics, vol. 9 of SIAM Studies in Applied Mathematics, Society for
Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1989.

[19] W. Han and M. Sofonea, Quasistatic contact problems in viscoelasticity and viscoplastic-
ity, vol. 30 of AMS/IP Studies in Advanced Mathematics, American Mathematical Society,
Providence, RI, 2002.

[20] A. Hansbo and P. Hansbo, A finite element method for the simulation of strong and
weak discontinuities in solid mechanics, Comput. Methods Appl. Mech. Engrg., 193 (2004),
pp. 3523–3540.

24
[21] J. Haslinger, I. Hlaváček, and J. Nečas, Handbook of Numerical Analysis (eds. P.G.
Ciarlet and J.L. Lions), vol. IV, North Holland, 1996, ch. 2. “Numerical methods for unilat-
eral problems in solid mechanics”, pp. 313–385.

[22] P. Heintz and P. Hansbo, Stabilized Lagrange multiplier methods for bilateral elastic
contact with friction, Comput. Methods Appl. Mech. Engrg., 195 (2006), pp. 4323–4333.

[23] P. Hild, Numerical implementation of two nonconforming finite element methods for uni-
lateral contact, Comput. Methods Appl. Mech. Engrg., 184 (2000), pp. 99–123.

[24] P. Hild and P. Laborde, Quadratic finite element methods for unilateral contact problems,
Appl. Numer. Math., 41 (2002), pp. 401–421.

[25] P. Hild and Y. Renard, A stabilized Lagrange multiplier method for the finite element
approximation of contact problems in elastostatics, Numer. Math., 115 (2010), pp. 101–129.

[26] , An improved a priori error analysis for finite element approximations of Signorini’s
problem, SIAM J. Numer. Anal., 50 (2012), pp. 2400–2419.

[27] S. Hüeber and B. I. Wohlmuth, An optimal a priori error estimate for nonlinear multi-
body contact problems, SIAM J. Numer. Anal., 43 (2005), pp. 156–173.

[28] N. Kikuchi and J. T. Oden, Contact problems in elasticity: a study of variational inequal-
ities and finite element methods, vol. 8 of SIAM Studies in Applied Mathematics, Society for
Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 1988.

[29] N. Kikuchi and Y. J. Song, Penalty-finite-element approximation of a class of unilateral

problems in linear elasticity, Quart. Appl. Math., 39 (1981), pp. 1–22.

[30] P. Laborde and Y. Renard, Fixed point strategies for elastostatic frictional contact prob-
lems, Math. Methods Appl. Sci., 31 (2008), pp. 415–441.

[31] T. A. Laursen, Computational contact and impact mechanics, Springer-Verlag, Berlin,

2002.

[32] J.-L. Lions, Quelques méthodes de résolution des problèmes aux limites non linéaires,
Dunod, Paris, 1969.

[33] M. Moussaoui and K. Khodja, Régularité des solutions d’un problème mêlé dirichlet-
signorini dans un domaine polygonal plan, Commun. Part. Diff. Eq., 17 (1992), pp. 805–826.

[34] J. Nitsche, Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung
von Teilräumen, die keinen Randbedingungen unterworfen sind, Abhandlungen aus dem
Mathematischen Seminar der Universität Hamburg, 36 (1971), pp. 9–15.

[35] A. Popp, B. Wohlmuth, M. Gee, and W. Wall, Dual quadratic mortar finite element
methods for 3D finite deformation contact, SIAM J. Sci. Comput., 34 (2012), pp. B421–B446.

[36] Y. Renard, Generalized Newton’s methods for the approximation and resolution of frictional
contact problems in elasticity, Comp. Meth. Appl. Mech. Engrg., 256 (2013), pp. 38–55.

25
[37] R. Stenberg, On some techniques for approximating boundary conditions in the finite ele-
ment method, Journal of Computational and Applied Mathematics, 63 (1995), pp. 139–148.
International Symposium on Mathematical Modelling and Computational Methods Mod-
elling 94 (Prague, 1994).

[38] V. Thomée, Galerkin finite element methods for parabolic problems, vol. 25 of Springer
Series in Computational Mathematics, Springer-Verlag, Berlin, 1997.

[39] B. Wohlmuth, Variationally consistent discretization schemes and numerical algorithms

for contact problems, Acta Numerica, (2011), pp. 569–734.

[40] B. I. Wohlmuth, A. Popp, M. W. Gee, and W. A. Wall, An abstract framework for a

priori estimates for contact problems in 3D with quadratic finite elements, Comput. Mech.,
49 (2012), pp. 735–747.

[41] P. Wriggers, Computational Contact Mechanics, Wiley, 2002.

[42] P. Wriggers and G. Zavarise, A formulation for frictionless contact problems using a
weak form introduced by Nitsche, Comput. Mech., 41 (2008), pp. 407–420.

26