Rothe method and numerical analysis for history-dependent hemivariational inequalities with applications to contact mechanics

In this paper an abstract evolutionary hemivariational inequality with a history-dependent operator is studied. First, a result on its unique solvability and solution regularity is proved by applying the Rothe method. Next, we introduce a numerical scheme to solve the inequality and derive error estimates. We apply the results to a quasistatic frictional contact problem in which the material is modeled with a viscoelastic constitutive law, the contact is given in the form of multivalued normal compliance, and friction is described with a subgradient of a locally Lipschitz potential. Finally, for the contact problem we provide the optimal error estimate.


Introduction
In this paper we are concerned with the existence and uniqueness of a solution to an abstract evolutionary hemivariational inequality which involve a history-dependent operator of the form for all v ∈ V , a.e. t ∈ (0, T ) with u(0) = u 0 . Here A and B are operators from a reflexive Banach space V to its dual V * , M is a linear, bounded operator, J 0 denotes the generalized gradient of a locally Lipschitz function, f : (0, T ) → V * and u 0 ∈ V are given, and R represents a history-dependent operator. The motivation to study the inequality of the form (1) comes from contact problems in solid mechanics. It is known that when the external forces and tractions evolve slowly in time in such a way that the acceleration in the system is rather small and negligible, then the inertial terms can be neglected. In such a way, we obtain the quasistatic approximation (equilibrium equation) for the equation of motion. Quasistatic contact models have been studied in several monographs and many papers dedicated to such phenomena, see [9,14,35,36] and the references therein.
In the first part of the paper, we deal with an abstract time-dependent hemivariational inequality of the form (1). The main results are delivered on existence, uniqueness and regularity of a solution to the abstract hemivariational inequality, see Theorem 11. We apply the Rothe method, see [17,18], combined with a surjectivity result for a multivalued and coercive operator. The hemivariational inequality (1) without a history-dependent operator has been recently investigated in [24] by using the vanishing acceleration method, where a local existence result was proved. In contrast to Theorem 17 of [24], here we provide a result on the global unique solvability to (1). Also, our proof is now based on the Rothe method and is simpler, since we have eliminated the additional space Z required in [24]. Moreover, being motivated by applications to contact mechanics in Section 6, the inequality (1) involves a historydependent operator. We recall that the notion of a history-dependent operator is quite recent and it was introduced in [39]. Various problems with history-dependent operators have been studied for the evolution variational and hemivariational inequalities in [1,6,10,11,20,21,22,27,28,29,30,34,43,44], and for the quasistatic problems in [19,26,37,38,40,41,45]. Furthermore, we study a fully discrete approximation for the problem (1) which consists in finite difference discretization in time and finite element approximation in the spatial variable. We prove in Theorem 13 the Céa type error estimate for the hemivariational inequality.
In the second part of the paper, we apply the abstract results to a quasistatic frictional contact model for viscoelastic materials. The process is described by multivalued versions of the nonmonotone normal compliance and friction boundary conditions. We provide the variational formulation of the contact problem for which we deliver a result on its unique weak global solvability. In this way we improve the local existence result of [24,Theorem 17]. Finally, for the frictionless contact we establish a result on an optimal error estimate for the fully discrete approximation scheme. Note that results on numerical anlaysis for hemivariational inequalities can be found in [3,12,15,16,37] and the references therein.
The outline of the paper is as follows. After recalling the basic notation in Section 2, in Section 3 we formulate the abstract hemivariational inequality with a history-dependent operator. In Section 4 we apply the Rothe method to deliver existence and uniquence result for this inequality. The error estimate of the Céa type for a fully discrete approximation is provided in Section 5. Finally, in Section 6, we illustrate the applicability of our results to the quasistatic frictional contact problem for viscoelastic material.

