Weimin Han
email: han@uiowa.edu

Mircea Sofonea
email: sofonea@univ-perp.fr

AND Mika Ël Barboteu
email: barboteu@univ-perp.fr

NUMERICAL ANALYSIS OF ELLIPTIC HEMIVARIATIONAL INEQUALITIES *

Keywords: elliptic hemivariational inequality, Clarke subdifferential, Galerkin approximation, finite element method, convergence, error estimates, contact problems AMS subject classifications. 65N30, 65N15, 74M10, 74M15

This paper is devoted to a study of the numerical solution of elliptic hemivariational inequalities with or without convex constraints by the finite element method. For a general family of elliptic hemivariational inequalities that facilitates error analysis for numerical solutions, the solution existence and uniqueness are proved. The Galerkin approximation of the general elliptic hemivariational inequality is shown to converge, and Céa's inequality is derived for error estimation. For various elliptic hemivariational inequalities arising in contact mechanics, we provide error estimates of their numerical solutions, which are of optimal order for the linear finite element method, under appropriate solution regularity assumptions. Numerical examples are reported on using linear elements to solve sample contact problems, and the simulation results are in good agreement with the theoretically predicted linear convergence.

1. Introduction. The mathematical theory of variational inequalities started in the sixties (cf. [START_REF] Lions | Variational inequalities[END_REF][START_REF] Brézis | Equations et inéquations non linéaires dans les espaces vectoriels en dualité[END_REF]). Since then, the area of variational inequalities has received a lot of attention. Results on mathematical theory and numerical solutions of variational inequalities are summarized in numerous books (e.g., [START_REF] Duvaut | Inequalities in Mechanics and Physics[END_REF][START_REF] Kinderlehrer | An Introduction to Variational Inequalities and Their Applications[END_REF][START_REF] Glowinski | Numerical Analysis of Variational Inequalities[END_REF][START_REF] Friedman | Variational Principles and Free-Boundary Problems[END_REF][START_REF] Baiocchi | Variational and Quasivariational Inequalities: Applications to Free-Boundary Problems[END_REF][START_REF] Chipot | Variational Inequalities and Flow in Porous Media[END_REF][START_REF] Glowinski | Numerical Methods for Nonlinear Variational Problems[END_REF][START_REF] Panagiotopoulos | Inequality Problems in Mechanics and Applications[END_REF][START_REF] Hlaváček | Solution of Variational Inequalities in Mechanics[END_REF]). Variational inequalities with particular emphasis on applications in contact mechanics and engineering are documented in [START_REF] Kikuchi | Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods[END_REF][START_REF] Han | Quasistatic Contact Problems in Viscoelasticity and Viscoplasticity[END_REF][START_REF] Shillor | Models and Analysis of Quasistatic Contact[END_REF][START_REF] Eck | Unilateral Contact Problems: Variational Methods and Existence Theorems[END_REF][START_REF] Sofonea | Analysis and Approximation of Contact Problems with Adhesion or Damage[END_REF][START_REF] Wriggers | Computational Contact Mechanics[END_REF][START_REF] Sofonea | Mathematical Models in Contact Mechanics[END_REF][START_REF] Han | Plasticity: Mathematical Theory and Numerical Analysis[END_REF]. Variational inequalities are mathematical problems with convex structures, and their analysis needs tools and techniques from convex analysis, including arguments of monotonicity, notion, and properties of the subdifferential of a convex function. In contrast, hemivariational inequalities are mathematical problems with nonconvex structures and are particularly useful for analyzing and solving some families of nonsmooth and nonconvex problems. The notion of hemivariational inequalities was first introduced by Panagiotopoulos in the early 1980s [START_REF] Panagiotopoulos | Nonconvex energy functions, hemivariational inequalities and substationary principles[END_REF] and is closely related to the development of the concept of the generalized gradient of a locally Lipschitz function provided by Clarke [START_REF] Clarke | Generalized gradients and applications[END_REF][START_REF] Clarke | Optimization and Nonsmooth Analysis[END_REF]. During the last three decades, hemivariational inequalities were shown to be very useful across a wide variety of subjects, ranging from nonsmooth mechanics, physics, and engineering to economics. For this reason, a large number of problems in applications lead to mathematical models expressed in terms of hemivariational inequalities. The mathematical literature dedicated to this field is growing rapidly. The theory and applications of hemivariational inequalities can be found in several books [START_REF] Panagiotopoulos | Hemivariational Inequalities[END_REF][START_REF] Naniewicz | Mathematical Theory of Hemivariational Inequalities and Applications[END_REF][START_REF] Haslinger | Finite Element Method for Hemivariational Inequalities[END_REF][START_REF] Goeleven | Variational and Hemivariational Inequalities: Theory, Methods and Applications, Volume I: Unilateral Analysis and Unilateral Mechanics[END_REF][START_REF] Carl | Nonsmooth Variational Problems and Their Inequalities: Comparison Principles and Applications[END_REF][START_REF] Migórski | Nonlinear Inclusions and Hemivariational Inequalities. Models and Analysis of Contact Problems[END_REF][START_REF] Han | Advances in Variational and Hemivariational Inequalities: Theory, Numerical Analysis, and Applications[END_REF] and the references cited therein.

Despite the substantial progress made on modeling and mathematical theories of hemivariational inequalities, only a handful of publications can be found in the literature on numerical analysis of hemivariational inequalities. The early monograph [START_REF] Haslinger | Finite Element Method for Hemivariational Inequalities[END_REF] provides discussions on the finite element method for solving hemivariational inequalities; however, no error estimation is given. In [START_REF] Barboteu | An analytical and numerical approach to a bilateral contact problem with nonmonotone friction[END_REF], a bilateral contact problem with nonmonotone friction law is discussed, and an error estimate is derived for the finite element solutions of the corresponding problem expressed in terms of an elliptic hemivariational inequality; yet, the convergence order presented there is suboptimal. More recently, in [START_REF] Han | A class of variational-hemivariational inequalities with applications to frictional contact problems[END_REF], a variational-hemivariational inequality is discussed theoretically and numerically, and an optimal first order error estimate is derived for the linear finite element solutions of the problem. This appears to be the first time in the literature where the optimal first order error estimate is established for the linear finite element method for solving a hemivariational inequality. Then, in [START_REF] Barboteu | Numerical analysis of a hyperbolic hemivariational inequality arising in dynamic contact[END_REF], numerical analysis is performed for solving a hyperbolic hemivariational inequality arising in dynamic frictional contact, and an optimal order error estimate is derived for the linear finite element solutions of the problem. In [START_REF] Kalita | Convergence of Rothe scheme for hemivariational inequalities of parabolic type[END_REF][START_REF] Kalita | Regularity and Rothe method error estimates for parabolic hemivariational inequality[END_REF], a temporal semidiscrete scheme based on the backward Euler difference approximation of the time derivative for solving a parabolic hemivariational inequality is analyzed. Convergence of the semidiscrete solutions is shown in [START_REF] Kalita | Convergence of Rothe scheme for hemivariational inequalities of parabolic type[END_REF], and an error estimate is derived in [START_REF] Kalita | Regularity and Rothe method error estimates for parabolic hemivariational inequality[END_REF] which is not of optimal order with respect to the time step-size. The convergence proof is extended in [START_REF] Kalita | Semidiscrete Variable Time-Step θ-Scheme for Nonmonotone Evolution Inclusion[END_REF] to a family of θ-schemes for the time discretization of the parabolic hemivariational inequality, where θ ∈ (0, 1], with θ = 1 for the backward Euler scheme and θ =1 /2 corresponding to the Crank-Nicolson scheme. In [START_REF] Bartosz | Rothe method for parabolic variational-hemivariational inequalities[END_REF], the numerical solution of parabolic variational-hemivariational inequalities is considered. In [START_REF] Bartosz | Numerical methods for evolution hemivariational inequalities[END_REF], numerical methods are discussed for solving evolutionary hemivariational inequalities, where the time derivatives are approximated by finite differences and the spatial discretization is done with the linear finite element. Convergence of the numerical methods is discussed, and optimal order error estimates are presented.

