Weimin Han
email: han@uiowa.edu

Mircea Sofonea
email: sofonea@univ-perp.fr

David Danan
email: david.danan@univ-perp.fr

Numerical analysis of stationary variational-hemivariational inequalities

Keywords: Variational-hemivariational inequality, Clarke subdifferential, Galerkin approximation, finite element method, convergence, error estimation, Contact Mechanics AMS subject classifications. 65N30, 65N15, 74M10, 74M15

come

1. Introduction. Variational-hemivariational inequalities are inequality problems where both convex and nonconvex functions are involved. They were introduced in the pioneering work [START_REF] Panagiotopoulos | Nonconvex problems of semipermeable media and related topics[END_REF], and were further studied in [START_REF] Naniewicz | Mathematical Theory of Hemivariational Inequalities and Applications[END_REF][START_REF] Panagiotopoulos | Hemivariational Inequalities, Applications in Mechanics and Engineering[END_REF]. Interest in variationalhemivariational inequalities arises in the study of various problems in mechanics and engineering applications. Their study requires arguments of Convex Analysis, including properties of the subdifferential of a convex function, and arguments of Nonsmooth Analysis, including properties of the subdifferential in the sense of Clarke, defined for locally Lipschitz functions which may be nonconvex.

Recently, a new variational-hemivariational inequality is studied in [START_REF] Han | A class of variational-hemivariational inequalities with applications to frictional contact problems[END_REF]. The inequality involves two nonlinear operators and two nondifferentiable functionals, of which at least one is convex. Solution existence, uniqueness and data continuous dependence are shown. Moreover, the finite element method is studied for solving the inequality problem. For the first time in the literature, an optimal order error estimate is derived for the linear element solution of a hemivariational inequality under appropriate solution regularity assumptions. A more general variational-hemivariational inequality is analyzed in [START_REF] Migórski | A class of variational-hemivariational inequalities in reflexive Banach spaces[END_REF]. Solution existence and uniqueness are proved, together with a result on the continuous dependence of the solution on the data. The variational-hemivariational inequalities studied in [START_REF] Han | A class of variational-hemivariational inequalities with applications to frictional contact problems[END_REF][START_REF] Migórski | A class of variational-hemivariational inequalities in reflexive Banach spaces[END_REF] are motivated by applications in Contact Mechanics.

The purpose of this paper is to consider numerical approximations of stationary variational-hemivariational inequality problems by the finite element method, substantially extending the relevant result found in [START_REF] Han | A class of variational-hemivariational inequalities with applications to frictional contact problems[END_REF]. For the general variationalhemivariational inequality studied in [START_REF] Migórski | A class of variational-hemivariational inequalities in reflexive Banach spaces[END_REF], we show the convergence of the numerical solution. Then, for some particular variational-hemivariational inequalities, we also derive error estimates, which are of optimal order for the linear elements. We provide numerical examples to illustrate the performance of the numerical method, including numerical convergence orders. One of the numerical examples also serves the purpose of showing the transition from variational inequalities to hemivariational inequalities.

The rest of the paper is organized as follows. In Section 2 we review some preliminary material needed later on in the study of variational-hemivariational inequalities, and recall an existence and uniqueness result for the general variationalhemivariational inequality from [START_REF] Migórski | A class of variational-hemivariational inequalities in reflexive Banach spaces[END_REF]. In Section 3 we introduce numerical methods for solving various variational-hemivariational inequalities, prove convergence and derive error estimates wherever feasible. In Section 4 we introduce several contact problems, in which the material's behavior is modeled with a nonlinear elastic constitutive law and the contact conditions are in a subdifferential form. We list the assumptions on the data and use the abstract result in [START_REF] Migórski | A class of variational-hemivariational inequalities in reflexive Banach spaces[END_REF] to prove the unique weak solvability to the problem. In addition, we apply our results in Section 3 in the numerical analysis of the contact model. Finally, in Section 5 we present numerical simulations which represent an evidence of our error estimates and convergence results. The simulations illustrate the transition from the variational to the hemivariational case and give rise to interesting mechanical interpretations.

Preliminaries.

All linear spaces in this paper are real. For a normed space X we denote by • X its norm, by 0 X its zero element, by X * its topological dual, and by •, • X * ×X the duality pairing of X and X * . When no confusion may arise, we simply write •, • instead of •, • X * ×X . Weak convergence is indicated by the symbol . For two normed spaces X and Y , L(X, Y ) denotes the space of all linear continuous operators from X to Y .

We recall that an operator A : X → X * is pseudomonotone if it is bounded and u n u in X together with lim sup Au n , u n -u X * ×X ≤ 0 imply

Au, u -v X * ×X ≤ lim inf Au n , u n -v X * ×X ∀ v ∈ X.
A function ϕ : K ⊂ X → R is lower semicontinuous (l.s.c.), if for any sequence {x n } ⊂ K and any x ∈ K, x n → x in X implies ϕ(x) ≤ lim inf ϕ(x n ). For a convex function ϕ, the set

∂ϕ(x) := {x * ∈ X * | ϕ(v) -ϕ(x) ≥ x * , v -x X * ×X ∀ v ∈ X} is called the subdifferential of ϕ at x ∈ X. If ∂ϕ(x) is non-empty, an element x * ∈ ∂ϕ(x) is called a subgradient of ϕ at x.
Let ψ : X → R be a locally Lipschitz function. The generalized (Clarke) directional derivative of ψ at x ∈ X in the direction v ∈ X is defined by

ψ 0 (x; v) := lim sup y→x, λ↓0 ψ(y + λv) -ψ(y) λ .
The generalized subdifferential of ψ at x is a subset of the dual space X * given by

∂ψ(x) := { ζ ∈ X * | ψ 0 (x; v) ≥ ζ, v X * ×X ∀ v ∈ X }.
We have the formula ( [START_REF] Clarke | Optimization and Nonsmooth Analysis[END_REF]).

(2.1)

ψ 0 (x; v) = max { ζ, v | ζ ∈ ∂ψ(x)} .
More details on the properties of the subdifferential mappings, both in the convex and Clarke sense, can be found in the books [START_REF] Clarke | Optimization and Nonsmooth Analysis[END_REF][START_REF] Denkowski | An Introduction to Nonlinear Analysis: Theory[END_REF][START_REF] Denkowski | An Introduction to Nonlinear Analysis: Applications[END_REF][START_REF] Migórski | Nonlinear Inclusions and Hemivariational Inequalities. Models and Analysis of Contact Problems[END_REF][START_REF] Naniewicz | Mathematical Theory of Hemivariational Inequalities and Applications[END_REF][START_REF] Panagiotopoulos | Inequality Problems in Mechanics and Applications[END_REF][START_REF] Sofonea | Mathematical Models in Contact Mechanics[END_REF].

We turn now to the study of variational-hemivariational inequalities. Let X, X ϕ , X j be normed spaces and K ⊂ X. Given operators A : X → X * , γ ϕ : X → X ϕ , γ j : X → X j , and functionals ϕ : K × K → R, j : X → R, we consider the following problem.

Problem (P) Find an element u ∈ K such that (2.2)

Au, v -u + ϕ(γ ϕ u, γ ϕ v) -ϕ(γ ϕ u, γ ϕ u) + j 0 (γ j u; γ j v -γ j u) ≥ f, v -u ∀ v ∈ K.
For the study of Problem (P), we introduce the following assumptions on the data.

(A 1 ) X is a reflexive Banach space, and K is a closed and convex subset of X with 0 X ∈ K.

(A 2 ) X ϕ is a Banach space and γ ϕ ∈ L(X, X ϕ ): for a constant c ϕ > 0,

(2.3) γ ϕ v Xϕ ≤ c ϕ v X ∀ v ∈ X.
(A 3 ) X j is a Banach space and γ j ∈ L(X, X j ): for a constant c j > 0,

(2.4)