Preliminaries
In this section we recall the basic notation and some results which are needed in the sequel, see [5,7,8,42]. We use the standard notation for the Lebesgue and Sobolev spaces of functions defined on a finite time interval [0, T ] with values in a Banach space. We denote by L(E, F ) the space of linear and bounded operators from a Banach space E to a Banach space F endowed with the usual norm · L(E,F ) . For a subset S of Banach space (E, · E ), we write S E = sup{ s E | s ∈ S}.
Let Y be a reflexive Banach space and ·, · denote the duality of Y and Y * .
converging weakly to y ∈ Y such that lim sup Ay n , y n − y ≤ 0, we have Ay, y − z ≤ lim inf Ay n , y n − z for all z ∈ Y.
Note that the operator A : Y → Y * is pseudomonotone if and only if the conditions y n → y weakly in Y and lim sup Ay n , y n − y ≤ 0 entail lim Ay n , y n − y = 0 and Ay n → Ay weakly in Y * . It is also easy to check that if A ∈ L(Y, Y * ) is nonnegative, then it is pseudomonotone.
We recall the notion of the pseudomonotonicity for a multivalued operator.
for every v ∈ Y , the set T v ⊂ Y * is nonempty, closed and convex, (b) T is upper semicontinuous from each finite dimensional subspace of Y to Y * endowed with the weak topology, (c) for any sequences {u n } ⊂ Y and {u * n } ⊂ Y * such that u n → u weakly in Y , u * n ∈ T u n for all n ≥ 1 and lim sup u * n , u n − u ≤ 0, we have that for every v ∈ Y , there exists u * (v) ∈ T u such that We recall the following fundamental surjectivity theorem, see [8,Theorem 1.3.70] or [42], which will be used to prove existence of a solution to a static hemivariational inequality in Section 4.
Theorem 2. Let Y be a reflexive Banach space and T : Y → 2 Y * be pseudomonotone and coercive. Then T is surjective, i.e., for every f ∈ Y * , there is u ∈ Y such that T u ∋ f .
We hereafter recall the definition of the Clarke subgradient.
Definition 3. Given a locally Lipschitz function J : E → R on a Banach space E, we denote by J 0 (u; v) the generalized (Clarke) directional derivative of J at the point u ∈ E in the direction v ∈ E defined by The generalized gradient of J : E → R at u ∈ E is defined by The following result provides an example of a multivalued pseudomonotone operator which is a superposition of the Clarke subgradient with a compact operator. The proof can be found in [3,Proposition 5.6].
Proposition 4. Let V and X be reflexive Banach spaces, M : V → X be a linear, bounded, and compact operator. We denote by M * : X * → V * the adjoint operator of M. Let J : X → R be a locally Lipschitz function such that We conclude this section with a discrete version of the Gronwall inequality whose proof can be found in [14,Lemma 7.25].

History-dependent hemivariational inequalities
In this section we introduce a class of history-dependent hemivariational inequalities. This class will be studied in Section 4 where the existence and uniqueness result for this class of inequalities will be provided. A fully discrete approximation for the inequalities in this class will be discussed in Section 5. We use the following standard notation, see [7,8,29,42] for details. Let V ⊂ H ⊂ V * be an evolution triple of spaces. Recall that this means that V is a reflexive and separable Banach space, H is a separable Hilbert space, and the embedding V ⊂ H is dense and continuous. Let i be the embedding operator between V and H which is assumed to be compact. It is known that the adjoint operator i * : H → V * is also linear, continuous and compact. The duality pairing between V * and V , and a norm in V , are denoted by ·, · and · , respectively. For the Hilbert space H, we denote its scalar product and a norm by (·, ·) and · H , respectively.
Given 0 < T < +∞, let V = L 2 (0, T ; V ) and H = L 2 (0, T ; H). It follows from the reflexivity of V that both V and its dual space V * = L 2 (0, T ; V * ) are reflexive Banach spaces as well. Identifying H = L 2 (0, T ; H) with its dual, we have the continuous embeddings V ⊂ H ⊂ V * .
The notation ·, · V * ×V stands for the duality pairing between V and V * . Moreover, by C(0, T ; V ) we denote the spave of continuous functions on [0, T ] with values in V .
Let X be a separable and reflexive Banach space. Given operators A, B : V → V * , M : V → X, the function J : X → R, f ∈ V * and u 0 ∈ V , we consider the following evolutionary hemivariational inequality involving a history-dependent operator.
Here R : where E : V → V * , α ∈ V and q : We impose the following assumptions on the data of Problem 6.

