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

Tomasz Janiczko
email: tomasz.janiczko@gmail.com

NUMERICAL ANALYSIS OF A HYPERBOLIC HEMIVARIATIONAL INEQUALITY ARISING IN DYNAMIC CONTACT *

Keywords: hyperbolic hemivariational inequality, dynamic contact, linearly viscoelastic material, nonmonotone friction law, finite element method, error estimate

In this paper a fully dynamic viscoelastic contact problem is studied. The contact is assumed to be bilateral and frictional, where the friction law is described by a nonmonotone relation between the tangential stress and the tangential velocity. A weak formulation of the problem leads to a second order nonmonotone subdifferential inclusion, also known as a second order hyperbolic hemivariational inequality. We study both semidiscrete and fully discrete approximation schemes and bound the errors of the approximate solutions. Under some regularity assumptions imposed on the true solution, optimal order error estimates are derived for the linear element solution. This theoretical result is illustrated numerically.

1. Introduction. This paper provides error analysis for numerical methods to solve a hyperbolic hemivariational inequality arising in a dynamic bilateral contact process for a viscoelastic material. The main mathematical difficulty in the study of the problem is due to the nonmonotonicity of the friction law, and hence, we cannot apply the standard techniques based on convex analysis. We formulate the contact condition corresponding to the friction law by means of an inclusion involving the Clarke subdifferential of a locally Lipschitz potential. Consequently, we deal with a second order evolutionary hemivariational inequality as a starting point to the numerical analysis of the contact problem. For approximation of the hemivariational inequality, we discuss both the spatially semidiscrete and fully discrete schemes. We use the finite element method for the spatial discretization and backward difference to approximate the time derivative. In both cases we derive error estimates that are of optimal order when the linear elements are used, if the true solution has certain regularity. Finally we present results of computer simulations on a two-dimensional contact problem, to show the performance of the numerical methods and to provide numerical evidence of the theoretically predicted optimal convergence order of the linear element solutions.

The problem is on the cutting edge of contact mechanics, the theory, numerical analysis, and computer simulations of hemivariational inequalities. The mathematical modeling of contact problems in mechanics has reached a mature level, as is witnessed by the recent large number of publications on its theory and applications in engineering and industry. For details concerning classical contact models and their analysis, we refer to [START_REF] Han | Quasistatic Contact Problems in Viscoelasticity and Viscoplasticity[END_REF][START_REF] Kikuchi | Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods[END_REF][START_REF] Sofonea | Analysis and Approximation of Contact Problems with Adhesion or Damage[END_REF], where numerical analysis involving error estimates is also conducted in the case of quasi-static and dynamic problems. For more recent mathematical results devoted to contact mechanics we refer also to [START_REF] Sofonea | Variational Inequalities with Applications: A Study of Antiplane Frictional Contact Problems[END_REF][START_REF] Sofonea | Mathematical Models in Contact Mechanics[END_REF]. The theory of hemivariational inequalities, which allows us to model nonmonotone and nonsmooth contact problems, is a relatively new approach. Early comprehensive references in the area are [START_REF] Naniewicz | Mathematical Theory of Hemivariational Inequalities and Applications[END_REF][START_REF] Panagiotopoulos | Inequality Problems in Mechanics and Applications,B i r k h äuser[END_REF][START_REF] Panagiotopoulos | Hemivariational Inequalities: Applications in Mechanics and Engineering[END_REF][START_REF] Panagiotopoulos | Modelling of nonconvex nonsmooth energy problems: Dynamic hemivariational inequalities with impact effects[END_REF]. For more recent work, we refer to [START_REF] Migórski | Nonlinear Inclusions and Hemivariational Inequalities[END_REF] and the references cited there. In [START_REF] Haslinger | Finite Element Method for Hemivariational Inequalities: Theory, Methods and Applications[END_REF], the finite element method is studied for solving some hemivariational inequalities. There are, however, still few publications devoted to the error estimates in the numerical solution of hemivariational inequalities. In [START_REF] Barboteu | An analytical and numerical approach to a bilateral contact problem with nonmonotone friction[END_REF], numerical approximation for a static hemivariational inequality is studied. In [START_REF] Han | A class of variational-hemivariational inequalities with applications to frictional contact problems[END_REF], a class of variational-hemivariational inequalities is studied, theoretically and numerically. The numerical analysis presented here is also motivated by techniques used in [ 7 ,8 ,9 ,4 1 ,5 ] .

The rest of the paper is structured as follows. In section 2 we introduce the notation as well as some preliminary material. In section 3 we present the classical formulation of the frictional contact problem, list assumptions on the data, and present variational formulations of the problem. In section 4 we introduce and analyze a spatially semidiscrete scheme for solving the problem, and in section 5 we study a fully discrete approximation scheme. For both schemes, we derive optimal order error estimates for the linear element solutions under certain solution regularity assumptions. In section 6 we present numerical results in simulations of a two-dimensional contact problem and provide numerical evidence of optimal order convergence for the linear element solutions.

Notation and preliminaries.

In this section we present the notation and some preliminary material to be used later. For further details we refer the reader to [START_REF] Duvaut | Inequalities in Mechanics and Physics[END_REF][START_REF] Han | Quasistatic Contact Problems in Viscoelasticity and Viscoplasticity[END_REF][START_REF] Ionescu | Functional and Numerical Methods in Viscoplasticity[END_REF][START_REF] Panagiotopoulos | Inequality Problems in Mechanics and Applications,B i r k h äuser[END_REF].

We denote by S d the space of second order symmetric tensors on R d (d ≤ 3i n applications), and use " • "a n d" |•|" for the inner product and the Euclidean norm on R d and S d , respectively,

u • v = u i v i , |v| =(v • v) 1 2 ∀ u, v ∈ R d , σ • τ = σ ij τ ij , |τ | =(τ • τ ) 1 2 ∀ σ, τ ∈ S d .
Here and below the indices i and j run between 1 and d, and the summation convention over repeated indices is adopted.

Let Ω ⊂ R d be a bounded domain with a Lipschitz boundary Γ. The unit outward normal vector ν is defined a.e. on Γ. We introduce the following function spaces:

H = L 2 (Ω; R d )={u =(u i ) | u i ∈ L 2 (Ω)}, Q = {σ =(σ ij ) | σ ij = σ ji ∈ L 2 (Ω)}, H 1 = {u ∈ H | ε(u) ∈ Q}, Q 1 = {σ ∈ Q | Div σ ∈ H}.
Here ε: H 1 → Q and Div : Q 1 → H are the deformation and divergence operators, defined by

ε(u)=(ε ij (u)),ε ij (u)= 1 2 (u i,j + u j,i ), Div σ =(σ ij,j ),
respectively, where the index following a comma indicates the partial derivative with respect to the corresponding component of the independent variable. The spaces H, Q, H 1 ,a n dQ 1 are real Hilbert spaces endowed with the canonical inner products given by

(u, v) H = Ω u i v i dx,( σ, τ ) Q = Ω σ ij τ ij dx, (u, v) H 1 =(u, v) H +(ε(u), ε(v)) Q ,( σ, τ ) Q 1 =(σ, τ ) Q +(Divσ, Div τ ) H .
The associated norms on these spaces are denoted by