γ j v Xj ≤ c j v X ∀ v ∈ X.
(A 4 ) A : X → X * is pseudomonotone and there exists a constant m A > 0 such that

(2.5) Av 1 -Av 2 , v 1 -v 2 ≥ m A v 1 -v 2 2 X ∀ v 1 , v 2 ∈ X. (A 5 ) ϕ : K ϕ × K ϕ → R is such that ϕ(z, •) : K ϕ → R
is convex and l.s.c. on K ϕ for all z ∈ K ϕ , and there exists a constant α ϕ ≥ 0 such that

ϕ(z 1 , z 4 ) -ϕ(z 1 , z 3 ) + ϕ(z 2 , z 3 ) -ϕ(z 2 , z 4 ) (2.6) ≤ α ϕ z 1 -z 2 Xϕ z 3 -z 4 Xϕ ∀ z 1 , z 2 , z 3 , z 4 ∈ K ϕ .
Here, K ϕ := γ ϕ (K) is a non-empty, closed and convex set of X ϕ .

(A 6 ) j : X j → R is locally Lipschitz, and there are constants c 0 , c 1 , α j ≥ 0 such that

∂j(z) X * j ≤ c 0 + c 1 z Xj ∀ z ∈ X j , (2.7) j 0 (z 1 ; z 2 -z 1 ) + j 0 (z 2 ; z 1 -z 2 ) ≤ α j z 1 -z 2 2 Xj ∀ z 1 , z 2 ∈ X j . (2.8) (A 7 ) (2.9) α ϕ c 2 ϕ + α j c 2 j < m A . (A 8 ) (2.10) f ∈ X * .
Note that in the statement of Problem (P) the function ϕ(u, •) is assumed to be convex whereas the function j is locally Lipschitz and is, in general, nonconvex. For this reason, we refer to the inequality (2.2) as a variational-hemivariational inequality. Moreover, note that Problem (P) contains as particular cases, various problems considered in the literature. We comment that the spaces X ϕ and X j are introduced to facilitate error analysis of numerical solutions of Problem (P) in later sections. For applications in contact mechanics, the functionals ϕ(•, •) and j(•) are integrals over the contact boundary Γ 3 . In such a situation, X ϕ and X j can be chosen to be L 2 (Γ 3 ) d and/or L 2 (Γ 3 ). The assumption (2.9) is a smallness assumption, which is a variant of a similar condition in [START_REF] Migórski | A class of variational-hemivariational inequalities in reflexive Banach spaces[END_REF], due to the use of the spaces X ϕ and X j in the problem setting. By slightly modifying the proof in [START_REF] Migórski | A class of variational-hemivariational inequalities in reflexive Banach spaces[END_REF], we have the following existence and uniqueness result. Note that for a locally Lipschitz function j :

X j → R, the inequality (2.8) is equivalent to (2.11) ∂j(z 1 ) -∂j(z 2 ), z 1 -z 2 X * j ×Xj ≥ -α j z 1 -z 2 2 Xj ∀ z 1 , z 2 ∈ X j ,
known as the relaxed monotonicity condition (cf. [START_REF] Migórski | Nonlinear Inclusions and Hemivariational Inequalities. Models and Analysis of Contact Problems[END_REF] and the references therein).

In addition, if j : X j → R is a convex function, then (2.8) or (2.11), equivalently, are satisfied with α j = 0, due to the monotonicity of the (convex) subdifferential. Since 0 X ∈ K, we derive the following relations from (2.5), (2.7) and (2.11), for some constant c:

Av, v X * ×X ≥ m A v 2 X -c v X ∀ v ∈ X, (2.12) ∂j(z), z X * j ×Xj ≥ -α j z 2 Xj -c 0 z Xj ∀ z ∈ X j . (2.13)
3. Numerical approximations. In this section, we consider numerical schemes for solving Problem (P). We keep assumptions (A 1 )-(A 8 ) so that Problem (P) has a unique solution u ∈ K.

Let X h ⊂ X be a finite dimensional subspace with h > 0 denoting a spatial discretization parameter. Let K h = X h ∩ K. Then K h is a closed and convex subset of X h and 0 X ∈ K h . We consider the following Galerkin approximation of Problem (P).

Problem (P h ) Find an element

u h ∈ K h such that Au h , v h -u h + ϕ(γ ϕ u h , γ ϕ v h ) -ϕ(γ ϕ u h , γ ϕ u h ) + j 0 (γ j u h ; γ j v h -γ j u h ) (3.1) ≥ f, v h -u h ∀ v h ∈ K h .
The arguments of the proof of Theorem 2.1 can be applied in the setting of the finite dimensional space X h , and we know that under assumptions (A 1 )-(A 8 ), Problem (P h ) has a unique solution

u h ∈ K h . Proposition 3.1. For some constant M > 0, u h X ≤ M ∀ h > 0. Proof. We let v h = 0 X in (3.1) to get (3.2) Au h , u h ≤ ϕ(γ ϕ u h , 0) -ϕ(γ ϕ u h , γ ϕ u h ) + j 0 (γ j u h ; -γ j u h ) + f, u h .
For any z ∈ K ϕ , take z 1 = z 3 = z and z 2 = z 4 = 0 X in (2.6),

(3.3) ϕ(z, 0 X ) -ϕ(z, z) ≤ α ϕ z 2 Xϕ -ϕ(0 X , z) + ϕ(0 X , 0 X ).
Use the lower bound ([1, p. 433])

ϕ(0 X , z) ≥ c 3 + c 4 z Xϕ ∀ z ∈ X ϕ
for some constants c 3 and c 4 , not necessarily positive. Then from (3.3), we have

(3.4) ϕ(z, 0 X ) -ϕ(z, z) ≤ α ϕ z 2 Xϕ + c z Xϕ + 1 .
Write z 1 = z and take z 2 = 0 X in (2.8), j 0 (z; -z) ≤ α j z 2 Xj -j 0 (0 X ; z).

Further, use (2.7) to get

(3.5) j 0 (z; -z) ≤ α j z 2 Xj + c 0 z Xj .
Use (2.12), (3.4), (3.5), (2.3) and (2.4) in (3.2) to obtain

m A -α ϕ c 2 ϕ -α j c 2 j u h 2 X ≤ c u h X + 1 .
Since m A -α ϕ c 2 ϕ -α j c 2 j > 0, we deduce from the above inequality that u h X is uniformly bounded with respect to h.

The uniform boundedness of the numerical solutions will be useful in error analysis of the numerical method.

The focus of this section is error analysis for the numerical solution defined by Problem (P h ). We will assume, in addition, that A : X → X * is Lipschitz continuous, i.e. for some constant L A > 0,

(3.6) Au -Av X * ≤ L A u -v X ∀ u, v ∈ X.
Note that under conditions (3.6) and (2.5), the operator A is pseudomonotone ([25, Proposition 27.6]). We also assume that for fixed z ∈ K ϕ , ϕ(z, •) is Lipschitz continuous: for some function L ϕ :

K ϕ → R + , (3.7) |ϕ(z, z 1 ) -ϕ(z, z 2 )| ≤ L ϕ (z) z 1 -z 2 Xϕ ∀ z 1 , z 2 ∈ K ϕ .
Finally, we assume {K h } approximates K in the following sense:

(3.8) ∀ v ∈ K, ∃ v h ∈ K h such that v h → v in X as h → 0.
All the additional assumptions are valid for a wide variety of contact conditions, cf. Section 4.

3.1. Convergence. We begin with an application of (2.5) with v 1 = u and v 2 = u h to obtain, for any