The main purpose of this paper is to study the numerical solution of general elliptic hemivariational inequality problems by the finite element method. We show the convergence of the numerical solution, and for some particular 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.

The rest of the paper is organized as follows. In section 2 we review some preliminary material needed in the study of hemivariational inequalities. In the literature (e.g., [37, section 4.3], [START_REF] Barboteu | An analytical and numerical approach to a bilateral contact problem with nonmonotone friction[END_REF]), there are solution existence and uniqueness results for some elliptic hemivariational inequalities. Nevertheless, we discuss in section 3 a general family of elliptic hemivariational inequalities in a form that facilitates numerical analysis. The family includes elliptic hemivariational inequalities without constraints, as well as those with convex constraints. The existence and uniqueness result for these elliptic hemivariational inequalities is somewhat more general than that in the existing literature. For example, in both [37, section 4.3] and [START_REF] Barboteu | An analytical and numerical approach to a bilateral contact problem with nonmonotone friction[END_REF], only elliptic hemivariational inequalities without constraints are discussed, the nonconvex functional is given in the form of a boundary integral, and moreover, a stronger growth condition on the nonconvex functional is assumed. In comparison, Theorem 3.1 in section 3 covers elliptic hemivariational inequalities with and without constraints; the nonconvex functional is not limited to a concrete form, and moreover, only a mild growth condition is used on the nonconvex functional. In section 4 we introduce Galerkin methods for solving the hemivariational inequalities, prove convergence, and derive Céa-type inequalities useful for error estimation. We note that in [START_REF] Han | A class of variational-hemivariational inequalities with applications to frictional contact problems[END_REF], Céa-type inequalities and optimal order error estimates for the linear element solutions were derived under the additional assumption that the nonconvex functional is Lipschitz continuous; in this paper, we remove this unsatisfactory assumption. Moreover, the inequality problem studied in [START_REF] Han | A class of variational-hemivariational inequalities with applications to frictional contact problems[END_REF] does not contain a constraint, whereas in this paper, inequality problems both with and without constraints are covered. In section 5 we introduce several contact problems, in which the material's behavior is modeled with a nonlinear elastic constitutive and contact conditions in subdifferential forms. We apply our results in earlier sections for the analysis and numerical approximations of the contact problems. Optimal first order error estimates are shown for the linear element method, under appropriate solution regularity assumptions. Finally, in section 6 we present simulation results which provide numerical evidence of our theoretical error estimates.

Preliminaries.

All linear spaces in this paper are assumed to be real. For a normed space X,w ed e n o t eb y • X its norm, by X * its topological dual, and by

•, • X * ×X the duality pairing of X and X * . Sometimes, when no confusion may arise, we simply write •, • instead of •, • X * ×X .W eu s e0 X to represent the zero element of X,2 X * to denote the collection of all the subsets of X * . Weak convergence is indicated by the symbol ⇀. Given two normed spaces X and Y , L(X, Y ) is the space of all linear continuous operators from X to Y . We deal with both single-valued and multivalued operators defined on a normed space X. We start by recalling several definitions for single-valued operators.

Definition 2.1. Consider an operator A : X → X * , which is generally nonlinear. It is bounded if it maps bounded sets of X to bounded sets of X * . It is monotone if Au -Av, uv ≥0 for all u, v ∈ X. It is maximal monotone if it is monotone and Auw, uv ≥0 for any u ∈ X implies that w = Av.I t i s c o e r c i v e i f there exists a function α : R + → R with lim r→+∞ α(r)=+ ∞ such that Au, u ≥ α( u X ) u X for all u ∈ X. It is pseudomonotone if it is bounded and u n ⇀uin X with lim sup n→∞ Au n ,u nu ≤0 imply Au, uv ≤lim inf n→∞ Au n ,u nv for all v ∈ X.

I tc a nb ep r o v e dt h a ta no p e r a t o rA : X → X * is pseudomonotone if and only if it is bounded and u n ⇀uin X together with lim sup n→∞ Au n ,u nu ≤0i m p l y Au n ⇀Auin X * and lim n→∞ Au n ,u n = Au, u .

For a multivalued operator T : X → 2 X * , its graph G(T )i s

G(T ):={(x, x * ) ∈ X × X * | x * ∈ Tx}. Definition 2.2. An operator T : X → 2 X * is monotone if u * -v * ,u-v ≥0 for all (u, u * ), (v, v * ) ∈G(T ).
It is maximal monotone if it is monotone and maximal in the sense of inclusion of graphs in the family of monotone operators from X to 2 X * . It is coercive if there exists a function α : R + → R with lim r→+∞ α(r)=+ ∞ such that u * ,u ≥α( u X ) u X for all (u, u * ) ∈G(T ).

Next, we recall the notions of pseudomonotonicity and generalized pseudomonotonicity for a multivalued operator. Definition 2.3. Let X be a reflexive Banach space. A multivalued operator T : X → 2 X * is pseudomonotone if (a) for every u ∈ X,t h es e tTu ⊂ X * is nonempty, closed, and convex;

(b) T is upper semicontinuous from each finite dimensional subspace of X to X * endowed with the weak topology;

(c) for any sequences {u n }⊂X and {u * n }⊂X * such that u n ⇀uin X, u * n ∈ Tu n for all n ≥ 1,a n dlim sup u * n ,u nu ≤0, we have that for every v ∈ X, there exists u * (v) ∈ Tu such that

u * (v),u-v ≤lim inf u * n ,u n -v .
Definition 2.4. Let X be a reflexive Banach space. A multivalued operator T : X → 2 The following surjectivity result (cf. [START_REF] Naniewicz | Mathematical Theory of Hemivariational Inequalities and Applications[END_REF]Theorem 2.11]) will be applied in studying the elliptic hemivariational inequalities.

Theorem 2.6. Let X be a reflexive Banach space, T 1 : X → 2 X * pseudomonotone and coercive, and