• H , • Q , • H 1 and • Q 1 , respectively. Let H Γ = H 1/2 (Γ; R d )a
n dl e tγ : H 1 → H Γ be the trace operator. For every element v ∈ H 1 ,w eu s et h es a m es y m b o lv to denote the trace γv of v on Γ, and we denote by v ν and v τ the normal and tangential components of v on the boundary Γ given by

v ν = v • ν, v τ = v -v ν ν.
Let H * Γ be the dual of H Γ and let

(σ, ε(v)) Q +(Divσ, v) H = σν, γv H * Γ ×H Γ ∀ v ∈ H 1 .
Moreover, if σ is a smooth (say, C 1 ) function, then

σν, γv H * Γ ×H Γ = Γ σν • v d Γ ∀ v ∈ H 1 .
We denote by σ ν and σ τ the normal and tangential traces of σ,

σ ν =(σν) • ν, σ τ = σν -σ ν ν.
Next recall the definitions of classical (one-sided) directional derivative and its generalization in the sense of Clarke. Let X be a Banach space and X * its dual. For a function ϕ : X → R,t h edirectional derivative of ϕ at x ∈ X in the direction v ∈ X is defined by

ϕ ′ (x; v) = lim λ↓0 ϕ(x + λv) -ϕ(x) λ
whenever this limit exists. The Clarke generalized directional derivative of a locally Lipschitz function ϕ : X → R at the point x ∈ X in the direction v ∈ X is defined by

ϕ 0 (x; v) = lim sup y→x,λ↓0 ϕ(y + λv) -ϕ(y) λ .
The Clarke subdifferential of ϕ at x is a subset of X * given by

∂ϕ(x)={ζ ∈ X * : ϕ 0 (x; v) ≥ ζ,v X * ×X ∀v ∈ X}.
A locally Lipschitz function ϕ : X → R is said to be regular (in the sense of Clarke) at x ∈ X if for all v ∈ X, the directional derivative ϕ ′ (x; v)existsandϕ ′ (x; v)=ϕ 0 (x; v). The function ϕ is regular (in the sense of Clarke) on X if it is regular at every point x ∈ X.