H(A):
The operator A : V → V * is linear, bounded, coercive and symmetric, i.e., H(B): The operator B : V → V * is linear, bounded and coercive, i.e., is Lipschitz continuous with respect to the first variable, i.e., there exists L q > 0 such that H(J): The functional J : X → R is such that (i) J is locally Lipschitz.
(ii) There exists c J > 0 such that ∂J(u) X * ≤ c J (1 + u X ) for all u ∈ X.
(iii) There exists m J ≥ 0 such that .
X , for all u, v ∈ X. In addition, examples of nonconvex functions which satisfy the relaxed monotonicity condition can be found in [23,37]. Particularly, it can be proved that for a convex function, condition H(J)(iii) holds with m J = 0.

Rothe method
In this section, we present a result on existence and uniqueness of solution for Problem 6. The technique of proof relies on the Rothe method (known also as a method of lines, see [17,18]). It consists in a time discretization in which we define an approximate sequence of functions by using the implicit (backward) Euler formula. Next, in each time step, we will solve a stationary hemivational inequality. Finally, we construct the piecewise constant and piecewise affine interpolants and prove a convergence result.
In the rest of the section, we denote by C > 0 a constant whose value may change from line to line.
Let N ∈ N be fixed and denote f k for all v ∈ V and for k = 1, 2, . . . , N, where x k τ ∈ V * is defined by First, we shall prove the existence and uniqueness of a solution to Problem 8. Proof. Let u 0 τ , u 1 τ , . . . , u k−1 τ be given. We will prove that there exists a unique element u k τ ∈ V which satisfies inclusion (4). To end this, we apply Theorem 2 to show that the operator L : is bounded, continuous and monotone for τ [23,Theorem 3.69], we conclude that the operator defined by (5) is pseudomonotone. On the other hand, taking into account assumptions H(J)(i)-(ii), H(M) and Proposition 4, it is clear that the operator v → M * ∂J(Mv) is pseudomonotone as well. Therefore, by using [23, Proposition 3.59(ii)], we infer that L is a pseudomonotone operator too.
Subsequently, we prove that the operator L is coercive. From hypothesis H(J) we derive the estimate (see [12]) Hence, we deduce that the operator L is coercive for all τ ∈ (0, τ 0 ). Therefore, by the use of Theorem 2, we obtain that L is surjective, i.e., Problem 8 has at least one solution u k τ ∈ V . For uniqueness part, we assume that u k τ and u k τ are two solutions in V of Problem 8, that is, where the elements x k τ and x k τ are defined by respectively. We take v = u k τ − u k τ in the first inequality and v = u k τ − u k τ in the second one. We add the resulting inequalities to get The smallness condition (H 0 ) guarantees that u k τ = u k τ , which completes the proof of this lemma.
Next, we establish the estimates for the solution of Problem 8.
Proof. We choose v = u k τ in (4), then use the hypotheses H(A), H(B) and the equality Next, the assumptions H(E) and H(q) imply Combining (10) and (11), and using the Cauchy inequality with ε > 0, we have Summing the above inequalities for k = 1, . . . , n, where 1 ≤ n ≤ N, and then applying H(A), we get Now, we use the discrete version of the Gronwall inequality in Lemma 5, to verify estimates (6) and (7). The estimate (8) follows directly from (6) and H(J)(ii). Denote The latter together with (6), (8), H(E), H(q) and the Cauchy inequality with ε > 0 implies So, we obtain the estimate (9), which completes the proof of this lemma.
Subsequently, for a given τ > 0, we define the piecewise affine function u τ and the piecewise constant interpolant functions u τ , ξ τ , f τ and w τ as follows Now, we rewrite Problem 8 in the following equivalent form The main results of this section is delivered in the following theorem.
Proof. The bound (6) ensures that {u τ } is bounded in V due to the following inequality It follows from the reflexivity of V that there exists a function u ∈ V such that, passing to a subsequence again indexed by τ , we have Also, from (6), we have that the sequence {u τ } is bounded in V, and therefore, there exists u 1 ∈ V such that Hence, we get u τ − u τ → u − u 1 weakly in V, as τ → 0. By the Hölder inequality and the boundedness of {u ′ τ } (see (9)) From estimate (15) we deduce that u = u 1 . On the other hand, by the boundedness of {u ′ τ } (see (9)), we also obtain (cf. [ In addition, using the boundedness of {ξ τ } (see (8)) and the reflexivity of the space X * , we conclude By virtue of the hypothesis H(q) and boundedness of {u τ } (see (6)), one has the following estimate for t for some C 0 > 0, which is independent of τ . Moreover, [4,Lemma 3.3] implies that Next, we shall show that u is a solution of Problem 6. To this end, we define the Nemytskii operators A, B : V → V * by (Av)(t) = A(v(t)) and (Bv)(t) = B(v(t)) for all v ∈ V and a.e. t ∈ (0, T ). From hypotheses H(A) and H(B), it is clear that A and B are both linear and bounded, so they are also weakly continuous. Thus from (16) and (13) we obtain Au ′ τ → Au ′ and B u τ → Bu both weakly in V * , as τ → 0, i.e., lim for all v ∈ V. Now, we consider the Nemitskii operators E, for all v ∈ V and a.e. t ∈ (0, T ). It is obvious that E is weakly continuous being bounded and linear. From the convergence (13), one has for all v ∈ V. Next, from H(E), H(q) and (18), we have for all v ∈ V.
Now, we introduce the Nemitskii operator M : V → X defined by (Mv)(t) = M(v(t)) for all v ∈ V and a.e. t ∈ (0, T ), so, from (17), we have for all v ∈ V. From (19)- (21), (23) and (24), we infer that for all v ∈ V with ξ(t) ∈ ∂J(Mu(t)) for a.e. t ∈ (0, T ). Furthermore, we shall show that u ∈ V with u ′ ∈ V is also a solution of Problem 6. Arguing by contradiction, we suppose that u is not a solution to Problem 6. This means there exist a measurable set I ⊂ [0, T ] with meas(I) > 0 and v * ∈ V such that Inserting v = v into (25) and taking account of (26), it follows from [23,Theorem 3.47] that This results a contradiction, so, u ∈ V with u ′ ∈ V is also a solution of Problem 6.
Finally, we will verify that the solution of Problem 6 is unique. Let u 1 and u 2 be two solutions of Problem 6. Then , v + J 0 (Mu 2 (t); Mv) ≥ 0 for all v ∈ V and a.e. t ∈ (0, T ). Taking v = u 2 (t) − u 1 (t) in the first inequality and v = u 1 (t) − u 2 (t) in the second one, we add the resulting inequalities to get We integrate this inequality on [0, t], where t ∈ [0, T ], and use H(A)(ii) and (H 0 ) to deduce Hence

A fully discrete approximation scheme
In this section, we study a fully discrete approximation scheme for the historydependent hemivariational inequality stated in Problem 6. In this method the time variable is discretized by finite difference and the spatial variable is approximated by finite elements. Assume that V h is a finite dimensional subspace of V and u h 0 ∈ V h is an approximation of the initial point u 0 ∈ V . For N ∈ N, N > 0 given, we denote the time step length by k = T N and t n = kn for n = 0, . . . , N. For a continuous function g defined on the interval [0, T ], in the sequel, we will write g n = g(t n ) for n = 0, . . . , N. In addition, for a sequence {u n } N n=0 , we use the notation δu n = u n − u n−1 k , n = 1, . . . , N.
For the history-dependent operator we introduce a modified trapezoidal approximation for R defined by for v = {v j } N j=1 . In addition, if w ∈ C(0, T ; V ), then the expression R k n w is understood as follows Subsequently, we consider the following fully discrete approximation problem for Problem 6.
We will provide an error analysis of the fully discrete approximation (28). Our goal is to prove the Céa type inequality for Problem 12.
First, exploiting the definition of δu hk n , the inequality (28) can be reformulated as follows This inequality represents a stationary hemivariational inequality. When k small enough, from Lemma 9, we know that under the hypotheses H Since A ∈ L(V, V * ) is coercive, in what follows, for a convenience, we introduce the norm · A by v 2 A = Av, v for all v ∈ V , which is equivalent to the norm · V . In the sequel, we denote by C > 0 a constant which may differ from line to line, but it is independent of h and k.
For an error analysis, we have from (1) at t = t n that for all v ∈ V , where R n u = (Ru)(t n ). Denote the errors δ n = δu n − u ′ n and e n = u n − u hk n for n = 1, 2, . . . , N. Taking v = u hk n in (30), one has Aδu n + Bu n + R n u, u hk n − u n + J 0 (Mu n ; Mu hk n − Mu n ) ≥ f n , u hk n − u n + Aδ n , u hk n − u n .
We add (31) and (28) to get Aδu n + Bu n + R n u, u hk n − u n + Aδu hk n + Bu hk Aδ(u n − u hk n ) + B(u n − u hk n ), u hk n − u n + R n u − R k n u hk , u hk n − u n + Aδu hk n + Bu hk n + R k n u hk , v h − u n + J 0 (Mu n ; Mu hk n − Mu n ) +J 0 (Mu hk n ; Mv h − Mu hk n ) ≥ f n , v h − u n + Aδ n , u hk n − u n for all v h ∈ V h . We use the fact that the function v → J 0 (Mu; Mv) is subadditive (see e.g., [23, Proposition 3.23(i)]), to obtain J 0 (Mu hk n ; Mv h − Mu hk n ) ≤ J 0 (Mu hk n ; Mv h − Mu n ) + J 0 (Mu hk n ; Mu n − Mu hk n ).

So, we have
Aδ(u n − u hk n ) + B(u n − u hk n ), u hk n − u n + R n u − R k n u hk , u hk n − u n + Aδu hk n + Bu hk Combining this inequality with the identity and using the hypotheses H(B) and H(J)(iii) (see Remark 7), it follows that for all v h ∈ V h . Furthermore, we introduce a residual type quantity by Using the fact that u ∈ H 1 (0, T ; V ) (see Theorem 11), we have From these inequalities, we obtain where where ξ n ∈ ∂J(Mu n ) and ξ hk n ∈ ∂J(Mu kh n ). Note that the hypothesis H(J)(ii) and u ∈ H 1 (0, T ; V ) imply that the sequence { ξ n X * } is uniformly bounded. It follows from Lemma 10 that { ξ hk n X * } is uniformly bounded as well. Hence, we have Applying (33) again, we obtain Combining (34)- (36) and applying the Cauchy inequality with ε > 0, we have Now, we replace n by l in the above inequality, and then sum it from 1 to n, where 1 ≤ n ≤ N to get This together with the following estimates It follows from the discrete Gronwall inequality, H(A), and Lemma 5, that for all v h n ∈ V h . We now summarize the results of the section in the form of a theorem.
Theorem 13. Suppose that assumptions of Lemma 9 are satisfied. Let u hk ∈ V h and u ∈ H 1 (0, T ; V ) be the solutions of Problems 12 and 6, respectively. Then, we have the estimate The inequality (39) is called the Céa type inequality of the fully discrete approximation problem, Problem 12.

A quasistatic viscoelastic contact problem
In this section we study the quasistatic contact problem between a viscoelastic body and a foundation. The volume forces and surface tractions are supposed to change slowly in time and therefore the acceleration in the system is negligible. Neglecting the inertial terms in the equation of motion leads to the quasistatic approximation for the process. We show that the variational formulation of the quasistatic contact problem is a time-dependent hemivariational inequality in Problem 6. For the latter, we apply the abstract result stated in Theorem 11 and prove a result on existence and uniqueness of weak solution. Further, we use the fully discrete approximation method discussed in Section 5 to study the numerical analysis of this contact problem and establish the result concerning optimal error estimate for the fully discrete scheme.

Mathematical model and its variational formulation
The physical setting of the contact problem is as follows. A deformable viscoelastic body occupies an open bounded subset Ω of R d , d = 2, 3 in applications. The volume forces of density f 0 act in Ω and surface tractions of density f N are applied on Γ 2 . They both can depend on time. We are interested in the quasi-static process of the mechanical state of the body on the time interval [0, T ] with 0 < T < +∞. The boundary Γ = ∂Ω of Ω is assumed to be Lipschitz continuous and it consists of three measurable parts Γ 1 , Γ 2 and Γ 3 which are mutually disjoint, and m(Γ 1 ) > 0. The unit outward normal vector ν exists a.e. on Γ. We suppose that the body is clamped on part Γ 1 , and the body may come in contact with an obstacle over the potential contact surface Γ 3 . We also put Q = Ω × (0, T ), Σ = Γ × (0, T ), Σ 1 = Γ 1 × (0, T ), Σ 2 = Γ 2 × (0, T ) and Σ 3 = Γ 3 × (0, T ). We often do not indicate explicitly the dependence of functions on the spatial variable x ∈ Ω.
Let S d denote the space of d×d symmetric matrices. The canonical inner products and norms on R d and S d are given by In what follows we always adopt the summation convention over repeated indices. Moreover, for a vector ξ ∈ R d , the normal and tangential components of ξ on the boundary are denoted by ξ ν = ξ · ν and ξ τ = ξ − ξ ν ν, respectively. The normal and tangential components of the matrix σ ∈ S d are defined on boundary by σ ν = (σν)·ν and σ τ = σν − σ ν ν, respectively.
We denote by u : Q → R d the displacement vector, by σ : Q → S d the stress tensor and by ε(u) = (ε ij (u)) the linearized (small) strain tensor, where i, j = 1, . . . , d. Recall that the components of the linearized strain tensor are given by ε(u) = 1/2(u i,j + u j,i ), where u i,j = ∂u i /∂x j .
The classical formulation of the contact problem reads as follows.
Problem P. Find a displacement field u : Q → R d and a stress field σ : Q → S d such that The relation (40) represents the equilibrium equation in which "Div" denotes the divergence operator for tensor valued functions defined by Divσ = (σ ij,j ). Equation (41) is the viscoelastic constitutive law with long memory, where A and B are linear viscosity and elasticity operators, and C denotes the relaxation operator. Next, conditions (42) and (43) represent the displacement and the traction boundary conditions. The multivalued relations (44) and (45) are the contact and friction conditions, respectively, in which ∂j ν and ∂j τ denote the Clarke generalized gradients of prescribed locally Lipschitz functions j ν and j τ . Finally, condition (46) represents the initial condition where u 0 denotes the initial displacement. For concrete examples of boundary conditions (44) and (45), we refer to [9,14,23,31,32,33]. Subsequently we introduce the spaces needed for the variational formulation. Let V be a closed subspace of H 1 (Ω; R d ) defined by and H = L 2 (Ω; R d ). Then (V, H, V * ) forms an evolution triple of spaces. Moreover, the trace operator is denoted by γ : V → L 2 (Γ; R d ). Given an element v ∈ V we use the same notation v for the trace of v on the boundary. The space V is equipped with the inner product and the corresponding norm given by and · are equivalent norms on V . In addition, we denote by Q ∞ the space of fourth-order tensor fields given by We assume that the viscosity and elasticity tensors have the usual properties of ellipticity and symmetry.

H(A ) :
A : Ω × S d → S d is a viscosity tensor, A = (a ijkl ) ∈ Q ∞ such that there exists m 1 > 0 satisfying A τ · τ ≥ m 1 τ 2 S d for all τ ∈ S d , a.e. in Ω.
The body forces, surface tractions and initial displacement satisfy The superpotentials satisfy In the hypotheses H(j ν ) and H(j τ ) the subdifferential is taken with respect to the last variables of j ν and j τ , respectively.
Next, we define the operators A, B ∈ L(V, V * ) by for u, v ∈ V , and the operator R : for all w ∈ V, v ∈ V , a.e. t ∈ (0, T ).
To obtain the weak formulation of the problem (40)-(46), we assume the sufficient smoothness of the functions involved, use the equilibrium equation (40) and the Green formula. We obtain for v ∈ V . Taking into account the boundary condition (42) and (43), we have where f ∈ V * is given by On the other hand, by the ortogonality relation, cf. (6.33) in [23], we get The contact and friction boundary conditions (44) and (45) can be equivalently formulated as follows Using (41), (48), (51) and (52), from (50), we obtain the following hemivariational inequality which is a weak formulation of the problem (40)-(46): find u : (53)

Existence and uniqueness for contact problem
Let X = L 2 (Γ 3 ; R d ) and consider the functional J : X → R defined by then the functional J defined by (54) satisfies (i) J is Lipschitz continuous on bounded subsets of X, for all v, w ∈ X, ξ ∈ ∂J(v) and η ∈ ∂J(w), we have with m 3 = m ν + m τ , (iv) for all v, w ∈ X, we have where J 0 (v; w) denotes the directional derivative of J at a point v ∈ X in the direction w ∈ X.
Under our notation we associate with the hemivariational inequality (53), the following inclusion: find u ∈ V such that u ′ ∈ V and    Note that if the hypotheses H(j ν ) and H(j τ ) hold, then every solution to (58) is a solution to (53). The converse holds provided j ν and j τ satisfy the regularity condition (55). These facts follow from the definition of the Clarke generalized gradient and Lemma 14.
The existence, uniqueness and regularity result for the hemivariational inequality (53) is given in the following result.  (49)) and hypothesis H(C ) that H(E) and H(q) are satisfied with E = I and q = C . Moreover, we put M = γ ∈ L(V, X), γ is the trace operator. It is a consequence of Lemma 14 that the functional J given by (54) satisfies H(J) with c J = c 1 and m J = m 3 (see Lemma 14). Also H(M) follows easily by the properties of the trace operator. The conclusion is a consequence of Theorem 11, which completes the proof of this theorem.
We say that a couple of functions (u, σ) which satisfies (41) and (53) is called a weak solution to Problem P. We conclude that, under the assumptions of Theorem 15, Problem P has a unique weak solution. Moreover, the weak solution has the following regularity u ∈ H 1 (0, T ; V ), σ ∈ L 2 (0, T ; L 2 (Ω, S d )) and Div σ ∈ V * .

Numerical analysis of contact problem
In this section, we will apply the results from Section 5 to establish an optimal order error estimate for the fully discrete solution of the contact problem in Problem P. Here, we consider the frictionless boundary condition on Γ 3 , i.e., the frictional boundary (45) will be reduced to σ τ (t) = 0 on Σ 3 .
In addition, without loss of generality, we may assume that u 0 = 0. We use the same the spaces as introduced in Section 6.1. Then, consider the trace operator γ : V → L 2 (Γ 3 ; R d ). It follows from the Sobolev trace theorem that for some constant c 0 > 0, which depends only on Ω, Γ 1 and Γ 3 . Let X = L 2 (Γ 3 ) and define the operators γ ν : We also consider the functional J : X → R defined by If either j ν (x, ·) or −j ν (x, ·) is regular and H(j ν ) holds, then by Lemma 14, (48) and (49), the contact problem (40)-(46) with j τ = 0 has following equivalent variational formulation.
Next, we pass to the numerical approximation of Problem 16. Likewise in Section 5, for an integer N > 0, let k = T N be the time step length. For simplicity, we suppose that Ω is a polygonal/polyhedral domain and express the three parts of the boundary, Γ k , k = 1, 2, 3, as a union of closed flat components with disjoint interiors Subsequently, we consider a regular family of meshes {T h } that partition Ω into triangles/tetrahedrons compatible with the splitting of the boundary ∂Ω into Γ j,i , 1 ≤ i ≤ i j , 1 ≤ j ≤ 3. This means that if the intersection of one side/face of an element with one set Γ j,i has a positive measure with respect to Γ j,i , then the side/face lies entirely in Γ j,i . Corresponding to the family {T h }, we define the linear element space where P 1 (U) d denotes a set of all linear functions whose domain of definition is U (cf. [16, p. 70]). Now, we are in a position to formulate the following fully discrete approximation problem for Problem 16. for all n = 1, 2, . . . , N.
This together with Hölder inequality implies that Hence, we have Recall that Π h u n,ν is the finite element interpolant of u n,ν on each component Γ 3,i . Combining (66) with the hypothesis (63), we get On the other hand, we estimate the residual quantity |S n (v)|. To this end, we use the fact (see (50), (51) and (59)) that to get S n (v) = Γ 3 σ n,ν + ξ n Π h u n,ν − u n,ν dΓ for some ξ n ∈ ∂J(Mu n ). This implies |S n (Π h u n )| ≤ C Π h u n,ν − u n,ν X , and, therefore, we have This estimate together with (65)-(67) implies the following optimal estimate for the fully discrete scheme (62).
Theorem 18. Assume that u and u hk are solutions to Problems 16 and 17, respectively, and the regularity condition (63) holds. Then, we have max 1≤n≤N u n − u hk n V ≤ C(k + h), where C > 0 is independent of k and h.
In the optimal error estimate of Theorem 18, the method is of first-order in spatial mesh size and in the time step.