T 2 : X → 2 X * maximal monotone with T 2 (0 X ) = ∅.I fe i t h e rT 1 or T 2 is bounded, then T 1 + T 2 is surjective.
Finally, we recall the definitions of the convex and the Clarke subdifferentials. Definition 2.7. Let ϕ : X → R ∪{+∞} be a proper, convex, and lower semicontinuous function. The mapping ∂ c ϕ : X → 2 X * defined by

∂ c ϕ(x):={x * ∈ X * | x * ,v-x ≤ϕ(v) -ϕ(x) ∀ v ∈ X} is cal led the (convex) subdifferential of ϕ.A ne l e m e n tx * ∈ ∂ c ϕ(x) (if any) is called a subgradient of ϕ at x. Definition 2.8. Let ψ : X → R be a locally Lipschitz functional. 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 gradient (subdifferential) of ψ at x is defined by

∂ψ(x):= ζ ∈ X * | ψ 0 (x; v) ≥ ζ,v ∀v ∈ X .
Details on the properties of the subdifferential mappings, in both the convex and the 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 | Hemivariational Inequalities[END_REF][START_REF] Sofonea | Mathematical Models in Contact Mechanics[END_REF]. In particular, knowing the generalized subdifferential, we can compute the generalized directional derivative through the formula [START_REF] Clarke | Optimization and Nonsmooth Analysis[END_REF] (1)

ψ 0 (x; v)=max{ ζ,v |ζ ∈ ∂ψ(x)} .
3. Analysis of a general elliptic hemivariational inequality. Let X be a reflexive Banach space, K ⊂ X,a n dX j a Banach space. Given an operator A : X → X * , a locally Lipschitz functional j : X j → R, a linear operator γ j : X → X j ,a n da linear functional f : X → R, we consider the following problem.

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

(2) Au, v -u + j 0 (γ j u; γ j v -γ j u) ≥ f, v -u ∀ v ∈ K.
For the study of Problem (P), we adopt 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 j is a Banach space, γ j ∈L(X, X j ): for a constant c j > 0,

(3)

γ j v X j ≤ c j v X ∀ v ∈ X.
(A 3 ) A : X → X * is pseudomonotone and strongly monotone: for a constant

m A > 0, (4) Av 1 -Av 2 ,v 1 -v 2 ≥m A v 1 -v 2 2 X ∀ v 1 ,v 2 ∈ X.
(A 4 ) j : X j → R is locally Lipschitz, and there are constants c 0 ,c 1 ,α j ≥ 0s u c h that

∂j(z) X * j ≤ c 0 + c 1 z X j ∀ z ∈ X j , (5)
j 0 (z 1 ; z 2 -z 1 )+j 0 (z 2 ; z 1 -z 2 ) ≤ α j z 1 -z 2 2 X j ∀ z 1 ,z 2 ∈ X j . (6) (A 5 ) (7) α j c 2 j <m A . (A 6 ) (8) f ∈ X * .
Note that Problem (P) contains as particular cases various problems considered in the literature. The space X j is introduced to facilitate derivation of optimal order error estimates for numerical solutions of Problem (P) in later sections. For applications in contact mechanics with a spatial dimension d, the functional j(•)isanintegral over the contact boundary Γ 3 ,a n dX j can be chosen to be 7) is a smallness assumption, which poses a limit on the size of j relative to the strong monotonicity of A. For a locally Lipschitz function j : X j → R,t h e hypothesis ( 6) is equivalent to the condition

L 2 (Γ 3 )o rL 2 (Γ 3 ) d .T h e relation (
(9) ∂j(z 1 ) -∂j(z 2 ),z 1 -z 2 X * j ×X j ≥-α j z 1 -z 2 2 X j ∀ z 1 ,z 2 ∈ X j .
This is known as the relaxed monotonicity condition and it was extensively used in the literature (e.g., [START_REF] Migórski | Nonlinear Inclusions and Hemivariational Inequalities. Models and Analysis of Contact Problems[END_REF]). In case j : X j → R is convex, then ( 6) and ( 9) hold with α j = 0, due to the monotonicity of the (convex) subdifferential. We call (2) an elliptic hemivariational inequality in view of the assumption (4). When K = X, (2) is an elliptic hemivariational inequality without constraint. When K is a proper convex subset of X, (2) is an elliptic hemivariational inequality with a convex constraint.

The assumption 0 X ∈ K is valid for all the contact problems we consider and is introduced to simplify the exposition. It is possible to replace this assumption by the non-emptiness of K; then instead of Theorem 2.6 that follows from [START_REF] Naniewicz | Mathematical Theory of Hemivariational Inequalities and Applications[END_REF]Theorem 2.11],

we apply [START_REF] Naniewicz | Mathematical Theory of Hemivariational Inequalities and Applications[END_REF]Theorem 2.12] in the existence proof for the elliptic hemivariational inequalities. Since 0 X ∈ K, we derive the following inequalities from ( 4), [START_REF] Brézis | Equations et inéquations non linéaires dans les espaces vectoriels en dualité[END_REF], and (5):

Av, v ≥m A v 2 X -c v X ∀ v ∈ X, (10) ∂j(z),z X * j ×X j ≥-α j z 2 X j -c 0 z X j ∀ z ∈ X j (11)
for some constant c ≥ 0.

We have the following existence and uniqueness result.

Theorem 3.1. Under the assumptions (A 1 )-(A 6 ), Problem (P) has a unique solution u ∈ K.

Proof. We first prove the existence. By making use of the indicator function

I K (v)= 0i f v ∈ K, +∞ if v ∈ X\K,
we can rewrite Problem (P) in the equivalent form: find u ∈ X such that

(12) Au, v -u + I K (v) -I K (u)+j 0 (γ j u; γ j v -γ j u) ≥ f, v -u ∀ v ∈ X.
Now consider the following problem: find u ∈ X such that

(13) Au + γ * j ∂j(γ j u)+∂ c I K (u) ∋ f,
where ∂j ⊂ X * j denotes the generalized gradient of j, ∂ c I K is the convex subdifferential of I K ,a n dγ * j ∈L(X * j ,X * ) is the adjoint of γ j . Corresponding to [START_REF] Clarke | Generalized gradients and applications[END_REF], introduce two multivalued operators T 1 , T 2 : X → 2 X * :

T 1 v = Av + γ * j ∂j(γ j v),T 2 v = ∂ c I K (v).
We know that I K is proper, convex, and lower semicontinuous with an effective domain K. It is well known (cf., e.g., [START_REF] Denkowski | An Introduction to Nonlinear Analysis: Applications[END_REF]Theorem 1.3.19]) that the operator

T 2 = ∂ c I K : X → 2 X * is maximal monotone with D(∂ c I K )=K.
We claim that the operator T 1 is bounded, coercive, and pseudomonotone. The boundedness of T 1 follows from that of A, γ j , γ * j , and the growth condition (5) on ∂j. For the coercivity, we use the inequalities [START_REF] Carl | Nonsmooth Variational Problems and Their Inequalities: Comparison Principles and Applications[END_REF] and [START_REF] Chipot | Variational Inequalities and Flow in Porous Media[END_REF], c, due to the smallness assumption [START_REF] Bartosz | Rothe method for parabolic variational-hemivariational inequalities[END_REF].