We will need the following discrete Gronwall inequality [ We assume that the body is clamped on Γ 1 and thus the displacement field vanishes there. A volume force of density f 0 acts in Ω and a surface traction of density f 2 acts on Γ 2 . The body is in frictional contact with an obstacle on Γ 3 . We assume the contact is bilateral, i.e., there is no loss of contact during the process. Thus, the normal displacement u ν vanishes on Γ 3 . We model the friction by a nonmonotone friction law. The dynamic process is considered.

The classical formulation of the mechanical problem is the following.

Problem P M . Find a displacement u :Ω× [0,T] → R d and a stress field σ : Ω × [0,T] → S d such that σ = Aε( u)+Bε(u)i n Ω × (0,T), (3.1) ρ ü =Divσ + f 0 in Ω × (0,T), (3.2) u = 0 on Γ 1 × (0,T), (3.3) σν = f 2 on Γ 2 × (0,T), (3.4) u ν =0 onΓ 3 × (0,T), (3.5) |σ τ |≤μ(0)S if uτ = 0, -σ τ = μ(| uτ |)S uτ | uτ | if uτ = 0 on Γ 3 × (0,T), (3.6) u(0) = u 0 , u(0) = u 1 in Ω. (3.7)
Here, (3.1) is the linearly viscoelastic constitutive law [START_REF] Duvaut | Inequalities in Mechanics and Physics[END_REF][START_REF] Han | Quasistatic Contact Problems in Viscoelasticity and Viscoplasticity[END_REF], (3.2) is the equation of motion, where ρ is the mass density, (3.3) is the homogeneous displacement boundary condition on Γ 1 , (3.4) is the traction boundary condition on Γ 2 ,( 3 . 5 )r e presents the bilateral contact condition, and (3.7) provides the initial displacement and velocity conditions. In (3.6), μ(| uτ |)S represents the magnitude of the limiting friction traction at which slip begins, S ≥ 0 being given. The friction coefficient μ is allowed to depend on the tangential speed | uτ |. The strict inequality in (3.6) holds in the stick zone and the equality holds in the slip zone. This physical model of slip-dependent friction was introduced in [START_REF] Rabinowicz | The nature of the static and kinetic coefficients of friction[END_REF] for geophysical context of earthquake modeling and it also was studied in [START_REF] Ionescu | Dynamic contact problems with slip dependent friction in viscoelasticity[END_REF][START_REF] Ionescu | Slip-dependent friction in dynamic elasticity[END_REF][START_REF] Ionescu | On the contact problem with slip displacement dependent friction in elastostatics[END_REF][START_REF] Kuttler | Dynamic contact with normal compliance wear and discontinuous friction coefficient[END_REF][START_REF] Kuttler | Dynamic contact with Signorini's condition and slip rate dependent[END_REF][START_REF] Migórski | Hemivariational inequality for viscoelastic contact problem with slip-dependent friction[END_REF][START_REF] Shillor | Models and Analysis of Quasistatic Contact[END_REF].

In the study of the contact problem we need the following assumptions on its data: (e) (A(x,t,ε 1 ) -A(x,t,ε 2 )) :

H(A): The viscosity operator A :Ω× [0,T] × S d → S d satisfies ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (a) A(•, •, ε)i sm e a s u r a b l eo nΩ× [0,T] ∀ ε ∈ S d ; (b) A(x,
(ε 1 -ε 2 ) ≥ m A |ε 1 -ε 2 | 2 ∀ ε 1 , ε 2 ∈ S d , a.e. (x,t) ∈ Ω × [0,T] with m A > 0; (f) |A(x,t,ε 1 ) -A(x,t,ε 2 )|≤L A |ε 1 -ε 2 | ∀ ε 1 , ε 2 ∈ S d , a.e. (x,t) ∈ Ω × [0,T] with L A > 0.

H(B):

The elasticity operator B :Ω× S d → S d is a bounded, symmetric, nonnegatively definite fourth order tensor, i.e.,

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ (a) B ijkl ∈ L ∞ (Ω), 1 ≤ i, j, k, l ≤ d; (b) Bσ • τ = σ •Bτ ∀ σ, τ ∈ S d , a.e. in Ω; (c) Bτ • τ ≥ 0 ∀ τ ∈ S d , a.e. in Ω.

H(f ):

The force and the traction densities satisfy

f 0 ∈ L 2 (0,T; L 2 (Ω; R d )), f 2 ∈ L 2 (0,T; L 2 (Γ 2 ; R d )). H(μ): The friction bound μ :[0, ∞) → R satisfies ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ (a) μ is continuous; (b) |μ(s)|≤c(1 + s) ∀ s ≥ 0,c>0; (c) μ(s 1 ) -μ(s 2 ) ≥-λ(s 1 -s 2 ) ∀ s 1 >s 2 ≥ 0 with λ>0. Remark 3.1. If A(•, •, ε) is linear in ε,t h e nH(A)(d)a n dH(A)(e)a r ee q u i v a l e n t with α = m A ,a n dH(A)(f ) implies H(A)(c) with a 0 =0anda 1 = L A .
Since μ corresponds to the physical resistance force, it is nonnegative. However, in mathematical analysis of the contact problem, we do not need to impose this condition. The condition (c) is the so-called one-side Lipschitz condition, which allows the function to decrease at a rate not faster than λ.

Using the Clarke subdifferential (cf. [START_REF] Clarke | Optimization and Nonsmooth Analysis[END_REF]), we can express the friction condition (3.6) in another form. Indeed, define a function j : R d → R by

(3.8) j(ξ)=S |ξ| 0 μ(s) ds ∀ ξ ∈ R d .
Then, assuming H(μ)(a) -(b), the condition (3.6) is equivalent to the following subdifferential inclusion:

-σ τ ∈ ∂j( uτ )o n Γ 3 × (0,T),
where ∂j(ξ) denotes the Clarke subdifferential of j at the point ξ ∈ R d .

Properties of the function j are summarized in the next lemma. Lemma 3.2. If the assumptions H(μ)(a)-(b) hold, then the function j defined by (3.8) is regular in the sense of Clarke, it is locally Lipschitz, and

|η|≤Sc(1 + |ξ|) ∀ ξ ∈ R d , η ∈ ∂j(ξ).
If furthermore the assumption H(μ)(c) holds, then we have

(3.9) (η 1 -η 2 ) • (ξ 1 -ξ 2 ) ≥-Sλ|ξ 1 -ξ 2 | 2 ∀ ξ 1 , ξ 2 ∈ R d , η i ∈ ∂j(ξ i ),i=1, 2.
Proof. We will show that j is regular in the sense of Clarke. First observe that for ξ = 0 we have ∂j(ξ)={Sμ(|ξ|)ξ/|ξ|} and so j is regular at ξ [12, Proposition 5.6.15]. Next, consider the case ξ = 0.L e tv ∈ R d .U s i n gH(μ)(a)w eh a v e

j ′ (0; v) = lim λ↓0 1 λ S |λv| 0 μ(t) dt = Sμ(0) |v|.
By definition,

j 0 (0; v) = lim sup ξ→0,λ↓0 S λ |ξ+λv| |ξ| μ(s) ds. Since μ ∈ C([0, ∞)), j 0 (0; v)=Sμ(0) lim sup ξ→0,λ↓0 |ξ + λv|-|ξ| λ = Sμ(0) lim sup ξ→0,λ↓0 2 ξ • v + λ |v| 2 |ξ + λv| + |ξ| = Sμ(0) lim sup ξ→0 ξ |ξ| • v = Sμ(0) |v|.
So j is regular at 0. The other properties then follow straightforwardly.

To introduce a weak formulation of the mechanical problem P M , we first define a closed subspace of H 1 ,

V = {v ∈ H 1 | v = 0 on Γ 1 ,v ν =0onΓ 3 }.
Since meas (Γ 1 ) > 0, Korn's inequality holds [34, p. 79]: for some constant C K > 0, depending only on Ω and Γ 1 ,

(3.10) ε(v) Q ≥ C K v H 1 ∀ v ∈ V.
On V , we use the inner product given by

(3.11) (u, v) V =(ε(u), ε(v)) Q ∀ u, v ∈ V
and let • V be the associated norm, i.e., (3.12)

v V = ε(v) Q ∀ v ∈ V.
It follows from (3.10) and (3.12) that • H 1 and • V are equivalent norms on V and therefore (V, • V ) is a real Hilbert space. The duality pairing between V and V * is denoted by •, • . Identifying H with its dual, we have an evolution triple V ⊂ H ⊂ V * with dense, continuous, and compact embeddings. We denote by i : V → H the identity mapping and by i * : V * → H its adjoint mapping. By the Sobolev trace theorem and by (3.10) there exists a constant C 0 depending only on the domain Ω, Γ 1 ,a n dΓ 3 such that

(3.13) v L 2 (Γ 3 ) d ≤ C 0 v V ∀ v ∈ V.
By (3.13) there exists a continuous trace operator γ : V → L 2 (Γ 3 ; R d ) and for the function v ∈ V we still denote by v its trace γv. In what follows we need the spaces V = L 2 (0,T; V ), H = L 2 (0,T; H), and W = {v ∈V| v ∈V * }, where the time derivative involved in the definition of W is understood in the sense of vector valued distributions. Equipped with the norm

v W = v 2 V + v 2 V *
1/2 the space W becomes a separable Hilbert space. We also have W⊂V⊂H⊂V * .I ti sw e l lk n o w n that the embeddings W⊂C([0,T]; H)a n d{w ∈V| ẇ ∈W }⊂C([0,T]; V )a r e continuous. Next we define operators A :(0,T)

× V → V * and B : V → V * by A(t, u), v =(A(t, ε(u)), ε(v)) Q for u, v ∈ V and t ∈ (0,T), (3.14) Bu, v =(Bε(u), ε(v)) Q for u, v ∈ V, (3.15) a functional J : L 2 (Γ 3 ; R d ) → R by (3.16) J(v)= Γ 3 j(v) dΓf o r v ∈ L 2 (Γ 3 ; R d ),
and a function f :(0,T) → V * by (3.17)

f (t), v = Ω f 0 (t) • v dx + Γ 2 f 2 (t) • v dΓf o r v ∈ V, a.e.t∈ (0,T).
Assuming H(A), we have the following properties for the operator

A :[ 0 ,T]× V → V * : ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (a) A(•, v)i sm e a s u r a b l eo n( 0 ,T) ∀v ∈ V ; (b) A(t, •) is pseudomonotone on V for a.e. t ∈ (0,T); (c) A(t, v) V * ≤ a 0 (t)+a 1 v V ∀v ∈ V, a.e. t ∈ (0,T) with a 0 ∈ L 2 (0,T),a 0 ≥ 0, and a 1 > 0; (d) A(t, v), v ≥α v 2 V ∀v ∈ V, a.e. t ∈ (0,T) with α>0; (e) A(t, v 1 ) -A(t, v 2 ), v 1 -v 2 ≥m A v 1 -v 2 2 V ∀ v 1 , v 2 ∈ V, a.e. t ∈ (0,T) with m A = m A > 0; (f) A(t, v 1 ) -A(t, v 2 ) V * ≤ L A v 1 -v 2 V ∀ v 1 , v 2 ∈ V, a.e. t ∈ (0,T) with L A = L A > 0.
Under the assumption H(B), the operator B ∈ L(V, V * ) is self-adjoint and monotone. Under the assumption H(μ), the functional J : L 2 (Γ 3 ; R d ) → R is locally Lipschitz, and we have the following inequalities:

η L 2 (Γ 3 ;R d ) ≤ SC(1 + v L 2 (Γ 3 ;R d ) ) ∀ η ∈ ∂J(v), (3.18) η 1 -η 2 , v 1 -v 2 L 2 (Γ 3 ;R d ) ≥-Sλ v 1 -v 2 2 L 2 (Γ 3 ;R d ) ∀ η i ∈ ∂J(v i ),i=1, 2, (3.19) where C = √ 2 c max 1, meas d-1 (Γ 3 ) . The assumption H(f ) implies f ∈V * .
For the initial values, we will assume the following:

H 0 : u 0 ∈ V , u 1 ∈ H.
Proceeding in a standard way [START_REF] Han | Quasistatic Contact Problems in Viscoelasticity and Viscoplasticity[END_REF][START_REF] Migórski | Nonlinear Inclusions and Hemivariational Inequalities[END_REF], we obtain the following variational formulation of the frictional Problem P M .

Problem P V . Find a displacement field u ∈Vwith u ∈Wand a friction density

ξ τ ∈ L 2 (0,T; L 2 (Γ 3 ; R d )) such that ρ ü(t)+A(t, u(t)) + Bu(t) -f (t), v = Γ 3 ξ τ (t)•v τ dΓ ∀ v ∈ V, a.e. t, (3.20) -ξ τ ∈ ∂j( uτ )a . e . o n Γ 3 × (0,T), (3.21) u(0) = u 0 , u(0) = u 1 .
Here and below, "a.e. t" means "a.e. t ∈ (0,T)." The above problem can be expressed equivalently as follows.

Problem P V,1 . Find a displacement field u ∈V with u ∈W such that

ρ ü(t)+A(t, u(t)) + Bu(t) -f (t), v + Γ 3 j 0 ( uτ (t); v τ ) dΓ ≥ 0 ∀ v ∈ V, a.e. t, u(0) = u 0 , u(0) = u 1 . Problem P V,
u(0) = u 0 , u(0) = u 1 ,
where γ is the trace operator on Γ 3 and γ * its adjoint, and γ uτ means (γ u) τ .

A function u ∈V is a solution of Problem P V,2 if and only if u ∈W and there exists η ∈ L 2 (0,T;

L 2 (Γ 3 ; R d )) such that ρ ü(t)+A(t, u(t)) + Bu(t)+η(t)=f (t) a.e. t ∈ (0,T), η(t) ∈ γ * ∂J(t, γ uτ (t)) a.e. t ∈ (0,T), u(0) = u 0 , u(0) = u 1 .
The hemivariational inequality corresponding to Problem P V,2 reads as follows.

Problem P V,3 . Find a displacement field u ∈V with u ∈W such that ρ ü(t)+A(t, u(t)) + Bu(t) -f (t), v + J 0 (γ uτ (t); γv τ ) ≥ 0 ∀ v ∈ V, a.e. t, u(0) = u 0 , u(0) = u 1 .
We complete this section with a result on solution existence and uniqueness for Problem P V,2 .

Theorem 3.3. Assume H(A), H(B), H(μ), H(f ), H 0 ,a n d

(3.22) α 2 >SCC 2 0 ,m A >SλC 2 0 .
Then Problem P V,2 has a unique solution u, and the following bound holds:

(3.23) u C(0,T ;V ) + u W ≤ C (1 + u 0 V + u 1 H + f V * )
with a positive constant C.

The proof of this result follows from the arguments used in the proof of Theorem 5.15 (for existence and uniqueness of a solution) and of Lemma 5.8 (for the bound (3.23)) of [START_REF] Migórski | Nonlinear Inclusions and Hemivariational Inequalities[END_REF].

Since j is regular (cf. Lemma 3.2), Problems P V , P V,1 , P V,2 ,a n dP V,3 are equivalent [START_REF] Migórski | Boundary hemivariational inequalities of hyperbolic type and applications[END_REF]Remark 4]. In particular, under the assumptions stated in Theorem 3.3, Problem P V has a unique solution.

4. Spatially semidiscrete approximation. In this section we introduce and analyze a spatially semidiscrete approximation for Problem P V . Let V h be a finite dimensional subspace of V ,w h e r eh>0 denotes a spatial discretization parameter. Let u h 0 , u h 1 ∈ V h be suitable approximations of u 0 and u 1 , characterized by

(4.1) (u h 0 -u 0 , v h ) V =0, (u h 1 -u 1 , v h ) H =0 ∀ v h ∈ V h .
It is easy to observe that (4.2)

u h 0 V ≤ u 0 and u h 1 H ≤ u 1 H .
Then we have the following semidiscrete approximation of Problem P V .

Problem P h V . Find a displacement field u h ∈ L 2 (0,T; V h ) with uh , üh ∈ L 2 (0,T; V h ), and a friction density

ξ h τ ∈ L 2 (0,T; L 2 (Γ 3 ; R d )) such that ρ üh (t)+A(t, uh (t)) + Bu h (t) -f (t), v h (4.3) = Γ 3 ξ h τ (t) • v h τ dΓ ∀ v h ∈ V h , a.e. t, -ξ h τ (t) ∈ ∂j( uh τ (t)) a.e. on Γ 3 × (0,T), (4.4) u h (0) = u h 0 , uh (0) = u h 1 . (4.5)
Under the assumptions of Theorem 3.3, we also have the existence and uniqueness of a solution to Problem P h V . Moreover, similar to (3.23), and thanks to (4.2), we have the bound

(4.6) u h C(0,T ;V ) + uh W ≤ C (1 + u 0 V + u 1 H + f V * )
with a positive constant C.

We provide a result on the error estimates between the solutions of Problems P V and P h V . Theorem 4.1. Assume that H(A), H(B), H(μ), H(f ), H 0 ,a n d (3.22) hold. Let u and u h be solutions of Problems P V and P h V , respectively. Then there exists a positive constant c depending only on the data of the problem, such that for any

v h ∈ L 2 (0,T; V h ) ∩W, u -u h 2 C(0,T ;V ) + u -uh 2 C(0,T ;H) + u -uh 2 V (4.7) ≤ c u 0 -u h 0 2 V + u 1 -u h 1 H u 1 -v h (0) H + u -v h 2 V + ü -vh 2 V * + uτ -v h τ L 2 (0,T ;L 2 (Γ 3 ;R d )) .
Proof. Let us define the functions w(t)= u(t)andw h (t)= uh (t) for all t ∈ [0,T]. Then,

u(t)=(Iw)(t)=u 0 + t 0 w(s) ds, (4.8) u h (t)=(I h w h )(t)=u h 0 + t 0 w h (s) ds,
ρ ẇ(t)+A(t, w(t)) + B(Iw)(t) -f (t), v (4.9) = Γ 3 ξ τ (t)•v τ d Γ ∀ v ∈ V, a.e. t,
ξ τ ∈ ∂j(w τ )a . e .o nΓ 3 × (0,T), (4.10)

w(0) = u 1 , (4.11) ρ ẇh (t)+A(t, w h (t)) + B(I h w h )(t) -f (t), v h (4.12) = Γ 3 ξ h τ (t)•v h τ d Γ ∀ v h ∈ V h , a.e. t,
ξ h τ ∈ ∂j(w h τ )a . e .o nΓ 3 × (0,T), (4.13)

w h (0) = u h 1 . (4.14)
For any v h ∈ V h , we have from (4.9) and (4.12) that for a.e. t ∈ (0,T),

ρ ẇ(t) -ẇh (t), v h + A(t, w(t)) -A(t, w h (t)), v h (4.15) + B(Iw)(t) -B(I h w h )(t), v h + Γ 3 (ξ h τ (t) -ξ τ (t)) • v h τ dΓ=0. Note that (4.16) ẇ(t) -ẇh (t), w(t) -w h (t) = 1 2 d dt w(t) -w h (t) 2 H .
From the strong monotonicity of A,w eh a v e (4.17)

m A w(t) -w h (t) 2 V ≤ Aw(t) -Aw h (t), w(t) -w h (t) .
By the symmetry of B,w eh a v e

B(Iw)(t) -B(I h w h )(t), w(t) -w h (t) (4.18) = Bu(t) -Bu h (t), u(t) -uh (t) = 1 2 d dt Bu(t) -Bu h (t), u(t) -u h (t) .
From (4.10), (4.13), (3.9), and (3.13), we obtain (4.19)

Γ 3 (ξ h τ (t) -ξ τ (t)) • (w τ (t) -w h τ (t)) dΓ ≤ SλC 2 0 w(t) -w h (t) 2 V . Denote c 0 = m A -SλC 2 0 . From (4.16)-(4.19) we obtain (4.20) 1 2 ρ d dt w(t) -w h (t) 2 H + c 0 w(t) -w h (t) 2 V + 1 2 d dt Bu(t) -Bu h (t), u(t) -u h (t) ≤ ρ ẇ(t) -ẇh (t), w(t) -w h (t) + A(t, w(t)) -A(t, w h (t)), w(t) -w h (t) + B(Iw)(t) -B(I h w h )(t), w(t) -w h (t) + Γ 3 (ξ h τ (t) -ξ τ (t)) • (w τ (t) -w h τ (t)) dΓ = ρ ẇ(t) -ẇh (t), w(t) -v h (t) + A(t, w(t)) -A(t, w h (t)), w(t) -v h (t) + B(Iw)(t) -B(I h w h )(t), w(t) -v h (t) + Γ 3 (ξ h τ (t) -ξ τ (t)) • (w τ (t) -v h τ (t)) dΓ,
where v h (t) ∈ V h for a.e. t ∈ (0,T) is arbitrary and the last equality follows from (4.15). For t ∈ (0,T), assuming v h ∈W we perform integration by parts [13, Proposition 8.4.14]:

t 0 ẇ(s) -ẇh (s), w(s) -v h (s) ds= w(t) -w h (t), w(t) -v h (t) H -w(0) -w h (0), w(0) -v h (0) H - t 0 w(s) -w h (s), ẇ(s) -vh (s) ds.
Thus,

t 0 ẇ(s) -ẇh (s), w(s) -v h (s) ds (4.21) ≤ w(t) -w h (t) H w(t) -v h (t) H + u 1 -u h 1 H u 1 -v h (0) H + t 0 w(s) -w h (s), ẇ(s) -vh (s) ds ≤ 1 4 w(t) -w h (t) 2 H + w(t) -v h (t) 2 H + u 1 -u h 1 H u 1 -v h (0) H + ε w -w h 2 V + 1 4ε ẇ -vh 2 V * .
Using the Lipschitz continuity of A,

t 0 A(s, w(s)) -A(s, w h (s)), w(s) -v h (s) ds (4.22) ≤ t 0 L A w(s) -w h (s) V w(s) -v h (s) V ds ≤ ε w -w h 2 V + L A 4ε w -v h 2 V .
Using the properties of B,

t 0 B(Iw)(s) -B(I h w h )(s), w(s) -v h (s) (4.23) ≤ t 0 B L(V,V * ) u(s) -u h (s) V w(s) -v h (s) V ≤ ε u -u h 2 V + B L(V,V * ) 4ε w -v h 2 V ≤ ε2T u 0 -u h 0 2 V + ε2T w -w h 2 V + B L(V,V * ) 4ε w -v h 2 V .
It remains to bound the last term of (4.20). From (3.18) and (3.23), we have (4.24)

t 0 Γ 3 (ξ h τ (s) -ξ τ (s)) • (w τ (s) -v h τ (s)) dΓ ds ≤ t 0 ξ h τ (s) L 2 (Γ 3 ;R d ) + ξ τ (s) L 2 (Γ 3 ;R d ) w τ (s) -v h τ (s) L 2 (Γ 3 ;R d ) ds ≤ t 0 SC 2+C 0 ( w(t) V + w h (t) V ) w τ (s) -v h τ (s) L 2 (Γ 3 ;R d ) ds ≤ 2SC √ T + C 0 ( w V + w h V ) w τ -v h τ L 2 (0,T ;L 2 (Γ 3 ;R d )) ≤ 2SC √ T +2 C(1 + u 0 V + u 1 H + f V * ) w τ -v h τ L 2 (0,T ;L 2 (Γ 3 ;R d )) .
Denote

r= u 0 -u h 0 2 V + u 1 -u h 1 H u 1 -v h (0) H + w -v h 2 V + ẇ -vh 2 V * + w τ -v h τ L 2 (0,T ;L 2 (Γ 3 ;R d )) .
We integrate (4.20) and apply (4.21)-(4.24) to get

(4.25) 1 2 ρ w(t) -w h (t) 2 H +2(c 0 -(2 + 2T )ε) t 0 w(s) -w h (s) 2 V ds ≤ c 1 r,
where the constant c 1 depends only on the data of the problem. Since t ∈ (0,T)i s arbitrary, with ε small enough, we obtain from (4.25) that

(4.26) w -w h 2 C(0,T ;H) + w -w h 2 V ≤ c 2 r
with c 2 > 0. For any t ∈ (0,T),

u(t) -u h (t) V ≤ u 0 -u h 0 V + T 0 w(t) -w h (t) V dt ≤ u 0 -u h 0 V + √ T w -w h V .
Thus,

(4.27) u -u h 2 C(0,T ;V ) ≤ 2 u 0 -u h 0 2 +2 √ T w -w h 2 V ≤ c 3 r
with a positive constant c 3 . From (4.26)-(4.27) we obtain the result (4.7). Remark 4.2. It follows from (3.23) and (4.6) that u, u h ∈ C(0,T; V ). Since the embedding W⊂C(0,T; H) is continuous, again from (3.23) and (4.6), we have u, uh ∈ C(0,T; H). Thus it is reasonable to estimate the norms uu h C(0,T ;V ) and uuh C(0,T ;H) in (4.7). Theorem 4.1 is valid for any finite dimensional subspace V h of V . In applications, V h is usually taken to be a finite element space. As a particular example, assume Ω is a polygonal/polyhedral domain and let {T h } be a regular family of finite element triangulations of Ω into triangles (d = 2) or tetrahedrons (d =3 ) . F o ra ne l e m e n t T ∈T h ,d en ot eb yP 1 (T ) the space of polynomials of a total degree less than or equal to one in T . Then we can use the linear element space of continuous piecewise affine functions:

(4.28) V h = {v h ∈ [C(Ω)] d : v h | T ∈ [P 1 (T )] d ∀ T ∈T h , v h = 0 on Γ 1 ,v h ν =0onΓ 3 }.
In the numerical simulations presented in section 6, this linear element space with d = 2 is used. Corollary 4.3. Keep the assumptions stated in Theorem 4.1. Assume Ω is a polygonal/polyhedral domain, and let {V h } be the family of linear element spaces defined by (4.28), corresponding to a regular family of finite element triangulations of Ω into triangles or tetrahedrons. Let u and u h be solutions of Problems P V and P h V , respectively. Assume

u 0 ∈ H 2 (Ω; R d ), u 1 ∈ H 1 (Ω; R d ),a n dt a k eu h 0 , u h 1 ∈ V h
to be projections of u 0 and u 1 , defined by (4.1). Under the regularity condition

u ∈ L 2 (0,T; H 2 (Ω; R d )), ü ∈ L 2 (0,T; H 2 (Ω; R d )), uτ ∈ L 2 (0,T; H 2 (Γ 3 ; R d )),
we have the optimal order error estimate

(4.29) u -u h C(0,T ;V ) + u -uh C(0,T ;H) + u -uh V ≤ ch for a constant c independent of h.
Proof. Note that under the stated regularity assumptions, for a.e. t ∈ [0,T], u(t), ü(t) are continuous on Ω, and uτ (t) is continuous on Γ 3 .L e tv h (t)=Π h u(t) ∈ V h be the finite element interpolant of u(t), a.e. t ∈ [0,T]. Note that v h τ (t)=(Π h u(t)) τ is the continuous piecewise linear interpolant of uτ (t)o nΓ 3 .M o r e o v e r ,v h (t)i st h e continuous piecewise linear interpolant of ü(t). Then by the standard 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 have the following approximation properties:

u(t) -v h (t) V ≤ ch u(t) H 2 (Ω;R d ) , ü(t) -vh (t) V * ≤ ch ü(t) H 2 (Ω;R d ) , uτ (t) -v h τ (t) L 2 (Γ 3 ;R d ) ≤ ch 2 uτ H 2 (Γ 3 ;R d )
and

u 0 -u h 0 V ≤ ch u 0 H 2 (Ω;R d ) , u 1 -u h 1 H ≤ ch u 1 H 1 (Ω;R d ) . It follows that u -v h V ≤ ch u L 2 (0,T ;H 2 (Ω;R d )) , ü -vh V * ≤ ch ü L 2 (0,T ;H 2 (Ω;R d )) , uτ -v h τ L 2 (0,T ;L 2 (Γ 3 ;R d )) ≤ ch 2 uτ L 2 (0,T ;H 2 (Γ 3 ;R d ))
. Then the error bound (4.29) follows from (4.7).

Note that for other choices of the finite element space V h , Theorem 4.1 can be applied similarly to derive error estimates of the finite element solutions, under certain corresponding regularity assumptions on the true solution u.

Fully discrete error estimates.

In this section we introduce a fully discrete approximation of Problem P V and bound the error of the fully discrete solutions. For simplicity in exposition, we assume (5.1)

A(•, v) ∈ C(0,T; V * ) ∀ v ∈ V, f ∈ C(0,T; V * ).
In addition to the finite dimensional subspace V h ⊂ V for spatial discretization, we need temporal discretization. We define a uniform partition of [0,T] denoted by 0 = t 0 <t 1 < ••• <t N = T .L e t k = T/N be a time step size and for a continuous function g we denote g n = g(t n ). For a sequence {z n } N n=0 ,w ed e n o t eb y δz n =( z n -z n-1 )/k for n =1 ,...,N the backward divided difference. With the backward Euler scheme for the time derivative, the fully discrete approximation of the Problem P V is the following.

Problem P kh V . Find a velocity field {w hk n } N n=0 ⊂ V h and a friction density

{ξ hk n } N n=0 ⊂ L 2 (Γ 3 ; R d ) such that ρδw hk n + A(t n , w hk n )+Bu hk n -f n , v h = Γ 3 ξ hk n τ • v h τ d Γ ∀ v h ∈ V h , (5.2)
-ξ hk n τ ∈ ∂j(w hk n τ ) a.e. on Γ 3 ,n=1,...,N, (

and

(5.4)

w hk 0 = u h 1 , where the discrete displacement field {u hk n } N n=0 ⊂ V h is given by (5.5) u hk n = u h 0 + n j=1 kw hk j .
Under the assumptions of Theorem 3.3, there exists a unique solution of Problem P kh V . The following boundedness property on the numerical solution will be needed in error estimation.

Theorem 5.1. Assume H(A), H(B), H(μ), H(f ), H 0 ,a n d (3.22). Then for some constant C>0,

(5.6) k N n=1 w hk n 2 ≤ C.
Proof. Taking v h = w hk n in (5.2) and using (5.5) we have

ρ(w hk n -w hk n-1 , w hk n ) H + k A(t n , w hk n ), w hk n + Bu hk n , u hk n -u hk n-1 (5.7)
= k f n , w hk n + k Γ 3 ξ hk nτ • w hk nτ d Γ.
Note that

(w hk n -w hk n-1 , w hk n ) H = 1 2 w hk n 2 H - 1 2 w hk n-1 2 H + 1 2 w hk n -w hk n-1 2
H .
From the property (d) of the operator A,w eg e t (5.8)

A(t n , w hk n ), w hk n ≥α w hk n 2 V .
From the properties of B,w eo b t a i n

Bu hk n , u hk n -u hk n-1 = 1 2 Bu hk n , u hk n - 1 2 Bu hk n-1 , u hk n-1
(5.9)

+ 1 2 B(u hk n -u hk n-1 ), u hk n -u hk n-1 ≥ 1 2 Bu hk n , u hk n - 1 2 Bu hk n-1 , u hk n-1 . Moreover, (5.10)
f n , w hk n ≤ f n V * w hk n V ≤ α 4 w hk n 2 V + 1 α f n 2 V * .
From (5.3) and (3.18), we get (5.11)

Γ 3 ξ hk nτ • w hk nτ d Γ ≤ 1 α SCC 2 0 + SCC 2 0 + α 4 w hk n 2 V .
Using (5.7)-(5.11) we obtain, with

c 0 = α/2 -SCC 2 0 , ρ 1 2 w hk n 2 H + ρ 1 2 w hk n -w hk n-1 2
H + c 0 k w hk n 2 V + 1 2 Bu hk n , u hk n ≤ k 1 α f n 2 V * + k 1 α SCC 2 0 + ρ 1 2 w hk n-1 2 H + 1 2 Bu hk n-1 , u hk n-1 .
Summing up the last inequality for n =1,...,N we obtain

ρ 1 2 w hk N 2 H + ρ 1 2 N n=1 w hk n -w hk n-1 2
H + c 0 k N n=1 w hk n 2 V ≤ 1 α f 2 V * + T 1 α SCC 2 0 + ρ 1 2 u h 1 2 H + 1 2 Bu h 0 , u h 0 .
From the last inequality and (3.22) we obtain (5.6). Now we state a result on error estimation. Theorem 5.2. Assume H(A), H(B), H(μ), H(f ), H 0 ,a n d (3.22),a n df o rt h e solution u of Problem P V ,

(5.12) u ∈ C 2 (0,T; H) ∩ C 1 (0,T; V ), uτ ∈ C(0,T; L 2 (Γ 3 ; R d )).
From (3.18) and (3.13) we also have (5.18)

Γ 3 (ξ hk n τ -ξ nτ )•(w nτ -v h τ ) dΓ ≤ C(1+ w n V + w hk n V ) w nτ -v h τ L 2 (Γ 3 ;R d ) .
In further estimations we use (5.16), strong monotonicity of A and (3.19) for the left-hand side of (5.15) and (5.17), (5.18), the properties of B, the inequalities ab ≤ 2a 2 +2b 2 or ab ≤ ǫa 2 + b 2 /4ǫ,f o ra, b, ǫ > 0, and (3.23) for its right-hand side. Thus we get, with

c 0 = m A -SλC 2 0 > 0a n dǫ<c 0 , 1 2k w n -w hk n 2 H -w n-1 -w hk n-1 2
H + c 0 w n -w hk n 2 V (5.19) ≤ C ẇn -δw n 2 H + w n -v h 2 V + u n -u hk n 2 V +(1+ w n V + w hk n V ) w nτ -v h τ L 2 (Γ 3 ;R d ) + ǫ w n -w hk n 2 V + ρ(δw n -δw hk n ), w n -v h H .
We replace n by j in the relation (5.19) and sum over j from 1 to n to obtain (5.20)

w n -w hk n 2 H +2k (c 0 -ǫ) n j=1 w j -w hk j 2 V ≤ w 0 -w hk 0 2 V + Ck n j=1 ẇj -δw j 2 H + w j -v h j 2 V + u j -u hk j 2 V + Ck n j=1 (1 + w j V + w hk j V ) w jτ -v h jτ L 2 (Γ 3 ;R d ) +2k n j=1 ρ(δw j -δw hk j ), w j -v h j H
for all {v h j } n j=1 ⊂ V h . We also have

n j=1 k ρ(δw j -δw hk j ), w j -v h j H (5.21) = n j=1 ρ(w j -w hk j -(w j-1 -w hk j-1 )), w j -v h j H ≤ ǫ w n -w hk n 2 H + C w n -v h n 2 H + c w 0 -w h 0 2 H + c w 1 -v h 1 2 H + n-1 j=1 ρ w j -w hk j H w j -v h j -(w j+1 -v h j+1 ) H ≤ ǫ w n -w hk n 2 H + C w n -v h n 2 H + c w 0 -w h 0 2 H + c w 1 -v h 1 2 H + n-1 j=1 4ρk w j -w hk j 2 H + 1 k n-1 j=1 w j -v h j -(w j+1 -v h j+1 ) 2 H .
Taking (4.8) at time t = t j and subtracting it from (5.5) we find that

(5.22) u j -u hk j V ≤ u 0 -u h 0 V + j l=1 k w l -w hk l V + I j ,
where I j is the integration error given by

I j = t j 0 w(s) ds - j l=1 kw l V .
We know that [START_REF] Han | Quasistatic Contact Problems in Viscoelasticity and Viscoplasticity[END_REF] I j ≤ k u H 2 (0,T ;V ) .

From (5.22) we get

u j -u hk j 2 V ≤ C u 0 -u h 0 2 V + j j l=1 k 2 w l -w hk l 2 V + k 2 u 2
H 2 (0,T ;V ) ; using inequality j ≤ n ≤ N and the fact that Nk = T we estimate

n j=1 k u j -u hk j 2 V ≤ CT u 0 -u h 0 2 V + k 2 u 2 H 2 (0,T ;V ) (5.23) + T n j=1 k j l=1 w l -w hk l 2 V .
Denote e n := w nw hk n 2

H + n j=1 k w j -w hk j 2
V and

g n := w 0 -w hk 0 2 V + k n j=1 ẇj -δw j 2 H + w j -v h j 2 V + Ck n j=1 1+ w j V + w hk j V w jτ -v h jτ L 2 (Γ 3 ;R d ) + u 0 -u h 0 2 V + k 2 u 2 H 2 (0,T ;V ) + w n -v h n 2 H + w 0 -w h 0 2 H + w 1 -v h 1 2 H + 1 k n-1 j=1 w j -v h j -(w j+1 -v h j+1 ) 2 H .
Then, from (5.20), (5. (5.25)

≤ T + T w C(0.T ;V ) + √ T k N j=1 w hk j 2 V ≤ C.
From (5.24), Lemma 2.1, and (5.25) we obtain (5.13) which completes the proof of the theorem. Similar to Theorem 4.1, Theorem 5.2 can be used to produce convergence order error estimates for the fully discrete approximations with particular choices of the finite dimensional subspace V h . As a sample result, we consider using the linear element spaces {V h } of (4.28).

Corollary 5.3. Keep the assumptions stated in Theorem 5.2. Assume Ω is a polygonal/polyhedral domain, and let {V h } be the family of linear element spaces defined by (4.28), corresponding to a regular family of finite element triangulations of Ω into triangles or tetrahedrons. Let u and {w hk n } N n=0 be solutions of Problems P V and P kh V , respectively. Assume u 0 ∈ H 2 (Ω; R d ), u 1 ∈ H 1 (Ω; R d ),a n dl e tu h 0 , u h 1 ∈ V h be defined by (4.1).L e t{u hk n } N n=0 be defined by (5.5). Under the regularity conditions u ∈ C 1 (0,T;

H 2 (Ω; R d )) ∩ H 3 (0,T; H), uτ ∈ C(0,T; H 2 (Γ 3 ; R d )),
we have the optimal order error estimate (5.26) max

1≤n≤N { u n -u hk n V + w n -w hk n H }≤c(h + k).
Proof.L e tv h j ∈ V h be the finite element interpolant of u j , t ∈ [0,T], 1 ≤ j ≤ N . Note that [START_REF] Han | Quasistatic Contact Problems in Viscoelasticity and Viscoplasticity[END_REF] k

N j=1 ẇj -δw j 2 H ≤ ck 2 u 2 H 2 (0,T ;H) , 1 k N -1 j=1 (w j -v h j ) -(w j+1 -v h j+1 ) 2 H ≤ ch 2 u 2 H 2 (0,T ;V ) .
Then similar to the proof of Corollary 4.3, we obtain (5.26) from (5.13).

6. Numerical simulations. The aim of this section is to present some numerical results to illustrate the behavior of the solution of the frictional contact problem Problem P V . We pay particular attention to the numerical convergence order.

The numerical solution of Problem P V is based on the backward Euler divided difference for the time discretization and the finite element approximation using the linear element space (4.28) for the spatial discretization. To solve the discrete problems, we use a "convexification" iterative procedure [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], which leads to a sequence of convex programming problems. For each "convexification" iteration, the coefficient of friction μ(| uτ |) is fixed to a given value depending on the tangential velocity solution uτ found in the previous iteration. Then, the resulting nonsmooth convex iterative problems are solved. The frictional bilateral condition is treated by using an augmented Lagrangian approach. For details about this numerical method, we refer the reader to [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF][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] Wriggers | Computational Contact Mechanics[END_REF]. For practical implementation of the method, we use additional fictitious nodes for the Lagrange multiplier in the initial mesh. Construction of these nodes depends on the contact elements used for the geometrical discretization of the interface Γ 3 . In our numerical example, the discretization is based on "node-to-rigid" contact element, which is composed of one node of Γ 3 and one Lagrange multiplier node. To keep this paper to a reasonable length, we skip the details of the numerical algorithms and implementation; details on the discretization step and computational contact mechanics, including algorithms similar to that used here, can be found in [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]. Different numerical methods in the study of such frictional problems, including the proximal bundle methods, also can be found in [START_REF] Haslinger | Finite Element Method for Hemivariational Inequalities: Theory, Methods and Applications[END_REF][START_REF] Mistakidis | Numerical treatment of problems involving nonmonotone boundary or stress-strain laws[END_REF][START_REF] Mistakidis | The search for substationary points in the unilateral contact problems with nonmonotone friction[END_REF][START_REF] Tzaferopoulos | Comparison of two methods for the solution of a class of nonconvex energy problems using convex minimization algorithms[END_REF]. Numerical example. We consider the physical setting shown in Figure 1. There, Ω=(0,L

1 ) × (0,L 2 ) ⊂ R 2 with L 1 ,L 2 > 0a n d Γ 1 = {0}×[0,L 2 ], Γ 2 =([0,L 1 ] ×{L 2 }) ∪ ({L 1 }×[0,L 2 ]), Γ 3 =[0,L 1 ] ×{0}.
The domain Ω represents the cross section of a three-dimensional linearly viscoelastic body subjected to the action of tractions in such a way that a plane stress hypothesis is valid. On Γ 1 = {0}×[0,L 2 ] the body is clamped, i.e., the displacement field vanishes there. Vertical compressions act on the part [0, Mechanical behavior of the solution. In Figure 2 we plot the deformed configuration as well as the interface forces on Γ 3 during the dynamic compression process at times t =0 .3 s, t =0 .5 s, t =0 .8 s,a n dt =1 .1 s. At the beginning of the process, the contact nodes are in status of stick, and then at the end of the process, on the right side of Γ 3 , a large proportion of contact nodes switches to status of slip when the compression of the domain is stronger. There, the friction bound has decreased with respect to the evolution of μ(| uτ |) and is reached.

In Figure 3 we plot the deformed meshes and the interface forces on Γ 3 for two different values of the coefficients a and b, respectively. Note that in the case a =1 and b =0.1 considered in Figure 2 the coefficient of friction is a nonmonotone function with respect to the slip rate, while in the cases a = b =0 .1a n da = b =1i ti sa constant. In the case a = b =0.1 we note that all the contact nodes are in slip contact since, there, the friction bound is low and, therefore, is reached. In contrast, in the case a = b = 1 the friction bound is higher and, as a consequence, all the contact nodes are in stick status. Numerical convergence order. In order to check the convergence of the discrete scheme and to illustrate the optimal error estimate obtained in section 5, we computed a sequence of numerical solutions by using uniform discretizations of the problem domain according to the spatial discretization parameter h and time step k. For instance, the deformed configuration and the interface forces plotted in Figure 2 correspond to the choices h =1/128 and k =1/128.

The numerical error u-u hk V is computed for several discretization parameters of h and k. Here, the boundary Γ of Ω is divided into 1/h equal parts. We start with h =1 /2a n dk =1 /2, which are successively halved. The numerical solution corresponding to h =1/256 and k =1/256 was taken as the "exact" solution, used to compute the errors of the numerical solutions; this fine discretization corresponds to a problem with 133, 896 degrees of freedom at each time level. The numerical results are presented in Figure 4, where the dependence of the error estimate uu hk V with respect to h and k is plotted. A first order convergence is clearly observed, providing numerical evidence of the theoretical optimal order error estimate obtained in section 5.

  and we can express (3.20)-(3.21) and (4.3)-(4.4) as follows:

  [START_REF] Ionescu | Functional and Numerical Methods in Viscoplasticity[END_REF], and (5.23),(5.24) e n ≤ Cg n + n j=1 ke j for n =1,...,N with C>0. Note that from (5.6),

Fig. 1 .

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

L 1 ]

 1 ×{L 2 } of the boundary and the part {L 1 }×[0,L 2 ] is traction free. No body forces are assumed to act on the elastic body during the process. The body is in frictional bilateral contact with an obstacle on the part Γ 3 =[ 0 ,L 1 ] ×{0} of the boundary. The friction follows a nonmonotone law which the friction coefficient depends on the tangential velocity | uτ |.F o rt h e coefficient of friction we choose a function μ : R d → R of the form (6.1) μ(| uτ |)=(a -b) e -α | uτ | + bwith a, b, α > 0, a ≥ b. Note that the friction law (3.6) with (6.1) describes the slip weakening phenomenon which appears in the study of geophysical problems; see[START_REF] Scholz | The Mechanics of Earthquakes and Faulting[END_REF] for details. The coefficient of friction decreases with the slip rate from the value a to the limit value b. For this reason, the corresponding friction law is nonmonotone. The compressible material response is governed by a linearly viscoelastic constitutive law in which the viscosity tensor A and the elasticity tensor B are given by(Aτ ) αβ = μ 1 (τ 11 + τ 22 )δ αβ + μ 2 τ αβ , 1 ≤ α, β ≤ 2, ∀ τ ∈ S 2 , (Bτ ) αβ = Eκ (1 + κ)(1 -2κ) (τ 11 + τ 22 )δ αβ + E 1+κ τ αβ , 1 ≤ α, β ≤ 2, ∀ τ ∈ S 2 ,where μ 1 and μ 2 are viscosity constants, E and κ are Young's modulus and Poisson's ratio of the material, and δ αβ denotes the Kronecker symbol.

Fig. 2 .

 2 Fig. 2. Evolution of deformed meshes and frictional contact forces during the dynamic compression process.

For

  computation we use the following data:L 1 =1m, L 2 =0.5 m, ρ = 1000 kg/m 3 ,T =1.1 s, μ 1 =50N/m 2 ,μ 2 = 100 N/m 2 ,E = 2000 N/m 2 ,κ =0.3, f 0 =(0, -10 -5 ) N/m 2 , f 2 = (0, 0) N/m on {L}×[0,L], (0, -600 t) N/m on [0,L 1 ] ×{L 2 }, a =1,b =0.1,α = 200.Our results are presented in Figures2, 3, and 4 and are explained below.

Fig. 3 .

 3 Fig. 3. Deformed meshes and interface forces on Γ 3 corresponding to different values of the coefficients a and b.

Fig. 4 .

 4 Fig. 4. Numerical errors.

  •, • H * Γ ×H Γ denote the duality pairing between H *

Γ and H Γ . For every σ ∈ Q 1 there exists an element σν ∈ H * Γ such that

  Assume that {g n } N n=1 and {e n } N n=1 are two sequences of nonnegative numbers satisfying e n ≤ cg n +c ⊂ R d with a Lipschitz boundary Γ that is partitioned into three parts Γ 1 , Γ 2 ,a n dΓ 3 with Γ 1 ,Γ 2 ,a n dΓ 3 being relatively open and mutually disjoint, and meas (Γ 1 ) > 0. Let [0,T] be a time interval of interest, T>0.

		n
		ke j ,n =1,...,N,
		j=1
	for a positive constant c independent of N or k. Then there exists a positive constant
	c, independent of N or k, such that
	max 1≤n≤N	e n ≤ c max 1≤n≤N	g n .
	3. Mechanical problem and variational formulations. We start with a
	description of the mechanical problem. A linearly viscoelastic body occupies an open
	bounded connected set Ω

16, Chapter 7].

Lemma 2.1. Let T>0 be given. For a positive integer N we define k = T/N.

Acknowledgment. We thank the three anonymous referees for their valuable comments on the original version of the paper.

* This research was supported by the Marie Curie International Research Staff Exchange Scheme Fellowship within the 7th European Community Framework Programme under grant agreement 295118, the National Science Center of Poland under Maestro Advanced Project DEC-2012/06/A/ ST1/00262, project Polonium "Mathematical and Nu-merical Analysis for Contact Problems with Friction" 2014/15 between Jagiellonian University and Universite´ de Perpignan Via Domitia, the National Science Center of Poland under grant N N201 604640, and the International Project cofinanced by the Ministry of Science and Higher Education of Republic of Poland under grant W111/7.PR/2012. The research of this author was also supported by grants from the Simons Foundation.

Let {w hk n } N n=0 be the solution of Problem P kh V and let {u hk n } N n=0 be given by (5.5). Then the following bound holds for all {v h j } N j=1 ⊂ V h :

Proof.T a k i n gt h es a m ev h ∈ V h in (3.20) and (5.2) we obtain for n =1,...,N,

From (5.14) we get

After some reformulation we obtain (5.15)

Using the formula 2(a -b, a

H ≤ ρ(δw n -δw hk n ), w nw hk n H .

By the Lipschitz continuity of A,

(5.17)