v h ∈ K h , m A u -u h 2 X ≤ Au -Au h , u -v h + Au, v h -u (3.9) + Au, u -u h + Au h , u h -v h . From (2.2) with v = u h ∈ K, Au, u -u h ≤ ϕ(γ ϕ u, γ ϕ u h ) -ϕ(γ ϕ u, γ ϕ u) (3.10) + j 0 (γ j u; γ j u h -γ j u) -f, u h -u . From (3.1), Au h , u h -v h ≤ ϕ(γ ϕ u h , γ ϕ v h ) -ϕ(γ ϕ u h , γ ϕ u h ) (3.11) + j 0 (γ j u h ; γ j v h -γ j u h ) -f, v h -u h .
Using (3.10) and (3.11) in (3.9), we have

(3.12) m A u -u h 2 X ≤ Au -Au h , u -v h + R(v h ) + I ϕ (v h ) + I j (v h ),
where

R(v h ) := Au, v h -u + ϕ(γ ϕ u, γ ϕ v h ) -ϕ(γ ϕ u, γ ϕ u) (3.13) + j 0 (γ j u; γ j v h -γ j u) -f, v h -u , I ϕ (v h ) := ϕ(γ ϕ u, γ ϕ u h ) + ϕ(γ ϕ u h , γ ϕ v h ) -ϕ(γ ϕ u, γ ϕ v h ) -ϕ(γ ϕ u h , γ ϕ u h ), (3.14) I j (v h ) := j 0 (γ j u; γ j u h -γ j u) + j 0 (γ j u h ; γ j v h -γ j u h ) -j 0 (γ j u; γ j v h -γ j u). (3.15)
Notice that for ε > 0 arbitrarily small, there is a constant c depending on ε such that

Au -Au h , u -v h ≤ L A u -u h X u -v h X ≤ ε u -u h 2 X + c u -v h 2 X .
We further deduce from (3.12) that

(3.16) (m A -ε) u -u h 2 X ≤ c u -v h 2 X + R(v h ) + I ϕ (v h ) + I j (v h ).
Let us bound the terms I ϕ (v h ) and I j (v h ). By (2.6), we have

I ϕ (v h ) ≤ α ϕ γ ϕ u -γ ϕ u h Xϕ γ ϕ u h -γ ϕ v h Xϕ ≤ α ϕ c 2 ϕ u -u h 2 X + u -u h X u -v h X .
Thus, (3.17)

I ϕ (v h ) ≤ α ϕ c 2 ϕ + ε u -u h 2 X + c u -v h 2 X
for another constant c depending on ε > 0. We will apply the subadditivity of the generalized directional derivative:

j 0 (z; z 1 + z 2 ) ≤ j 0 (z; z 1 ) + j 0 (z; z 2 ) ∀ z, z 1 , z 2 ∈ X j . Using j 0 (γ j u h ; γ j v h -γ j u h ) ≤ j 0 (γ j u h ; γ j u -γ j u h ) + j 0 (γ j u h ; γ j v h -γ j u),
we have

I j (v h ) ≤ j 0 (γ j u; γ j u h -γ j u) + j 0 (γ j u h ; γ j u -γ j u h ) (3.18) + j 0 (γ j u h ; γ j v h -γ j u) -j 0 (γ j u; γ j v h -γ j u) . By (2.8), j 0 (γ j u; γ j u h -γ j u) + j 0 (γ j u h ; γ j u -γ j u h ) ≤ α j γ j u -γ j u h 2 Xj ; by (2.7), j 0 (γ j u h ; γ j v h -γ j u) ≤ c 0 + c 1 γ j u h Xj γ j u -γ j v h Xj , j 0 (γ j u; γ j v h -γ j u) ≤ c 0 + c 1 γ j u Xj γ j u -γ j v h Xj .
Thus, from (3.18), we have

(3.19) I j (v h ) ≤ α j γ j u -γ j u h 2 Xj + c γ j u -γ j v h Xj
for some constant c > 0 independent of h, where we used the fact that γ j u h Xj is uniformly bounded (cf. Proposition 3.1). Combine (3.16), (3.17), and (3.19),

m A -α ϕ c 2 ϕ -α j c 2 j -2 ε u -u h 2 X ≤ c u -v h 2 X + c γ j u -γ j v h Xj + R(v h ). Since α ϕ c 2 ϕ + α j c 2 j < m A , we can choose ε = (m A -α ϕ c 2 ϕ -α j c 2 j )/4 > 0 and get (3.20) u -u h 2 X ≤ c u -v h 2 X + γ j u -γ j v h Xj + R(v h ) ∀ v h ∈ K h .
The residual term (3.13) can be bounded as follows:

R(v h ) ≤ Au X * + L ϕ (γ ϕ u) c ϕ + c 0 + c 1 γ j u Xj c j + f X * u -v h X .
Thus, from (3.20), we have two constants c 1 , c 2 > 0, depending on the data of the problem and the solution, such that

(3.21) u -u h 2 X ≤ c 1 u -v h 2 X + c 2 u -v h X ∀ v h ∈ K h .
By (3.8), we can choose a sequence {v h } ⊂ K h that converges to u. From (3.21), we then conclude the convergence

u -u h X → 0 as h → 0.
Note that while (3.21) implies the convergence of the numerical method, it does not lead to optimal convergence order for error estimation.

3.2.

Error estimation for the particular case K = X. In the special case K = X, we have K h = X h , and the original problem (2.2) and its approximation (3.1) become

Au, v -u + ϕ(γ ϕ u, γ ϕ v) -ϕ(γ ϕ u, γ ϕ u) (3.22) + j 0 (γ j u; γ j v -γ j u) ≥ f, v -u ∀ v ∈ X and Au h , v h -u h + ϕ(γ ϕ u h , γ ϕ v h ) -ϕ(γ ϕ u h , γ ϕ u h ) (3.23) + j 0 (γ j u h ; γ j v h -γ j u h ) ≥ f, v h -u h ∀ v h ∈ X h . We replace v by 2u -v in (3.22), Au, u -v + ϕ(γ ϕ u, 2γ ϕ u -γ ϕ v) -ϕ(γ ϕ u, γ ϕ u) + j 0 (γ j u; γ j u -γ j v) ≥ f, u -v ∀ v ∈ X. Thus, Au, v h -u ≤ ϕ(γ ϕ u, 2γ ϕ u -γ ϕ v h ) -ϕ(γ ϕ u, γ ϕ u) + j 0 (γ j u; γ j u -γ j v h ) -f, u -v h ∀ v h ∈ X h .
Using this inequality in (3.12), we have

(3.24) m A u -u h 2 X ≤ Au -Au h , u -v h + Ĩϕ (v h ) + Ĩj (v h ),
where (3.14) and it is bounded in (3.17). We use (3.7) to find

Ĩϕ (v h ) := I ϕ (v h ) + ϕ(γ ϕ u, 2γ ϕ u -γ ϕ v h ) + ϕ(γ ϕ u, γ ϕ v h ) -2 ϕ(γ ϕ u, γ ϕ u), (3.25) Ĩj (v h ) := j 0 (γ j u; γ j u h -γ j u) + j 0 (γ j u; γ j u -γ j v h ) + j 0 (γ j u h ; γ j v h -γ j u h ). (3.26) Recall that I ϕ (v h ) is defined in
ϕ(γ ϕ u, 2γ ϕ u -γ ϕ v h ) + ϕ(γ ϕ u, γ ϕ v h ) -2 ϕ(γ ϕ u, γ ϕ u) ≤ 2 L ϕ (u) γ ϕ u -γ ϕ v h Xϕ .
This implies, combined with (3.17), that

(3.27) Ĩϕ (v h ) ≤ α ϕ c 2 ϕ + ε u -u h 2 X + c u -v h 2 X + 2 L ϕ (u) γ ϕ u -γ ϕ v h Xϕ .
Applying the inequality

j 0 (γ j u h ; γ j v h -γ j u h ) ≤ j 0 (γ j u h ; γ j u -γ j u h ) + j 0 (γ j u h ; γ j v h -γ j u),
we have