T 1 v, v = Av, v + ∂j(γ j v),γ j v ≥ m A v 2 X -c v X -α j γ j v 2 X j -c 0 γ j v X j ≥ m A -α j c 2 j v 2 X -c v X for all v ∈ X.T h u s ,T 1 is coercive with α(v): = m A -α j c 2 j v X -
We now prove that the operator T 1 is pseudomonotone. Observe that for all v ∈ X,t h es e tAv + γ * j ∂j(γ j v) is nonempty, closed, and convex in X * . According to Proposition 2.5, we only need to show that T 1 is generalized pseudomonotone. From (4), [START_REF] Bartosz | Numerical methods for evolution hemivariational inequalities[END_REF], and (7), we know the operator T 1 is strongly monotone, ( 14)

T 1 v 1 -T 1 v 2 ,v 1 -v 2 ≥ m A -α j c 2 j v 1 -v 2 2 X ∀ v 1 ,v 2 ∈ X. Let u n ∈ X, u n ⇀uin X, u * n ∈ T 1 u n , u * n ⇀u * in X * , and lim sup n→∞ u * n ,u n -u ≤ 0. Use (14), m A -α j c 2 j u n -u 2 X ≤ u * n ,u n -u -T 1 u, u n -u ,
and we deduce that u n → u in X.S i n c eu

in X * for some u * a ,u * b ∈ X * . Letting n →∞in u * n = u * n,1 + u * n,2 ,w eh a v eu * = u * a + u * b .
Exploiting the equivalent condition for the pseudomonotonicity of A,w e have Au n ⇀A uin X * , implying u * a = Au.S i n c eX ∋ v → γ * j ∂j(γ j v) ∈ 2 X * has a closed graph with respect to the strong topology in X and weak topology in X * , we infer that u * b ∈ γ * j ∂j(γ j u). Hence,

u * = u * a + u * b ∈ Au + γ * j ∂j(γ j u)=T 1 u.F r o m u * n ⇀u * in X * and u n → u in X,w eha v e u * n ,u n → u * ,u .
Hence, T 1 is generalized pseudomonotone and is then also pseudomonotone.

We then apply Theorem 2.6 and deduce that there exists a solution u ∈ X to the inclusion [START_REF] Clarke | Generalized gradients and applications[END_REF]. This solution also solves the problem [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]. Indeed, let u ∈ X be such that [START_REF] Denkowski | An Introduction to Nonlinear Analysis: Theory[END_REF] Au

+ y * + z * = f with y * ∈ ∂ c I K (u)a n dz * ∈ γ * j ∂j(γ j u). For all v ∈ X,w eh a v e y * ,v-u ≤I K (v) -I K (u), z * ,v ≤j 0 (γ j u; γ j v).
Combining [START_REF] Denkowski | An Introduction to Nonlinear Analysis: Theory[END_REF] with these inequalities, we obtain

Au, v -u + I K (v) -I K (u)+j 0 (γ j u; γ j v -γ j u) ≥ f, v -u ∀ v ∈ X.
So u ∈ X satisfies [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] and we conclude that Problem (P) has at least one solution u ∈ K.

The solution uniqueness is proved by a standard approach and is hence omitted. Note that the smallness assumption [START_REF] Bartosz | Rothe method for parabolic variational-hemivariational inequalities[END_REF] is needed in the uniqueness proof.

4. Galerkin approximation. In this section, we consider numerical schemes for solving Problem (P). We keep assumptions (A 1 )-( A 6 ) so that Problem (P) has a unique solution u ∈ K. In the rest of the paper, we will use c to represent a generic positive constant that is independent of the meshsize h and the solution u and whose value may vary from one place to another.

Let X h ⊂ X be a finite dimensional subspace with h>0d e n o t i n gas p a t i a l discretization parameter. We use K h := X h ∩ K to approximate the convex set K. The Galerkin approximation of Problem (P) is the following.

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

(16) Au h ,v h -u h + j 0 (γ j u h ; γ j v h -γ j u h ) ≥ f, v h -u h ∀ v h ∈ K h .
The arguments of the proof of Theorem 3.1 can be applied in the setting of the finite dimensional set K h , and we know that under the assumptions (A 1 )-( A 6 ), Problem (P h ) has a unique solution u h ∈ K h . The focus of this section is convergence analysis and error estimation for the numerical solution of Problem (P h ). We assume

(17) ∀ v ∈ K, ∃ v h ∈ K h such that v h → v in X as h → 0.
We also assume A : X → X * is Lipschitz continuous, i.e., for some constant

L A > 0, (18) Au -Av X * ≤ L A u -v X ∀ u, v ∈ X.
Note that conditions ( 18) and ( 4) imply the pseudomonotonicity of the operator A [46, Proposition 27.6].

Proposition 4.1. The solution u h ∈ K h of Problem (P h ) is uniformly bounded independent of h.

Proof. We let v h = 0 in [START_REF] Denkowski | An Introduction to Nonlinear Analysis: Applications[END_REF] to get [START_REF] Friedman | Variational Principles and Free-Boundary Problems[END_REF] Au h ,u h ≤j 0 (γ j u h ; -γ j u h )+ f, u h .

From ( 6) and ( 5),

j 0 (γ j u h ; -γ j u h ) ≤ α j γ j u h 2 X j -j 0 (0; γ j u h ) (20) ≤ α j γ j u h 2 X j + c 0 + c 1 γ j u h X j ≤ α j c 2 j u h 2 X + c 1+ u h X .
Apply ( 10) and ( 20) in [START_REF] Friedman | Variational Principles and Free-Boundary Problems[END_REF],

m A -α j c 2 j u h 2 X ≤ c 1+ u h X .
We then use the elementary implication for real numbers

a, b, x ≥ 0a n dx 2 ≤ ax+ b ⇒ x 2 ≤ a 2 +2b
to conclude that u h ∈ K h is uniformly bounded independent of h.

4.1. Convergence. We begin with an application of (4) 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 (21) + Au, u -u h + Au h ,u h -v h . From (2) with v = u h , (22) Au, u -u h ≤j 0 (γ j u; γ j u h -γ j u) -f, u h -u .
From ( 16),

(23) Au h ,u h -v h ≤j 0 (γ j u h ; γ j v h -γ j u h ) -f, v h -u h .
Using ( 22) and ( 23) in ( 21), we have

m A u -u h 2 X ≤ Au -Au h ,u-v h + Au, v h -u -f, v h -u (24)
+ j 0 (γ j u; γ j u hγ j u)+j 0 (γ j u h ; γ j v hγ j u h ).

For error estimation, it will be more convenient to rewrite (24) as ( 25)

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

R(v h ):= Au, v h -u + j 0 (γ j u; γ j v h -γ j u) -f, v h -u , (26)
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). (27)
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 ( 25) that

(28) (m A -ε) u -u h 2 X ≤ c u -v h 2 X + R(v h )+I j (v h ).
We will apply repeatedly the subadditivity of the generalized directional derivative:

