A nonsmooth optimization approach for hemivariational inequalities with applications in Contact Mechanics

In this paper we introduce an abstract nonsmooth optimization problem and prove existence and uniqueness of its solution. We present a numerical scheme to approximate this solution. The theory is later applied to a sample static contact problem describing an elastic body in frictional contact with a foundation. This contact is governed by a nonmonotone friction law with dependence on normal and tangential components of displacement. Finally, computational simulations are performed to illustrate obtained results.


Introduction
In the literature we can find examples of many models describing displacement of deformable body that is partly in contact with another object, the so-called foundation. In various contact models boundary conditions enforced on the part of the body contacting the foundation appear. Functions that occur in these conditions model response of the foundation in direction normal to the contact boundary and in direction tangential to the boundary (friction law). In many cases these functions are monotone, such as when Coulomb's law of dry friction is considered, but in applications this may not always be the case. What is more, the friction bound may change as the penetration of the foundation by body increases. Nonmonotonicity of functions describing contact laws and influence of normal displacement of the body on friction law cause some difficulties in analytical and numerical treatment of considered problems.
In this paper we introduce an abstract framework that can be used to numerically approximate a solution to a class of mechanical contact problems. We present a nonsmooth optimization problem and prove existence and uniqueness of a solution to this problem. Next we present a numerical scheme approximating this solution and provide numerical error estimation. We apply this theory to a static contact problem describing an elastic body in contact with a foundation. This contact is governed by a nonmonotone friction law with dependence on normal and tangential components of displacement. Weak formulation of introduced contact problem is presented in the form of hemivariational inequality. In the end we show results of computational simulations and describe the numerical algorithm that was used to obtain these results.
Let us now briefly present references in the literature. The definition and properties of Clarke subdifferential and tools used to solve optimization problems were introduced in [6]. Comparison of nonsmooth and nonconvex optimization methods can be found in [1], and details on computational contact mechanics is presented in [16]. The theory of hemivariational inequalities was developed in [15], and the idea to use Finite Element Method to solve these inequalities was presented in [11]. Another early study of vector-valued hemivariational problems in the context of FEM can be found in [12]. More recent analysis of hemivariational and variational-hemivariational inequalities was presented in [13], [14], whereas numerical analysis of such problems can be found for example in papers [2], [3], [4], [8], [9], [10].
A similar mechanical model to the one described in the paper was already considered in [14], where the authors prove only existence of a solution using surjectivity result for pseudomonotone, coercive multifunction without requiring any smallness assumption.
An error estimation concerning stationary variational-hemivariational inequalities was presented in [8]. In our case variational part of inequality is not present and the inequality is not constrained, however error estimations had to be generalized to reflect dependence of friction law on normal component of the displacement.
A numerical treatment of mechanical problem leading to hemivariational inequality using two approaches -nonsmooth and nonconvex optimization and quasi-augmented Lagrangian method is presented in [2]. As the smallness assumption is not required, this once again does not guarantee uniqueness and leads to a nonconvex optimization problem. There, the authors assume contact to be bilateral and consider friction law which does not depend on normal component of the displacement. This paper is organized as follows. Section 2 contains a general differential inclusion problem and an optimization problem. We show that under introduced assumptions both problems are equivalent and have a unique solution. In Section 3 we proceed with a discrete scheme that approximates solution to introduced optimization problem and we prove theorem concerning numerical error estimation. An application of presented theory in the form of mechanical contact model is indicated in Section 4, along with its weak formulation. Finally, in Section 5, we describe computational algorithm used to solve mechanical contact problem and present simulations for a set of sample data.