Ĩj (v h ) ≤ j 0 (γ j u; γ j u h -γ j u) + j 0 (γ j u h ; γ j u -γ j u h ) + j 0 (γ j u; γ j u -γ j v h ) + j 0 (γ j u h ; γ j v h -γ j u) .
This is similar to (3.18) and we have the following analogue of (3.19):

(3.28) Ĩj (v h ) ≤ α j c 2 j u -u h 2 X + c γ j u -γ j v h Xj .
Combining (3.24), (3.27) and (3.28), we get

m A u -u h 2 X ≤ Au -Au h X u -v h X + α ϕ c 2 ϕ + α j c 2 ϕ + ε u -u h 2 X + c u -v h 2 X + 2 L ϕ (u) γ ϕ u -γ ϕ v h Xϕ + c γ j u -γ j v h Xj .
Using the smallness assumption (2.9) and the Lipschitz condition (3.6), we deduce from the above inequality that

u -u h 2 X (3.29) ≤ c u -v h 2 X + γ ϕ u -γ ϕ v h Xϕ + γ j u -γ j v h Xj ∀ v h ∈ X h .
This is a basis for deriving error estimates.

4. Error analysis for contact problems. We illustrate applications of the framework developed in Section 3 on convergence and error estimation for numerical solutions of a number of static contact problems with elastic materials. Let Ω be the reference configuration of the elastic body, assumed to be an open, bounded, connected set in R d (d = 2, 3). The boundary Γ = ∂Ω is assumed Lipschitz continuous and is partitioned into three disjoint and measurable parts Γ 1 , Γ 2 and Γ 3 such that meas (Γ 1 ) > 0. The body is in equilibrium under the action of a total body force of density f 0 in Ω and a surface traction of density f 2 on Γ 2 , is fixed on Γ 1 , and is in contact on Γ 3 with a foundation. Different contact conditions lead to different contact problems, as discussed below.

We use S d for the space of second order symmetric tensors on R d . Also, "•" and • will represent the canonical inner product and the Euclidean norm on the spaces R d and S d . We denote by u : Ω → R d and σ : Ω → S d the displacement field and the stress field, respectively. In addition, we use ε(u) to denote the linearized strain tensor. Let ν be the unit outward normal vector, defined a.e. on Γ. For a vector field v, we use v ν := v • ν and v τ := v -v ν ν for the normal and tangential components of v on Γ. Similarly, for the stress field σ, its normal and tangential components on the boundary are defined as σ ν := (σν) • ν and σ τ := σν -σ ν ν, respectively. Then for the contact problems under consideration, we have the elastic constitutive law

(4.1) σ = Fε(u) in Ω, the equilibrium equation (4.2) Div σ + f 0 = 0 in Ω, the displacement boundary condition (4.3) u = 0 on Γ 1 ,
and the traction boundary condition

(4.4) σν = f 2 on Γ 2 .
The relations (4.1)-(4.4) will be supplemented by a set of boundary conditions on Γ 3 . Note that in (4.1)-(4.4) and sometimes below, we do not indicate explicitly the dependence of various functions on the spatial variable x ∈ Ω ∪ Γ. In (4.1), F : Ω × S d → S d is the elasticity operator and is assumed to have the following properties:

(4.5)                (a) there exists L F > 0 such that for all ε 1 , ε 2 ∈ S d , a.e. x ∈ Ω, F(x, ε 1 ) -F(x, ε 2 ) ≤ L F ε 1 -ε 2 ; (b) there exists m F > 0 such that for all ε 1 , ε 2 ∈ S d , a.e. x ∈ Ω, (F(x, ε 1 ) -F(x, ε 2 )) • (ε 1 -ε 2 ) ≥ m F ε 1 -ε 2 2 ; (c) F(•, ε) is measurable on Ω for all ε ∈ S d ; (d) F(x, 0) = 0 for a.e. x ∈ Ω.
To study the contact problems, we need some function spaces. For the stress and strain fields, we use the space Q = L 2 (Ω; S d ), which is a Hilbert space with the canonical inner product

(σ, τ ) Q := Ω σ ij (x) τ ij (x) dx, σ, τ ∈ Q;
the associated norm is denoted by • Q . The displacement fields will be sought in the space

V = v = (v i ) ∈ H 1 (Ω; R d ) | v = 0 a.e. on Γ 1
or its subset. Since meas (Γ 1 ) > 0, it is known that V is a Hilbert space with the inner product

(u, v) V := Ω ε(u) • ε(v) dx, u, v ∈ V
and the associated norm • V . For v ∈ H 1 (Ω; R d ) we use the same symbol v for the trace of v on Γ. By the Sobolev trace theorem we have

v L 2 (Γ3;R d ) ≤ γ v V ∀ v ∈ V, γ being the norm of the trace operator γ : V → L 2 (Γ 3 ; R d ).
We consider several choices of the boundary conditions on the contact boundary Γ 3 , leading to different contact problems which are examples of Problem (P) or its special cases. a) Frictional contact with normal compliance and unilateral constraint. The first set of contact boundary conditions is ( [START_REF] Migórski | A class of variational-hemivariational inequalities in reflexive Banach spaces[END_REF])

u ν ≤ g, σ ν + ξ ν ≤ 0, (u ν -g) (σ ν + ξ ν ) = 0, ξ ν ∈ ∂j ν (u ν ) on Γ 3 , (4.6) σ τ ≤ F b (u ν ), -σ τ = F b (u ν ) u τ u τ if u τ = 0 on Γ 3 . (4.7)
Here g > 0, ∂j ν is the Clarke subdifferential of a function j ν , and the friction bound F b is a positive function. The condition (4.6) models the contact with a foundation made of a rigid body covered by a layer made of elastic material, say asperities. The relation u ν ≤ g restricts the allowed penetration, where g represents the thickness of the elastic layer. When there is penetration, as long as the normal displacement does not reach the bound g, the contact is described with a multivalued normal compliance condition -σ ν = ξ ν ∈ ∂j ν (u ν ). Thus, the unknown ξ ν may be interpreted as the opposite of the normal stress on the contact surface. Examples of normal compliance contact conditions in the subdifferential form -σ ν ∈ ∂j ν (u ν ) can be found in [START_REF] Migórski | Nonlinear Inclusions and Hemivariational Inequalities. Models and Analysis of Contact Problems[END_REF]. Thus, the contact condition (4.6) represents a combination of the Signorini contact condition (which models the contact with a rigid foundation) and the normal compliance condition (which models the contact with a deformable foundation). Details on the normal compliance and Signorini contact conditions can be found in [START_REF] Han | Quasistatic Contact Problems in Viscoelasticity and Viscoplasticity[END_REF][START_REF] Migórski | Nonlinear Inclusions and Hemivariational Inequalities. Models and Analysis of Contact Problems[END_REF][START_REF] Shillor | Models and Analysis of Quasistatic Contact[END_REF][START_REF] Sofonea | Mathematical Models in Contact Mechanics[END_REF]. The contact is assumed to be frictional and is described with a version of Coulomb's law of dry friction, (4.7). The friction bound F b may depend on the normal displacement u ν , cf. [START_REF] Sofonea | Numerical analysis of history-dependent variational inequalities with applications to contact problems[END_REF] for explanation.

On the potential function j ν : Γ 3 × R → R, we assume (4.8)

                     (a) j ν (•, r) is measurable on Γ 3 for all r ∈ R and there exists e ∈ L 2 (Γ 3 ) such that j ν (•, e(•)) ∈ L 1 (Γ 3 ); (b) j ν (x, •) is locally Lipschitz on R for a.e. x ∈ Γ 3 ; (c) |∂j ν (x, r)| ≤ c 0 + c 1 |r| for a.e. x ∈ Γ 3 , for all r ∈ R with c 0 , c 1 ≥ 0; (d) j 0 ν (x, r 1 ; r 2 -r 1 ) + j 0 ν (x, r 2 ; r 1 -r 2 ) ≤ α jν |r 1 -r 2 | 2 for a.e. x ∈ Γ 3 , all r 1 , r 2 ∈ R with α jν ≥ 0.