j 0 (z; w 1 + w 2 ) ≤ j 0 (z; w 1 )+j 0 (z; w 2 ) ∀ z, w 1 ,w 2 ∈ X j . Using j 0 (γ j u; γ j u h -γ j u) ≤ j 0 (γ j u; γ j u h -γ j v h )+j 0 (γ j u; γ j v h -γ j u),
we first bound the term I j (v h )a s

I j (v h ) ≤ j 0 (γ j u; γ j u h -γ j v h )+j 0 (γ j u h ; γ j v h -γ j u h ).
Then, j 0 (γ j u; γ j u hγ j v h ) ≤ j 0 (γ j u; γ j uγ 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 h ; γ j v hγ j u)+j 0 (γ j u h ; γ j uγ j u h ), and 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 ) (29) + j 0 (γ j u; γ j u -γ j v h )+j 0 (γ j u h ; γ j v h -γ j u) .
By [START_REF] Bartosz | Numerical methods for evolution hemivariational inequalities[END_REF],

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 X j .
Thus, (30)

I j (v h ) ≤ α j c 2 j u -u h 2 X + j 0 (γ j u; γ j u -γ j v h )+j 0 (γ j u h ; γ j v h -γ j u).
From (1) and ( 5),

j 0 (γ j u; γ j u -γ j v h ) ≤ c 0 + c 1 γ j u X j γ j u -γ j v h X j , j 0 (γ j u h ; γ j v h -γ j u) ≤ c 0 + c 1 γ j u h X j γ j u -γ j v h X j .
Note that γ j u h X j is bounded by a constant independent of h (Proposition 4.1). Thus, (31)

I j (v h ) ≤ α j c 2 j u -u h 2 X + c γ j u -γ j v h X j .
Combine ( 28) and [START_REF] Kalita | Semidiscrete Variable Time-Step θ-Scheme for Nonmonotone Evolution Inclusion[END_REF],

m A -α j c 2 j -ε u -u h 2 X ≤ c u -v h 2 X + c γ j u -γ j v h X j + R(v h ).
Since α j c 2 j <m A ,w ec a nc h o o s eε =(m Aα j c 2 j )/2 > 0 and get

(32) u -u h 2 X ≤ c u -v h 2 X + γ j u -γ j v h X j + R(v h ) ∀ v h ∈ K h .
This is a starting point for convergence statement (Theorem 4.2 below) and error estimation (cf. subsection 5.3).

Theorem 4.2. Under the assumptions (A 1 )-(A 6 ), ( 17),a n d( 18), we have the convergence of the numerical solution defined by Problem (P h ) to the solution of Problem (P).

Proof. B y( 1 7 ) ,w eh a v eas e q u e n c e{v h }, v h ∈ K h ,t h a tc o n v e r g e st ou in X as h → 0. The residual term R(v h ) can be bounded as follows:

R(v h ) ≤ Au X * u -v h X + c 0 + c 1 γ j u X j γ j u -γ j v h X j + f X * u -v h X .
Note that

γ j u -γ j v h X j ≤ c j u -v h X . Then, lim h→0 R(v h )=0,
and we conclude the convergence uu h X → 0a sh → 0 from (32).

4.2.

Error estimation for numerical solutions of the problem without constraint. In the special case K = X,w eh a v eK h = X h , and the original problem (2) and its approximation ( 16) become [START_REF] Kikuchi | Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods[END_REF] Au, v + j 0 (γ j u;

γ j v) ≥ f, v ∀ v ∈ X and (34) Au h ,v h + j 0 (γ j u h ; γ j v h ) ≥ f, v h ∀ v h ∈ X h .
We replace v by uv in [START_REF] Kikuchi | Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods[END_REF],

Au, u -v + j 0 (γ j u; γ j u -γ j v) ≥ f, u -v ∀ v ∈ X. Thus, Au, v h -u ≤j 0 (γ j u; γ j u -γ j v h ) -f, u -v h ∀ v h ∈ X h .
Using this inequality in [START_REF] Han | Advances in Variational and Hemivariational Inequalities: Theory, Numerical Analysis, and Applications[END_REF], we have

(35) m A u -u h 2 X ≤ Au -Au h ,u-v h + Ĩj (v h ),
where (36) Ĩ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 ).

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 write Ĩ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) .

The right side of the above inequality is the same as that of [START_REF] Kalita | Regularity and Rothe method error estimates for parabolic hemivariational inequality[END_REF]. Then, by ( 31), ( 37)

Ĩj (v h ) ≤ α j c 2 j u -u h 2 X + c γ j u -γ j v h X j .
Combining ( 35) and ( 37), we get

m A u -u h 2 X ≤ Au -Au h X u -v h X + α j c 2 j u -u h 2 X + c γ j u -γ j v h X j .
Using the smallness assumption [START_REF] Bartosz | Rothe method for parabolic variational-hemivariational inequalities[END_REF], similar to (32), we deduce from the above inequality that

(38) u -u h 2 X ≤ c u -v h 2 X + γ j u -γ j v h X j ∀ v h ∈ X h .
This is a Céa's inequality and is a basis for deriving error estimates (cf. subsections 5.1 and 5.2).

5.

Error analysis for contact problems. We illustrate applications of the framework developed in section 4 on convergence and error estimation for the finite element solutions of three sample 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 ,a n dΓ 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 ,i sfi x e do nΓ 1 , and is in potential contact on Γ 3 with a foundation. Different contact conditions will lead to different contact problems, as discussed below.

We use S d for the space of second order symmetric tensors on R d and " • "t h e canonical inner product on the spaces R d and S d . W ed e n o t eb yu :Ω→ R d and σ :Ω→ S d the displacement field and the stress field, respectively. The linearized strain tensor associated with u is denoted by ε(u). Let ν be the unit outward normal vector, defined a.e. on Γ. For a vector field v,w eusev ν := v • ν and v τ := vv ν ν for the normal and tangential components of v on Γ. Similarly, for the stress field σ,i t s 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 [START_REF] Panagiotopoulos | Nonconvex energy functions, hemivariational inequalities and substationary principles[END_REF] σ = Fε(u)i n Ω , the equilibrium equation Here and below we do not always indicate explicitly the dependence of a quantity on the spatial variable x. In [START_REF] Panagiotopoulos | Nonconvex energy functions, hemivariational inequalities and substationary principles[END_REF], F :Ω× S d → S d is the elasticity operator and is assumed to have the following properties:

(43)

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (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 S d )=0 S d for a.e. x ∈ Ω.
The relations ( 39)-( 42) will be supplemented by a set of boundary conditions on Γ 3 .

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;
the associated norm is denoted by • Q . The displacement fields will be sought in the space

(44) V = v ∈ 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 .F o rv ∈ H 1 (Ω; R d )w eu s et h es a m es y m b o lv for its trace on Γ. By the Sobolev trace theorem, we have a constant c>0 such that

v L 2 (Γ 3 ;R d ) ≤ c v V ∀ v ∈ V.
We assume the densities of body forces and surface tractions satisfy