A general optimization problem
Let us start with basic notation used in this paper. For a normed space X, we denote by · X its norm, by X * its dual space and by ·, · X * ×X the duality pairing of X * and X. By c > 0 we denote a generic constant (value of c may differ in different equations).
Let us now assume that j : X → R is locally Lipschitz continuous. The generalized directional derivative of j at x ∈ X in the direction v ∈ X is defined by The generalized subdifferential of j at x is a subset of the dual space X * given by If j : X n → R is a locally Lipschitz function of n variables, then we denote by ∂ i j and j 0 i the Clarke subdifferential and generalized directional derivative with respect to i-th variable of j, respectively.
Let now V be a reflexive Banach space and X be a Banach space. Let γ ∈ L(V, X) be linear and continuous operator from V to X, and c γ := γ L(V,X) . We denote by γ * : X * → V * the adjoint operator to γ. Let A : V → V * , J : X × X → R and f ∈ V * . We formulate the differential inclusion problem as follows.
In the study of Problem P incl we make the following assumptions.
H(J) : The functional J : X × X → R satisfies (a) J is locally Lipschitz continuous with respect to its second variable, We remark that condition H(J)(c) is a more general form of a relaxed monotonicity condition, i.e. for all w 1 , Moreover, in a special case when J does not depend on the first variable (i.e. w 1 = w 2 = w), condition H(J)(c) is equivalent to a relaxed monotonicity condition, i.e. for all w, v 1 , v 2 ∈ X (2.1) We start with a uniqueness result for Problem P incl .

Lemma 1 Assume that H(A), H(J)
, H(f ) and (H s ) hold. If Problem P incl has a solution u ∈ V , then it is unique and satisfies with a positive constant c.
Proof. Let u ∈ V be a solution to Problem P incl . This means that there exists z ∈ ∂ 2 J(γu, γu) such that Let us now assume that Problem P incl has two different solutions u 1 and u 2 . For a solution Adding the above inequalities, we obtain Hence, H(A)(c) and H(J)(c) yield and finally Under assumption (H s ), we obtain that if Problem P incl has a solution, it is unique.
Using H(J)(b) and (c), we get Combining (2.4) and (2.5), we have From (H s ) we obtain required estimation.
We now consider an optimization problem, which will be equivalent to Problem P incl under introduced assumptions. To this end, let the operator L : The next lemma collects some properties of the operator L.
Proof. The proof of (i) is immediate since for a fixed w ∈ V the operator L(w, ·) is locally Lipschitz continuous as a sum of locally Lipschitz continuous functions with respect to v.
For the proof of (ii), we observe that from H(A) and H(f ), the functions are strictly differentiable and we calculate Now, using the sum and the chain rules for generalized subgradient (c.f. Propositions 3.35 and 3.37 in [14]), we obtain which concludes (ii).
In order to prove (iii), let us fix Hence, using H(A)(c) and (2.1), we obtain From (H s ) we see that ∂ 2 L(w, ·) is strongly monotone for every w ∈ V . This is equivalent to the fact that L(w, ·) is strongly convex for every w ∈ V (see Theorem 3.4 in [7]), which implies that it is strictly convex.
The problem under consideration reads as follows.
We are now in a position to prove the existence and uniqueness result for the above optimization problem. Proof. We introduce operator Λ : V → V defined for all w ∈ V as follows From Lemma 2 (iii) we see that operator Λ is well defined. Now we prove that the operator Λ is a contraction. Let u i = Λu i for u i ∈ V fixed, i = 1, 2. Because of strict convexity of L(w, ·) we have (see Theorem 1.23 in [11]). From similar arguments to those used in proofs of Lemmata 1 and 2 with fixed first argument of operator L, we have for all Because of (H s ), we can rearrange these terms to get Using assumption (H s ) once more, we obtain that the operator Λ is a contraction. From the Banach fixed point theorem we know that there exists a unique u * ∈ V such that Λu * = u * , so 0 ∈ ∂ 2 L(u * , u * ).
Let us conclude the results from Lemmata 1, 2 and 3 in the following theorem.
with a positive constant c.
Proof. Lemma 2 (ii) implies that every solution to Problem P opt solves Problem P incl . Using this fact, Lemmata 1 and 3 we see that a unique solution to Problem P opt is also a unique solution to Problem P incl . Because of the uniqueness of the solution to Problem P incl we get that Problems P incl and P opt are equivalent. The estimation in the statement of the theorem follows from Lemma 1.