On the penetration bound g : Γ 3 → R and the friction bound

F b : Γ 3 × R → R + , we assume (4.9) g ∈ L 2 (Γ 3 ), g(x) ≥ 0 a.e. on Γ 3 , (4.10)
       (a) there exists L F b > 0 such that |F b (x, r 1 ) -F b (x, r 2 )| ≤ L F b |r 1 -r 2 | ∀ r 1 , r 2 ∈ R, a.e. x ∈ Γ 3 ; (b) F b (•, r) is measurable on Γ 3 , for all r ∈ R; (c) F b (x, r) = 0 for r ≤ 0, F b (x, r) ≥ 0 for r ≥ 0, a.e. x ∈ Γ 3 .
On the densities of body forces and surface tractions we assume (4.11)

f 0 ∈ L 2 (Ω; R d ), f 2 ∈ L 2 (Γ 2 ; R d ). Define f ∈ V * by (4.12) f , v V * ×V = (f 0 , v) L 2 (Ω;R d ) + (f 2 , v) L 2 (Γ2;R d ) ∀ v ∈ V.
Corresponding to the constraint u ν ≤ g on Γ 3 in (4.6), we introduce a subset of the space V :

(4.13)

U := {v ∈ V | v ν ≤ g on Γ 3 } .
By a standard approach (cf. [START_REF] Han | Quasistatic Contact Problems in Viscoelasticity and Viscoplasticity[END_REF][START_REF] Migórski | Nonlinear Inclusions and Hemivariational Inequalities. Models and Analysis of Contact Problems[END_REF]), the following weak formulation of the first contact problem can be derived.

Problem (P 1 ). Find a displacement field u ∈ U such that

(F(ε(u)), ε(v -u)) Q + Γ3 F b (u ν ) ( v τ -u τ ) dΓ (4.14) + Γ3 j 0 ν (u ν ; v ν -u ν ) dΓ ≥ f , v -u V * ×V ∀ v ∈ U.
To apply the theory presented in the previous sections, we let

X = V , K = U , X ϕ = L 2 (Γ 3 ) d with γ ϕ the trace operator from V to X ϕ , X j = L 2 (Γ 3 ) with γ j v = v ν for v ∈ V . Then, α ϕ = L F b and α j = α jν . Let λ 1,V > 0 be the smallest eigenvalue of the eigenvalue problem u ∈ V, Ω ε(u)•ε(v) dx = λ Γ3 u•v dΓ ∀ v ∈ V,
and let λ 1ν,V > 0 be the smallest eigenvalue of the eigenvalue problem

u ∈ V, Ω ε(u)•ε(v) dx = λ Γ3 u ν v ν dΓ ∀ v ∈ V.
Then we may take

c ϕ = λ -1/2 1,V , c j = λ -1/2 1ν,V .
Applying Theorem 2.1, we know that under the assumptions (4.5), (4.8), (4.10), and

(4.15) L F b λ -1 1,V + α jν λ -1 1ν,V < m F ,
Problem (P 1 ) has a unique solution u ∈ U . We now consider the finite element method of solving Problem (P 1 ). For simplicity, assume Ω is a polygonal/polyhedral domain and express the three parts of the boundary, Γ k , 1 ≤ k ≤ 3, as unions of closed flat components with disjoint interiors:

Γ k = ∪ i k i=1 Γ k,i , 1 ≤ k ≤ 3.
Let {T h } be a regular family of partitions of Ω into triangles/tetrahedrons that are compatible with the partition of the boundary

∂Ω into Γ k,i , 1 ≤ i ≤ i k , 1 ≤ k ≤ 3,
in the sense that if the intersection of one side/face of an element with one set Γ k,i has a positive measure with respect to Γ k,i , then the side/face lies entirely in Γ k,i . Construct the linear element space corresponding to T h :

V h = v h ∈ C(Ω) d | v h | T ∈ P 1 (T ) d , T ∈ T h , v h = 0 on Γ 1 ,
and the related finite element subset U h = V h ∩ U . Assume g is a concave function. Then

U h = v h ∈ V h | v h ν ≤ g at node points on Γ 3 .
Note that 0 ∈ U h . Define the following numerical method for Problem (P 1 ).

Problem (P h 1 ). Find a displacement field u h ∈ U h such that (F(ε(u h )), ε(v h -u h )) Q + Γ3 F b (u h ν ) v h τ -u h τ dΓ (4.16) + Γ3 j 0 ν (u h ν ; v h ν -u h ν ) dΓ ≥ f , v h -u h V * ×V ∀ v h ∈ U h .
For an error analysis, we assume

(4.17) u ∈ H 2 (Ω) d , σν ∈ L 2 (Γ 3 ) d .
Note that for many application problems, σν ∈ L 2 (Γ 3 ) d follows from u ∈ H 2 (Ω) d ; e.g., this is the case where the material is linearly elastic with suitably smooth coefficients, or where the elasticity operator F depends on x smoothly. We apply (3.20) to derive an error estimate. For this purpose, we need to bound the residual term defined in (3.13):

R(v h ) = (F(ε(u)), ε(v h -u)) Q + Γ3 F b (u ν ) v h τ -u τ dΓ + Γ3 j 0 ν (u ν ; v h ν -u ν ) dΓ -f , v h -u V * ×V .
We follow the procedure found in [START_REF] Han | Quasistatic Contact Problems in Viscoelasticity and Viscoplasticity[END_REF]. Take v = u ± w with w in the subset Ũ of U defined by

Ũ := w ∈ C ∞ (Ω) d | w = 0 on Γ 1 ∪ Γ 3 ,
and derive from (4.14) that

(F(ε(u)), ε(w)) Q = f , w V * ×V ∀ w ∈ Ũ . Therefore, Div F(ε(u)) + f 0 = 0 in Ω, (4.18) σν = f 2 on Γ 2 . (4.19)
Then multiply (4.18) by vu with v ∈ U , integrate over Ω, and integrate by parts,

∂Ω σν•(v -u) dΓ - Ω F(ε(u))•ε(v -u) dx + Ω f 0 •(v -u) dx = 0, i.e., (4.20) Ω F(ε(u))•ε(v -u) dx = f , v -u V * ×V + Γ3 σν•(v -u) dΓ. Thus, R(v h ) = Γ3 σν•(v h -u) + F b (u ν ) v h τ -u τ + j 0 ν (u ν ; v h ν -u ν ) dΓ,
and then,

(4.21) R(v h ) ≤ c u -v h L 2 (Γ3) d .
Finally, from (3.20), we derive the inequality

(4.22) u -u h 2 V ≤ c u -v h 2 V + u -v h L 2 (Γ3) d ∀ v h ∈ U h .
Under additional solution regularity assumption (4.23)

u ν | Γ3,i ∈ H 2 (Γ 3,i ; R d ), 1 ≤ i ≤ i 3 ,
we have the optimal order error bound

(4.24) u -u h V ≤ c h.
We comment that similar results hold for the frictionless version of the model, i.e., where the friction condition (4.7) is replaced by

σ τ = 0 on Γ 3 .
Then the problem is to solve the inequality (4.14) without the term

Γ3 F b (u ν ) ( v τ -u τ ) dΓ.
The condition (4.15) reduces to

α jν λ -1 1ν,V < m F .
The inequality (4.22) and the error bound (4.24) still hold for the linear finite element solution.

b) Frictional contact with normal compliance. Instead of (4.6)-(4.7), the following contact boundary conditions were considered in [START_REF] Han | A class of variational-hemivariational inequalities with applications to frictional contact problems[END_REF]:

(4.25) -σ ν ∈ ∂j ν (u ν ), σ τ ≤ F b (u ν ), -σ τ = F b (u ν ) u τ u τ if u τ = 0 on Γ 3 ,
where the potential function j ν is assumed to satisfy (4.8), whereas the friction bound F b is assumed to satisfy (4.10). Note that in contrast with the contact boundary condition (4.6) which involves a unilateral constraint on the displacement field, the contact condition in (4.25) does not involve such a restriction. It can be viewed as a limiting case of the (4.6) when g → ∞. The weak formulation of the corresponding contact problem is the following.

Problem (P 2 ). Find a displacement field u ∈ V such that (F(ε(u)), ε(v -u)) Q + Γ3 F b (u ν ) ( v τ -u τ ) dΓ (4.26) + Γ3 j 0 ν (u ν ; v ν -u ν ) dΓ ≥ f , v -u V * ×V ∀ v ∈ V.
This problem and its numerical approximations were studied in [START_REF] Han | A class of variational-hemivariational inequalities with applications to frictional contact problems[END_REF], where the functional j ν was assumed to be Lipschitz continuous. We can apply the theory developed in Section 3 for error analysis of numerical solutions of Problem (P 2 ), with X ϕ = L 2 (Γ 3 ) d and X j = L 2 (Γ 3 ), without the need of assuming j ν to be Lipschitz continuous. -σ ν ∈ ∂j ν (u ν ), σ τ = 0 on Γ 3 .

These conditions can be obtained from (4.25) when the friction bound vanishes, i.e. when F b ≡ 0. Such a condition represents an idealization of the process, since even completely lubricated surfaces generate shear resistance to tangential motion. The weak formulation of the corresponding contact problem is the following.

Problem (P 3 ). Find a displacement field u ∈ V such that (4.28) (F(ε(u)), ε(v -u)) Q + Γ3 j 0 ν (u ν ; v ν -u ν ) dΓ ≥ f , v -u V * ×V ∀ v ∈ V.
Note that, in contrast with the inequality (4.26) which involves both a convex and a nonconvex function, the inequality (4.28) is governed by a nonconvex function only; it is an example of a pure hemivariational inequality. In applying the theory of Section 3 in the study of this inequality, we choose X j = L 2 (Γ 3 ) and there is no X ϕ . d) Frictionless contact with subdifferential boundary conditions and unilateral constraint. We turn now to a new model of frictional contact described with two subdiferential boundary conditions, associated with a unilateral constraint for the normal displacement field. The boundary conditions are formulated as follows:

σ ν = σ 1 ν + σ 2 ν + σ 3 ν , -σ 1 ν ∈ ∂ϕ ν (u ν ), -σ 2 ν ∈ ∂j ν (u ν ), (4.29) u ν ≤ g, σ 3 ν ≤ 0, σ 3 ν (u ν -g) = 0, on Γ 3 , σ τ = 0 on Γ 3 . (4.30)
Here g is a prescribed bound, ϕ ν and j ν are given functions, and ∂ϕ ν , ∂j ν denote the convex subdifferential of ϕ ν and the Clarke subdifferential of the j ν , respectively. A relevant example of the contact condition which can be cast in the general subdifferential framework (4.29) follows.

Example 4.1. We adopt assumptions a)-e) below in which equalities and inequalities hold on the contact surface Γ 3 .

a) The foundation is made of a rigid body covered by a layer of soft material, say asperities. Therefore, the penetration is restricted, i.e. (4.31) u ν ≤ g, where g > 0 represents the thickness of the soft layer. We consider the nonhomogeneous case, i.e., g is allowed to be a function of the spatial variable x ∈ Γ 3 .

b) The normal stress has an additive decomposition of the form (4.32)

σ ν = σ D ν + σ R ν ,
where the term σ D ν describes the reaction of the soft layer and σ R ν describes the reaction of the rigid body.

c) The part σ D ν has, in turn, an additive decomposition of the form (4.33)

σ D ν = σ 1 ν + σ 2 ν ,
where

(4.34) -σ 1 ν = q ν (u ν ), -σ 2 ν = p ν (u ν ).
Here q ν and p ν are continuous positive functions, vanishing for a negative argument. Moreover, q ν is an increasing function and p ν is allowed to be nonmonotone. Equalities (4.33), (4.34) show that the function σ D ν follows a normal compliance contact condition of the form

(4.35) -σ D ν = q ν (u ν ) + p ν (u ν ).
d) The part σ R ν satisfies the Signorini condition in a form with a gap function, i.e.

(4.36) σ R ν ≤ 0, σ R ν (u ν -g) = 0.
The contact conditions (4.34) can be expressed in terms of the subdifferential operators. Introduce functions ϕ ν : R → R and j ν : R → R by (4.37)

ϕ ν (r) = r 0 q ν (s) ds, j ν (r) = r 0 p ν (s) ds ∀ r ∈ R.
Then, as explained in [START_REF] Migórski | Nonlinear Inclusions and Hemivariational Inequalities. Models and Analysis of Contact Problems[END_REF][START_REF] Sofonea | Mathematical Models in Contact Mechanics[END_REF] we have ∂ϕ(r) = q ν (r), ∂ ν j(r) = p ν (r) for all r ∈ R. Thus, (4.34) leads to a subdifferential condition of the form

(4.38) -σ 1 ν ∈ ∂ϕ ν (u ν ), -σ 2 ν ∈ ∂j(u ν ).
Denote σ R ν = σ 3 ν . Then, gathering relations (4.32), (4.33), (4.38), (4.31) and (4.36) we obtain that u ν and σ ν satisfy the contact condition (4.29). We conclude that the contact model based on assumptions (4.31)-(4.36) can be cast in this abstract subdifferential setting, as claimed.

Remark 4.2. We can write the contact condition in Example 4.1 in an equivalent form, which shall be useful for the mechanical interpretation of the numerical results we present in Section 5. We use equalities (4.33) and (4.34) to see that σ D ν = -q ν (u ν ) -p ν (u ν ). We replace this equality in (4.32) to find that σ R ν = σ ν + q ν (u ν ) + p ν (u ν ). Then we substitute this equality in (4.36) and use condition (4.31) to obtain that

(4.39) u ν ≤ g, σ ν + q ν (u ν ) + p ν (u ν ) ≤ 0, (σ ν + q ν (u ν ) + p ν (u ν ))(u ν -g) = 0. Denote (4.40) k ν (r) = q ν (r) + p ν (r), r ∈ R.
Then, (4.39) becomes

(4.41) u ν ≤ g, σ ν + k ν (u ν ) ≤ 0, (σ ν + k ν (u ν ))(u ν -g) = 0.
This represents the normal compliance contact condition with unilateral constraint.

It was first introduced in [START_REF] Jarušek | On the solvability of dynamic elastic-visco-plastic contact problems[END_REF] under the assumption that k ν is a Lipschitz continuous increasing functions. Here, k ν is only assumed to be locally Lipschitz continuous and it may be non-monotone.

Note that various examples of contact conditions can be considered in the form (4.29), some of them involving multivalued functions. For this reason, we proceed our analysis by considering the contact conditions in the general form (4.29). There, the function j ν is assumed to satisfy (4.8) and ϕ ν : Γ 3 × R → R satisfies the following conditions:

(4.42)    (a) ϕ ν (•, r) is measurable on Γ 3 for all r ∈ R and there exists e ∈ L 2 (Γ 3 ) such that ϕ ν (•, e(•)) ∈ L 1 (Γ 3 ). (b) ϕ ν (x, •) is convex on R for a.e. x ∈ Γ 3 .
By a standard procedure, we derive the following variational formulation of the problem (4.1)-(4.4), (4.29), (4.30).

Problem (P 4 ). Find a displacement field u ∈ U such that Ω F(ε(u)) • ε(v -u) dx (4.43) + Γ3 ϕ ν (v ν ) dΓ - Γ3 ϕ ν (u ν ) dΓ + Γ3 j 0 ν (u ν ; v ν -u ν ) dΓ ≥ Ω f 0 • (v -u) dx + Γ2 f 2 • (v -u) dΓ ∀ v ∈ U.
The unique solvability of Problem (P 4 ) is given by the following existence and uniqueness result.