(45) f 0 ∈ L 2 (Ω; R d ), f 2 ∈ L 2 (Γ 2 ; R d ) and define f ∈ V * by (46) f , v V * ×V =(f 0 , v) L 2 (Ω;R d ) +(f 2 , v) L 2 (Γ 2 ;R d ) ∀ v ∈ V.
As examples of Problem (P), we consider three choices of the boundary conditions on the contact boundary Γ 3 , leading to different contact problems. For each problem, convergence of the numerical solutions follows from Theorem 4.2. Thus, we will focus on the derivation of error bounds.

A bilateral contact problem with friction.

The first set of contact boundary conditions we consider is

(47) u ν =0, -σ τ ∈ ∂j τ (u τ )o n Γ 3 .
The feature of bilateral contact is reflected by the condition u ν = 0. Note that a problem similar to the one described by ( 39)-( 42) and (47) has already been studied in [26, p. 177], where F is a linear operator, and in [44, p. 144], where F satisfies conditions [START_REF] Sofonea | Analysis and Approximation of Contact Problems with Adhesion or Damage[END_REF]. In these references, the function j τ is assumed to be convex. The novelty of our results below in this subsection is that we extend the numerical analysis of those models to the nonconvex case. To this end, on the potential function

j τ :Γ 3 × R d → R, we assume (48) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (a) j τ (•, z) is measurable on Γ 3 for all z ∈ R d and j τ (•, z 0 (•)) ∈ L 1 (Γ 3 )f o rs o m ez 0 ∈ L 2 (Γ 3 ) d ; (b) j τ (x, •) is locally Lipschitz on R d for a.e. x ∈ Γ 3 ; (c) |∂j τ (x, z)|≤c 0 + c 1 z for a.e. x ∈ Γ 3 , for all z ∈ R d with c 0 , c 1 ≥ 0; (d) j 0 τ (x, z 1 ; z 2 -z 1 )+j 0 τ (x, z 2 ; z 1 -z 2 ) ≤ α j τ z 1 -z 2 2 for a.e. x ∈ Γ 3 , all z 1 , z 2 ∈ R d with α j τ ≥ 0.
We comment that in the case where j τ : R d → R is independent of x, the condition (48)(a) can be dropped.

To apply Theorem 3.1, we let

X = K = V 1 := {v ∈ V | v ν =0a.e. on Γ 3 } , (49)
X j = L 2 (Γ 3 ) d ,γ j v = v τ for v ∈ V 1 , (50)
and define 1,V , λ 1,V > 0b e i n gt h e smallest eigenvalue of the eigenvalue problem

Au, v =(F(ε(u h )), ε(v h )) Q , u, v ∈ V, (51)
j(z)= Γ 3 j τ (•, z(•)) dΓ, z ∈ X j . (52)
u ∈ V 1 , Ω ε(u)•ε(v) dx = λ Γ 3 u τ •v τ dΓ ∀ v ∈ V 1 .
Following 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 ∈ V 1 such that

(54) (F(ε(u)), ε(v)) Q + Γ 3 j 0 τ (u τ ; v τ ) dΓ ≥ f , v V * ×V ∀ v ∈ V 1 .
By Theorem 3.1, we know that assuming additionally (55)

α j τ <λ 1/2 1,V m F , there is a unique element u ∈ V 1 satisfying (F(ε(u)), ε(v)) Q + j 0 (u τ ; v τ ) ≥ f , v V * ×V ∀ v ∈ V 1 .
Then the contact condition -σ τ ∈ ∂j τ (u τ ) from ( 47) is equivalent to

(63) σ τ ≤µ(0) if u τ = 0, -σ τ = µ( u τ ) u τ u τ if u τ = 0 on Γ 3 .
Here, µ(t) can be interpreted as a friction bound function. This function is assumed to be measurable from [0, ∞)toR, µ(0+) > 0, and with two positive constants c 1 ,c 2 ,

0 ≤ µ(t) ≤ c 1 (1 + t) ∀ t ≥ 0, µ(t 2 ) -µ(t 1 ) ≥-c 2 (t 2 -t 1 ) ∀ t 2 >t 1 ≥ 0.
Then (48) holds with α j τ = c 2 .

5.2.

A frictionless normal compliance contact problem. Here, the contact boundary conditions are

(64) -σ ν ∈ ∂j ν (u ν ), σ τ = 0 on Γ 3 .
The first relation in ( 64) is a normal compliance contact condition, whereas the second relation indicates that the contact is frictionless. The problem described by ( 39)-( 42) and ( 64) represents the frictionless version of a nonlinear elastic contact model studied in [37, p. 202]. In that reference, the unique solvability of the model is provided but there is no numerical analysis of the problem. Here, we fill this gap by providing both the numerical analysis and numerical simulations in the study of such a contact model. To this end, we assume the following properties on the potential function

j ν :Γ 3 × R → R: (65) ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (a) j ν (•,z) is measurable on Γ 3 for all z ∈ R and there exists z 0 ∈ L 2 (Γ 3 ) such that j ν (•,z 0 (•)) ∈ L 1 (Γ 3 ); (b) j ν (x, •) is locally Lipschitz on R for a.e. x ∈ Γ 3 ; (c) |∂j ν (x,z)|≤c 0 + c 1 |z| for a.e. x ∈ Γ 3 , for all z ∈ R with c 0 , c 1 ≥ 0; (d) j 0 ν (x,z 1 ; z 2 -z 1 )+j 0 ν (x,z 2 ; z 1 -z 2 ) ≤ α j ν |z 1 -z 2 | 2 for a.e. x ∈ Γ 3 , all z 1 ,z 2 ∈ R with α j ν ≥ 0.
Again, in case j ν : R → R is independent of x, the condition (65) (a) can be dropped. The weak formulation of the contact problem is the following.

Problem (P 2 ). Find a displacement field u ∈ V such that

(66) (F(ε(u)), ε(v)) Q + Γ 3 j 0 ν (u ν ,v ν ) dΓ ≥ f , v V * ×V ∀ v ∈ V.
We use the finite element space V h of (56) and introduce the following approximation of Problem (P 2 ).

Problem (P h 2 ). Find a displacement field u h ∈ V h such that

(67) (F(ε(u h )), ε(v h )) Q + Γ 3 j 0 ν (u h ν ,v h ν ) dΓ ≥ f , v h V * ×V ∀ v h ∈ V h .
Discussion of Problem (P 2 ) and error analysis for Problem (P h 2 ) are similar to that in subsection 5.1, with the following modifications:

X = K = V, X j = L 2 (Γ 3 ),γ j v = v ν for v ∈ V, (68) j(z)= Γ 3 j ν (•,z(•)) dΓ,z ∈ X j , (69)
Finally, from [START_REF] Kalita | A class of subdifferential inclusions for elastic unilateral contact problems[END_REF], we have Céa's inequality

(84) u -u h 2 V ≤ c u -v h 2 V + u ν -v h ν L 2 (Γ 3 ) ∀ v h ∈ U h .
Recalling the solution regularity (72), an assumption made earlier in order to derive the pointwise equations ( 80) and (81), we conclude the optimal order error bound (85) uu h V ≤ ch.