Numerical scheme
Let V h ⊂ V be a family of finite dimensional subspaces with a discretization parameter h > 0. We present the following discrete scheme of Problem P opt .
We remark that existence of a unique solution to Problem P h opt and equivalence to the discrete version of Problem P incl follow from application of Theorem 4 in this new setting. Now let us present the following main theorem concerning error estimation of introduced numerical scheme.
where a residual quantity is given by Proof. Let u be a solution to Problem P opt and u h be a solution to Problem P h opt . Then they are solutions to corresponding inclusion problems and satisfy respectively We observe that by subadditivity of generalized directional derivative (cf. [14], Proposition 3.23(i)) and H(J)(c), we have From the statement of Lemma 1 applied to discrete version of Problem P incl we get that γu h X ≤ c γ u h V ≤ c (1 + f V * ) is uniformly bounded with respect to h. Hence, returning to (3.5) and using (3.6), we obtain for all

By assumption H(A) and definition (3.2), we get for all
Finally, the elementary inequality ab ≤ εa 2 + b 2 4ε with ε > 0 yields

This is equivalent for all
Taking sufficiently small ε and using (H s ) we obtain the desired conclusion.

Application to Contact Mechanics
In this section we apply the results of previous sections to a sample mechanical contact problem. Let us start by introducing the physical setting and notation useful in the problem. An elastic body occupies a domain Ω ⊂ R d , where d = 2, 3 in application. We assume that its boundary Γ is divided into three disjoint measurable parts: Γ D , Γ C , Γ N , where the part Γ D has a positive measure. Additionally Γ is Lipschitz continuous, and therefore the outside normal vector ν to Γ exists a.e. on the boundary. The body is clamped on Γ D , i.e. its displacement is equal to 0 on this part of boundary. A surface force of density f N acts on the boundary Γ N and a body force of density f 0 acts in Ω. The contact phenomenon on Γ C is modeled using general subdifferential inclusions. We are interested in finding the displacement of the body in a static state.
Let us denote by "·" and · the scalar product and the Euclidean norm in R d or S d , respectively, where S d = R d×d sym . Indices i and j run from 1 to d and the index after a comma represents the partial derivative with respect to the corresponding component of the independent variable. Summation over repeated indices is implied. We denote the divergence operator by Div σ = (σ ij,j ). The standard Lebesgue and Sobolev spaces L 2 (Ω) d = L 2 (Ω; R d ) and H 1 (Ω) d = H 1 (Ω; R d ) are used. The linearized (small) strain tensor for displacement u ∈ H 1 (Ω) d is defined by ε(u) = (ε ij (u)), ε ij (u) = 1 2 (u i,j + u j,i ).
Let u ν = u · ν and σ ν = σν · ν be the normal components of u and σ, respectively, and let u τ = u − u ν ν and σ τ = σν − σ ν ν be their tangential components, respectively. In what follows, for simplicity, we sometimes do not indicate explicitly the dependence of various functions on the spatial variable x. Now let us introduce the classical formulation of considered mechanical contact problem.
Problem P : Find a displacement field u : Ω → R d and a stress field σ : Ω → S d such that in Ω (4.1) Div Here, equation (4.1) represents an elastic constitutive law and A is an elasticity operator. Equilibrium equation (4.2) reflects the fact that problem is static. Equation (4.3) represents clamped boundary condition on Γ D and (4.4) represents the action of the traction on Γ N . Inclusion (4.5) describes the response of the foundation in normal direction, whereas the friction is modeled by inclusion (4.6), where j ν and j τ are given superpotentials, and h τ is a given friction bound. We consider the following Hilbert spaces endowed with the inner scalar products respectively. The fact that space V equipped with the norm · V is complete follows from Korn's inequality, and its application is allowed because we assume that meas(Γ D ) > 0. We consider the trace operator γ : V → L 2 (Γ C ) d = X. By the Sobolev trace theorem we know that γ ∈ L(V, X) with the norm equal to c γ . Now we present the hypotheses on data of Problem P .
x ∈ Γ C , (c) there exist c ν0 , c ν1 ≥ 0 such that We remark that condition H(j τ )(b) is equivalent to the fact that j τ (x, ·) is locally Lipschitz continuous and there exists c τ ≥ 0 such that ∂ 2 j τ (x, ξ) ≤ c τ for all ξ ∈ R d and a.e. x ∈ Γ C .
Using the standard procedure, the Green formula and the definition of generalized subdifferential, we obtain a weak formulation of Problem P in the form of hemivariational inequality.
Here, the operator A : V → V * and f ∈ V * are defined for all w, v ∈ V as follows It is easy to check that under assumptions H(A) and (H 0 ), the operator A and the functional f satisfy H(A) and H(f ), respectively. We also define the functional J : Below we present some properties of the functional J. The proof of the following lemma is similar to the proof of Corollary 4.15 in [14] and is skipped here.
With the above properties, we have the following existence and uniqueness result for Problem P hvi .
with a positive constant c.
Proof. We notice that the assumptions of Theorem 4 are satisfied. This implies that Problem P incl has a unique solution. By (2.3) and Corollary 4.15 (iii) in [14] we get that every solution to Problem P incl solves Problem P hvi . Using similar technique as in the proof of Lemma 1, we can show that if Problem P hvi has a solution, it is unique. Combining these facts we obtain our assertion.
We conclude this section by providing a sample error estimate under additional assumptions on the solution regularity. We consider a polygonal domain Ω and a space of continuous piecewise affine functions V h . We introduce the following discretized version of Problem P hvi .
Proof. We denote by Π h u ∈ V h the finite element interpolant of u. By the standard finite element interpolation error bounds (see [5]) we have for all η ∈ H 2 (Ω) d We now bound the residual term defined by (3.2) using similar procedure to one described in [8]. Let v = ±w in inequality (4.7), where the arbitrary function w ∈ V is such that w ∈ C ∞ (Ω) d and w = 0 on Γ D ∪ Γ C . Then we obtain the identity From this identity, using fundamental lemma of calculus of variations, we can deduce that Div A(ε(u)) + f 0 = 0 in Ω, (4.13) (4.14) We multiply equation (4.13) by v h − u and obtain Using the homogenous Dirichlet boundary condition of v h −u on Γ D and the traction boundary condition given by (4.14) we have Using this and (3.2) we obtain From inequalities (3.1), (4.11), (4.12) and (4.17) we get and we obtain required estimation.