Theorem 4.3. Assume (4.5), (4.8), (4.9), (4.11), (4.42) and, in addition, assume the smallness condition

(4.44) α jν λ -1 1ν,V < m F .
Then Problem (P 4 ) has a unique solution u ∈ U . Proof. We apply Theorem 2.1 with X = V , K = U being the set defined by (4.13),

A : V → V * , Au, v = Ω Fε(u) • ε(v) dx for u, v ∈ V, ϕ : V × V → R, ϕ(u, v) = Γ3 ϕ ν (v ν ) dΓ for u, v ∈ V, j : V → R, j(v) = Γ3 j ν (v ν ) dΓ for v ∈ V, f ∈ V * , f , v = Ω f 0 • v dx + Γ3 f 2 • v dx for v ∈ V.
Observe that hypothesis (4.5) implies that the operator A satisfies condition (2.5) with α A = m A = m F . Moreover, condition (4.42) implies that ϕ satisfies (2.6) with α ϕ = 0. Next, hypothesis (4.8)(a) guaranties that the function j is well defined, and conditions (4.8)(b) and (c) imply that (2.7) holds. For more details in proving these properties we refer the reader to Theorem 3.47 in [START_REF] Migórski | Nonlinear Inclusions and Hemivariational Inequalities. Models and Analysis of Contact Problems[END_REF]. Condition (2.8) is a consequence of the relation

j 0 (u; v) ≤ Γ3 j 0 ν (x, u ν (x); v ν (x)) dΓ ∀ u, v ∈ V,
combined with the hypotheses (4.8)(d). Thus, we conclude that j satisfies condition (2.7) with α j = α jν γ 2 . Moreover, using (4.9), it is easy to see that the set (4.13) is a nonempty, closed, convex set in V . In addition, (4.12) implies (2.10). Finally, we see from above that assumption (4.44) implies the smallness condition (2.9). Theorem 4.3 is now a direct consequence of Theorem 2.1.

The results in Section 3 apply in the study the numerical approximation of the contact problems (P 2 )-(P 4 ) as was done for the contact problem (P 1 ) earlier in the section. In particular, under the solution regularity assumptions (4.17) and (4.23), we have the optimal order error bound (4.24).

Numerical simulations.

In this section we present numerical simulation results. We restrict ourselves to provide numerical simulations in the study of Problems (P 4 ) and (P 2 ).

The numerical algorithm we use is described in [START_REF] Barboteu | An analytical and numerical approach to a bilateral contact problem with nonmonotone friction[END_REF][START_REF] Barboteu | Analysis of a contact problem with normal compliance, finite penetration and nonmonotone slip dependent friction[END_REF][START_REF] Barboteu | Analysis of quasistatic viscoplastic contact problems with normal compliance[END_REF][START_REF] Barboteu | On the behavior of the solution of a viscoelastic contact problem[END_REF]. The numerical solution is based on an iterative procedure which leads to a sequence of convex programming problems. For each "convexification" iteration, the value of the normal compliance function is fixed to a given value depending on the normal displacement solution u ν found in the previous iteration. Then, the resulting nonsmooth convex iterative problems are solved by classical numerical methods. Furthermore, the frictional contact conditions are treated by using a numerical approach based on the combination of the penalized method and the augmented Lagrangian method. We consider additional fictitious nodes for the Lagrange multiplier in the initial mesh. The construction of these nodes depends on the contact element used for the geometrical discretization of the interface Γ 3 . In the examples presented below, the discretization is based on "node-to-rigid" contact element, which is composed by one node of Γ 3 and one Lagrange multiplier node. For more details on the discretization step and Computational Contact Mechanics, we refer to [START_REF] Khenous | On the discretization of contact problems in elastodynamics[END_REF][START_REF] Khenous | Hybrid discretization of the Signorini problem with Coulomb friction. Theoretical aspects and comparison of some numerical solvers[END_REF][START_REF] Laursen | Computational Contact and Impact Mechanics[END_REF][START_REF] Wriggers | Computational Contact Mechanics[END_REF].

Physical setting and values of parameters. Let Ω = (0, L) × (0, L) ⊂ R 2 with L > 0 and

Γ 1 = ({0} × [0, L]) ∪ ({L} × [0, L]), Γ 2 = [0, L] × {L}, Γ 3 = [0, L] × {0}.
The domain Ω represents the cross section of a three-dimensional elastic body subjected to the action of tractions in such a way that the plane stress hypothesis is valid. On Γ 1 = ({0} × [0, L]) ∪ ({L} × [0, L]), the body is clamped and, therefore, the displacement field vanishes there. Vertical tractions act on Γ 2 = [0, L] × {L}. No body forces are assumed to act on the body during the process. The body is in contact with an obstacle on Γ 3 = [0, L]×{0}. The contact conditions used correspond both to Problems (P 4 ) and (P 2 ) and will be described below. The material response is governed by a linear constitutive law defined by the elasticity tensor F given by

(Fτ ) ij = Eκ (1 + κ)(1 -2κ) (τ 11 + τ 22 )δ ij + E 1 + κ τ ij , 1 ≤ i, j ≤ 2, ∀ τ ∈ S 2 .
Here, E and κ are Young's modulus and Poisson's ratio of the material and δ ij denotes the Kronecker delta. For the computation below, we use the following data:

L = 1 m, E = 70 GP a, κ = 0.3, f 0 = (0, 0) GP a, f 2 = (0, -4) GP a•m on Γ 2 .
For the numerical simulation we use linear finite elements on uniform triangulations of the domain Ω. The contact boundary Γ 3 is divided in 1/h parts, h being the spatial discretization parameter. The numerical solutions presented below correspond to the case h = 1/64 where the spatial domain is discretized into 16384 elements for a total number of degrees of freedom equal to 16770.

Numerical simulations for Problem (P 4 ). The contact is frictionless; it follows a normal compliance condition as far as the penetration is less than the bound g and, when this bound is reached, it follows a unilateral constraint. The behavior of the foundation is an elastic-rigid one and corresponds to an obstacle made of a hard material covered by a layer composed of a soft material, say asperities, with thickness g, as depicted in Figure 5.1. For the numerical simulations we choose g = 0.02 m. In addition, we choose the contact condition (4.39) which represents a particular case of the contact condition (4.29), as explained in Example 4.1 and Remark 4.2. For simplicity, we consider only the homogeneous case, i.e., we skip the dependence on the function on the spatial variable x. Let (5.1) q ν (r) = αβr + , r ∈ R, If β ≥ 1, then k ν is a continuous increasing function and, therefore, ψ ν is a convex function. Moreover, it can be proved that u is a solution of the inequality (5.5) if and only if

u ∈ U, Ω F(ε(u)) • ε(v -u) dx + Γ3 [ψ ν (v ν ) -ψ ν (v ν )] dΓ (5.7) ≥ Ω f 0 • (v -u) dx + Γ3 f 2 • (v -u) dx ∀ v ∈ U.
Since ψ ν is a convex function, the inequality (5.7) is a purely variational inequality.

Next, if 0 ≤ β < 1, then k ν is not a monotone function and, therefore, ψ ν is not convex. Nevertheless, the function ψ ν is locally Lipschitz. It can be proved that u is a solution of the inequality (5.5) if and only if

u ∈ U, Ω F(ε(u)) • ε(v -u) dx + Γ3 ψ 0 ν (u ν ; v ν -u ν ) dΓ (5.8) ≥ Ω f 0 • (v -u) dx + Γ3 f 2 • (v -u) dx ∀ v ∈ U,
where ψ 0 ν denotes generalized derivative in the sense of Clarke. Since ψ is nonconvex, the inequality (5.8) is a purely hemivariational inequality.