Remark 5.1. The same technique can be applied to the problems without constraints that are studied in subsections 5.1a n d5 .2. However, in deriving Céa's inequalities (59) and (71), there is no need to assume the solution regularity. In other words, under solution regularity assumptions weaker than (61) or (72), we are still able to derive corresponding error estimates for the numerical solutions of the problems in subsections 5.1a n d5 .2. On the other hand, to derive Céa's inequality for the contact problem in this subsection, we need to first assume the solution regularity (72). Thus, it is meritorious to consider the particular case K = X in section 4.

Numerical examples.

In this section we report simulation results for three numerical examples, corresponding to the contact problems (P 1 )-(P 3 ). The solution of the discrete problems is based on numerical methods presented 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]. The main ingredient of these methods is a "convexification" iterative procedure which approximates the solution of a nonconvex problem by solutions of a sequence of convex problems. The nonsmooth convex problems are solved by classical numerical methods that can be found, for instance, in [START_REF] Laursen | Computational Contact and Impact Mechanics[END_REF][START_REF] Wriggers | Computational Contact Mechanics[END_REF].

For the three numerical examples, we use the same physical setting as depicted in Figure 1, with different contact boundary conditions. The domain Ω represents the cross section of a three-dimensional linear elastic body such that the plane stress hypothesis is valid. For simulations, we take Ω to be the unit square: Ω = (0, 1) × (0, 1) ⊂ R 2 ,a n dl e t

Γ 3 contact interface Γ 2 f 2 f 2 Ω deformable body
Γ 1 =[0, 1] ×{1}, Γ 2 =({0}×(0, 1)) ∪ ({1}×(0, 1)), Γ 3 =[0, 1] ×{0}.
T h eb o d yi sc l a m p e do nΓ 1 and is subject to the actions of a vertical body force of constant density and of horizontally compressive forces on the part ({0}×[0.5, 1)) ∪ ({1}×[0.5, 1)) of the boundary Γ 2 . The part ({0}×(0, 0.5))∪({1}×(0, 0.5)) is traction free. The body is in contact with an obstacle on Γ 3 . For numerical simulations, linear finite elements on uniform triangulations of the domain Ω are used. The boundary of the spatial domain is divided into 1/h equal parts, and h is used as the discretization parameter.

The mechanical response of the material is described by a linear elastic constitutive law. The components of the elasticity tensor F are

(Fτ ) ij = Eκ 1 -κ 2 (τ 11 + τ 22 ) δ ij + E 1+κ τ ij , 1 ≤ i, j ≤ 2, ∀ τ ∈ S 2 ,
E, κ,a n dδ αβ being the Young modulus, the Poisson ratio of the material and the Kronecker symbol, respectively. For the numerical simulations, we use

E = 2000N/m 2 ,κ =0.4, f 0 =(0, -0.5 × 10 -3 )N/m 2 , f 2 = (8 × 10 -3 , 0) N/m on {0}×[0.5, 1), (-8 × 10 -3 , 0) N/m
on {1}×[0.5, 1).

Example 6.1. This is an example of Problem (P 1 ). The contact is bilateral and is frictional. We use the friction law (63) in which Figure 2 shows a typical deformed mesh and interface forces on Γ 3 .W eo b s e r v et h a t a large proportion of contact nodes situated at the extremities of Γ 3 are in the status of slip since the friction bound is reached there. In addition, the nodes situated in the center of Γ 3 are in the status of stick. The numerical solution shown in Figure 2 corresponds to a meshsize h =1 /64: the spatial domain is discretized into 16449 elements with 64 contact elements; the total number of degrees of freedom is equal to 16772. The average iterations number of the "convexification" procedure for the solution of discrete problem was 22. The simulation was completed in 586 CPU time (expressed in seconds) on an IBM computer equipped with Intel Dual core processors (Model 5148, 2.33 GHz).

In Table 1 and Figure 3, we report relative errors of the numerical solutions in the energy norm, u refu h E / u ref E , where the energy norm is defined by the formula Note that the energy norm v E is equivalent to the norm v V , and the error bound (60) predicts an optimal first order convergence of the numerical solutions measured in the energy norm, under the regularity assumptions (61). Since the true solution u is not available, we use the numerical solution corresponding to a fine discretization of Ω as the "reference" solution u ref in computing the solution errors.

v E := 1 √ 2 (F(ε(v)), ε(v))
Here, the numerical solution with h =1 /256 is taken to be the "reference" solution u ref . This fine discretization corresponds to a problem with 132612 degrees of freedom and 131329 elements. We clearly observe the theoretically predicted optimal linear convergence of the numerical solutions.

the numerical simulation, g = -0.15. For the conditions (89) we use the multivalued normal compliance response (87) from Example 2.

As in the two previous examples, we plot in Figure 5 the deformed mesh as well as the interface forces on Γ 3 for h =1/64. On the extremities of the boundary Γ 3 ,w e can see that the contact nodes are in multivalued normal compliance status. At the center of Γ 3 , the nodes are in unilateral contact status since, there, the penetration reached the maximum value g. In this case, the average iteration number of the "convexification" procedure for the solution of the discrete problems was 9 and the simulation runs in 267 CPU time.

Table 2 and Figure 6 provide relative errors of numerical solutions in the energy norm, similar to that in Example 1. Again, we clearly observe the theoretically predicted optimal linear convergence of the numerical solutions.

( 40 )

 40 Div σ + f 0 = 0 in Ω, the displacement boundary condition (41) u = 0 on Γ 1 , and the traction boundary condition (42) σν = f 2 on Γ 2 .

  z(•); w(•)) dΓ, z, w ∈ X j . Then (A 1 )-( A 4 )a n d( A 6 ) are satisfied with m A = m F from (43)(b), α j = α j τ from (48)(d). Note that the pseudomonotonicity of A follows from (43)(a) and (43)(b) [46, Proposition 27.6]. The inequality (3) holds for any c j ≥ λ -1/2

Fig. 1 .

 1 Fig. 1. Reference configuration of the two-dimensional body.

  (86) µ(t)=(ab) e -βt + b with a =0 .4, b =0 .2, and β = 2000. The friction bound decreases with the slip from the value a to the limit value b. The corresponding friction law is nonmonotone.

Fig. 2 .

 2 Fig. 2. Deformed meshes and interface forces on Γ 3 for Example 1.

Fig. 3 .

 3 Fig. 3. Relative errors in energy norm for Example 1.

Fig. 5 .

 5 Fig. 5. Deformed meshes and interface forces on Γ 3 for Example 3.