Simulations
In this section we present results of our computational simulations. From Theorems 4 and 7 we know that Problems P hvi and P opt are equivalent. Hence, we can apply numerical scheme P h opt and use Theorem 5 to approximate solution of P hvi . We employ Finite Element Method and use space V h of continuous piecewise affine functions as a family of approximating subspaces. The idea for algorithm used to calculate solution of discretized problem is based on the proof of Lemma 3 and is described by Algorithm 1.

Algorithm 1 Iterative optimization algorithm
Let ε > 0 and u h 0 be given In order to minimize not necessarily differentiable function L(w h , ·) we use Powell's conjugate direction method. For a starting point u h 0 we take a solution to problem with σν = 0 on Γ C , although it can be chosen arbitrarily. The elasticity operator A is defined by A(τ ) = 2ητ + λtr(τ )I, τ ∈ S 2 .
Both functions j ν and j τ are nondifferentiable and nonconvex. In Figure 1 we present output obtained for chosen data. We push the body down and to the left with a force f 0 . As a result the body penetrates the foundation, but because of frictional forces it is squeezed to the left more in the higher part than in the lower part. In order to illustrate the error estimate obtained in Section 4, we present a comparison of numerical errors u − u h V computed for a sequence of solutions to discretized problems. We use a uniform discretization of the problem domain according to the spatial discretization parameter h. The boundary Γ C of Ω is divided into 1/h equal parts. We start with h = 1, which is successively halved. The numerical solution corresponding to h = 1/128 was taken as the exact solution u. The numerical results are presented in Figure 2, where the dependence of the error estimate u − u h V with respect to h is plotted on a log-log scale. A first order convergence can be observed, providing numerical evidence of the theoretical optimal order error estimate obtained at the end of Section 4.