We conclude that our contact model leads to a variational formulation (5.5) whose intrinsic nature depends on the value of the parameter β. For β ≥ 1, the function k ν is increasing and, therefore, using (4.39), (4.40) it follows that we have the case of a monotone normal compliance contact condition. In contrast, for 0 ≤ β < 1 the function k ν is not increasing and using (4.39), (4.40) we have the case of a socalled non monotone normal compliance contact condition. Our aim in what follows it to perform simulations in the two cases above and to provide the corresponding mechanical interpretations.

a) The monotone case. In Figure 5.2 we plot the deformed configuration as well as the interface forces on Γ 3 , for α = 40 while in Figure 5.3 we plot that for α = 10. In both cases we take β = 2 which guarantees that we are in the monotone case.

We note that all nodes in Figure 5.2 are in status of normal compliance, i.e. 0 ≤ u ν < 0.02. Moreover, the interface forces increase with respect to the penetration. This numerical result corresponds with the theoretic one. Indeed, in this case there is no contact with the rigid foundation, the normal stress reduces to its component provided by the normal compliance condition and, therefore -σ ν = k ν (u ν ). Since k ν is an increasing function we deduce that the magnitude of σ ν is increasing with the penetration.

In Figure 5.3 part of the nodes are in status of normal compliance (i.e. 0 ≤ u ν < 0.02) and part of them are in unilateral contact (i.e. u ν = 0.02). This situation arises since the stiffness of the deformable foundation considered in Figure 5.3 is α = 10 which is lower than the one used in Figure 5.2, where α = 40. As a result, there is a complete flattening of the asperities in the center of the contact boundary. Moreover, note that the interface forces increase with respect to the penetration and this agrees with the theory, since we are in the monotone case. In addition, we note that the interface forces are more important on the unilateral contact zone since the component σ 3 ν of the stress is active there. b) The nonmonotone case. In Figure 5.4 we plot the deformed configuration as well as the interface forces on Γ 3 for α = 150 and β = 0.5 while in Figure 5.5 that for α = 40 and β = 0.5. Here we are in the nonmonotone case, since β < 1.

We note that in Figure 5.4 all nodes are in status of normal compliance, i.e. 0 ≤ u ν < 0.02. Nevertheless for part of the nodes we have 0 ≤ u ν < 0.01 and for the other part we have 0.01 ≤ u ν < 0.02. We note that for 0 ≤ u ν < 0.01 the normal forces are increasing with respect to the penetration and for 0.01 ≤ u ν < 0.02 they decrease. This behaviour represent the softening property of the deformable layer. It arises since there -σ ν = k ν (u ν ) and k ν is an increasing function on [0, 0.01], and it is decreasing on [0, 0.02].

In Figure 5.5 for part of the nodes we have 0 ≤ u ν ≤ 0.01, for other part we have 0.01 ≤ u ν ≤ 0.02 and, finally, for the remainder part we have u ν = 0.02. We note that for 0 ≤ u ν < 0.01 the normal forces are increasing with respect to the penetration and for 0.01 ≤ u ν < 0.02 they decrease. Moreover, the interface forces increase when u ν = 0.02 since, there, the component σ 3 of the stress is active. Finally, in Figure 5.6 we present our numerical results in the case α = 40 and β = 0. In this case, the non monotonicity reached its peak and, therefore, the average iterations number of the "convexification" procedure to solve the problem was larger that the number of iterations we needed in the previous cases. Note that the value of σ ν on the contact boundary decreases when 0.01 ≤ u ν ≤ 0.02 and converge to zero for a node located in the transition area between the normal compliance zone and the unilateral constraint zone. c) Numerical solution errors. We report relative numerical solution errors in the energy norm u ref -u h E / u ref E in Table 5.1 and Figure 5.7, where the energy norm is defined by

v E := (F(ε(v)), ε(v)) 1/2 Q
which is equivalent to the norm v V . Since the true solution u is not available, we use instead the numerical solution corresponding to a fine discretization of Ω The curve of the relative numerical solution errors is asymptotically linear, which is consistent with the theoretically predicted optimal linear convergence of the numerical solution established in Section 4. Numerical solution errors. We report relative numerical solution errors in the energy norm in Table 5.2 and Figure 5.9. Again, we use the numerical solution with h = 1/256 as the "reference" solution u ref in computing the solution errors. This fine discretization for the "reference" solution u ref corresponds to a problem with 132098 degrees of freedom, 131072 elements and was computed in 5700 CPU time (expressed in seconds) on an IBM computer equipped with Intel Dual core processors (Model 5148, 2.33 GHz). The curve of the relative numerical solution errors is asymptotically linear, which is consistent with the theoretically predicted optimal linear convergence of the numerical solution established in Section 4.

Theorem 2 . 1 .

 21 Under assumption (A 1 )-(A 8 ), Problem (P) has a unique solution u ∈ K.

  c) Frictionless contact with normal compliance. The frictionless version of the contact boundary conditions (4.25) are (4.27)

Fig. 5 . 1 .

 51 Fig. 5.1. Reference configuration of the two-dimensional body

Fig. 5 . 2 .

 52 Fig. 5.2. Problem (P 4 ): Deformed mesh and interface forces for a normal compliance case for α = 40 and β = 2

Fig. 5 . 3 .

 53 Fig. 5.3. Problem (P 4 ): Deformed mesh and interface forces for a normal compliance and unilateral constraint case for α = 10 and β = 2

Fig. 5 . 4 .

 54 Fig. 5.4. Problem (P 4 ): Deformed mesh and interface forces for a normal compliance case for α = 150 and β = 0.5

Fig. 5 . 5 . 5 Fig. 5 . 6 .

 55556 Fig. 5.5. Problem (P 4 ): Deformed mesh and interface forces for a normal compliance and unilateral constraint case for α = 40 and β = 0.5

Fig. 5 . 7 .

 57 Fig. 5.7. Problem (P 4 ): Energy norm errors for frictional contact with normal compliance and unilateral constraint

Fig. 5 . 8 . 5 h 1 2 Fig. 5 . 9 .

 5851259 Fig. 5.8. Problem (P 2 ): Deformed mesh and interface forces for a normal compliance case with friction for α = 150 and β = 0.5 h 1/4 1/8 1/16 1/32 1/64 relative error 0.212 0.104 0.0519 0.0261 0.0128Table 5.2 Problem (P 2 ): Errors in energy norm for a normal compliance case with friction

(5.2)

if r ∈ [0, 0.01], 0.02 -r if r ∈ (0.01, 0.02], r -0.02 if r > 0.02.

Here and below r + represents the positive part of r, i.e. r + = max {r, 0} and α > 0, β ≥ 0 are stiffness coefficients of the foundation.

Obviously, q ν is a Lipschitz continuous increasing function; p ν is Lipschitz continuous and is not monotone. Define the functions ϕ ν : R → R and j ν : R → R by equalities (4.37). It follows that ϕ ν (r) = αβ(r + ) 2 /2 for all r ∈ R which, clearly, satisfies condition (4.42). Note also that j ν is a not a convex function. Assumption (4.8)(a) is obviously satisfied. It is also easy to see that (4.8)(b) and (4.8)(c) hold. By

and is thus also equivalent to the statement that the function

is nondecreasing for some α jν ≥ 0. It is easy to show that for p ν given by (5.2) and α jν = α the function (5.4) is nondecreasing. It follows from here that j ν satisfies the hypothesis (4.8)(d) with α jν = α. Assume now that α is chosen to be sufficiently small. Then it follows that the smallness condition (4.44) holds, too. We conclude from above that the choice (5.1)-(5.2) leads to a a model of contact for which the weak formulation is given by the variational-hemivariational inequality (4.43) and, in addition, the assumptions of Theorem 4.3 are satisfied. Moreover, using (5.3) we deduce that Problem (P 4 ) can be written, in an equivalent form as follows: find u such that

We examine in what follows the feature of the inequality (5.