Fig. 6 .

 6 Fig.6. Relative errors for Example 3.

  X * is generalized pseudomonotone if for any sequences {u n }⊂X and {u * n }⊂X * such that u n ⇀uin X, u * n ∈ Tu n for n ≥ 1, u * n ⇀u * in X * ,a n d lim sup u * n ,u nu ≤0, we have u * ∈ Tu and Proposition 2.5. Let X be a reflexive Banach space and T : X → 2 X * . (a) If T is pseudomonotone, then it is generalized pseudomonotone. (b) If T is a bounded, generalized pseudomonotone operator such that for all u ∈ X, Tu is a nonempty, closed, and convex subset of X * ,t h e nT is pseudomonotone.

	lim n→∞	u * n ,u n = u * ,u .
	The following result displays relations between the two notions (cf. [16, Proposi-
	tions 1.3.65 and 1.3.66]).

  * n ∈ T 1 u n , we can write u * n = u * n,1 +u * n,2 with u * n,1 = Au n and u * n,2 ∈ γ * j ∂j(γ j u n ). Since A and γ * j ∂j(γ j •) are bounded operators, by passing to a subsequence if necessary, we may assume that u * n,1 ⇀u * a and u * n,2 ⇀u *

	b

Table 1

 1 Relative errors in energy norm for Example 1.

	h	1/4	1/8	1/16	1/32	1/64
	error	48.01%	20.77%	9.755%	4.745%	2.347%

Table 2

 2 Relative errors in energy norm for Example 3.

	h	1/4	1/8	1/16	1/32	1/64
	error	51.78%	23.17%	11.18%	5.443%	2.682%

Acknowledgment. We thank the two anonymous referees for their valuable comments on the first version of the paper.

The work of the authors was partially supported by the National Science Center of Poland under the Maestro 3 project DEC-2012/06/A/ST1/00262. The work of the first author was also supported by NSF grant DMS-1521684.

The element u is also a solution of Problem (P 1 ) due to the property (53). Uniqueness of the solution u can be shown similar to the argument in the uniqueness part of the proof of Theorem 3.1. In the rest of this subsection, we assume (55).

We introduce a finite element method to solve 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:

Let {T h } be a regular family of partitions of Ω into triangles/tetrahedrons that are compatible with the partition of the boundary ∂Ωi n t oΓ 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 .T h e n construct linear element spaces corresponding to T h :

The finite element approximation of Problem (P 1 ) is the following.

Under the same assumptions, Problem (P h 1 ) has a unique solution u h ∈ V h 1 .T o apply the theory developed in section 4, we note that [START_REF] Duvaut | Inequalities in Mechanics and Physics[END_REF]

1 defined by (56). The Lipschitz condition [START_REF] Eck | Unilateral Contact Problems: Variational Methods and Existence Theorems[END_REF] follows from (43)(a). Then we have the convergence by applying Theorem 4.2:

From [START_REF] Naniewicz | Mathematical Theory of Hemivariational Inequalities and Applications[END_REF], we have Céa's inequality

Then, applying finite element interpolation error estimates [START_REF] Atkinson | Theoretical Numerical Analysis: A Functional Analysis Framework[END_REF][START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF][START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF], we conclude the optimal order error bound

under the regularity assumptions

Note that we cannot derive optimal order error estimates from (59) for higher order elements. For example, if the conforming quadratic element is used in defining the finite element space V h 1 , then under higher solution regularity assumptions, we can only get a suboptimal error estimate

As a concrete example of j τ , let (cf. [START_REF] Barboteu | An analytical and numerical approach to a bilateral contact problem with nonmonotone friction[END_REF])

with λ 2,V > 0 the smallest eigenvalue of the eigenvalue problem

We have the following results. With the additional assumption (70)

Problem (P 2 ) has a unique solution u ∈ V . For the numerical solution u h defined by Problem (P h 2 ), we have Céa's inequality

From this inequality, we conclude the convergence of the method:

Moreover, under the regularity assumptions

we have the optimal order error bound

A frictionless unilateral contact problem. The third set of contact boundary conditions we consider is

These conditions model a frictionless contact with a foundation made of a rigid body covered by a layer made of elastic material. Penetration is restricted by the relation u ν ≤ g,whereg represents the thickness of the elastic layer. When there is penetration and the normal displacement does not reach the bound g, the contact is described by a multivalued normal compliance condition: -σ ν = ξ ν ∈ ∂j ν (u ν ). The problem described by ( 39)-( 42) and ( 74)-(75) represents the frictionless version of a nonlinear elastic contact model studied in [START_REF] Kalita | A class of subdifferential inclusions for elastic unilateral contact problems[END_REF], where the unique solvability of the model is provided by using both a primal and a dual variational formulation of the model. Here we continue the study of the model by providing the numerical analysis of the frictionless problem in [START_REF] Kalita | A class of subdifferential inclusions for elastic unilateral contact problems[END_REF]. To this end, for the potential function j ν :Γ 3 × R → R, we again assume (65).

Corresponding to the constraint u ν ≤ g on Γ 3 , we introduce a subset of the space V of (44):

The weak formulation of the contact problem is the following. Problem (P 3 ). Find a displacement field u ∈ U such that

To apply Theorem 3.1, we let

and use the operator A defined in (51) and the functional j defined in (69). Assuming additionally (70), Problem (P 3 ) has a unique solution u ∈ U . For its numerical solution, we use a related finite element subset of the space V h defined in (56):

(78)

Assume g is a concave function. Then U h ⊂ U . We define the following numerical method for Problem (P 3 ). Problem (P h 3 ). Find a displacement field u h ∈ U h such that

We apply [START_REF] Kalita | A class of subdifferential inclusions for elastic unilateral contact problems[END_REF] to derive an error estimate. The key step is to bound the residual term defined in [START_REF] Han | Quasistatic Contact Problems in Viscoelasticity and Viscoplasticity[END_REF]. We assume the regularity properties (72). Then,

Letting v be the sum of u and an arbitrary function from the subset Ũ ,w ed e r i v e from (77) that

With an argument similar to that in [26, section 8.1], we deduce from the above relation the following equalities:

Multiply (80) by vu with v ∈ U , integrate over Ω, and integrate by parts, Example 6.2. This is an example of Problem (P 2 ) . T h ebod yi sa s s u m e dt obei n frictionless contact with normal compliance. For simulation, we let (87)

-

ν on Γ 3 ,w i t hc 1 ν = 100, c 2 ν = -100, c 3 ν = 400, u 1 ν =0 .1, and u 2 ν =0 .15. In this case, -σ ν has a nonmonotone behavior with respect to the normal displacement u ν .

The deformed mesh and interface forces on Γ 3 corresponding to h =1 /64 are plotted in Figure 4. The zoom illustrates the nonmonotone relationship between the normal stress and the normal displacement on the contact zone. The average iteration number of the "convexification" procedure for the solution of discrete problems was 16 and the simulation was completed in 329 CPU time. Numerical solution error results are similar to that reported in Figure 3 for Example 1 and are omitted here. Example 6.3. This is an example of Problem (P 3 ). The contact boundary conditions on Γ 3 are characterized by a frictionless multivalued normal compliance contact in which the penetration is restricted by unilateral constraint. For simulations, we let u ν ≤ g, σ ν + ξ ν ≤ 0, (u νg)(σ ν + ξ ν )=0, (88)

σ τ = 0 (90) on Γ 3 . Note that in the conditions (88), g represents the maximum value of the allowed penetration. When this value of penetration is reached, the contact is unilateral; in