Solving Quasistatic Contact Problems Using Nonsmooth Optimization Approach
Michal Jureczka Anna Ochal Weimin Han
math.NA cs.NA 35Q74, 49J40, 65K10, 65M60, 74S05, 74M15, 74M10, 74G15







Solving quasistatic contact problems
 using nonsmooth optimization approach





Michal Jureczka[Jagiellonian University in Krakow, Faculty of Mathematics and Computer Science,
Lojasiewicza 6, 30-348 Krakow, Poland. Email: michal.jureczka@uj.edu.pl],  Anna Ochal[Jagiellonian University in Krakow, Faculty of Mathematics and Computer Science,
Lojasiewicza 6, 30-348 Krakow, Poland. Email: ochal@ii.uj.edu.pl]  and   Weimin Han[Program in Applied Mathematical and Computational Sciences (AMCS) &
Department of Mathematics, University of Iowa, Iowa City, IA 52242, USA.
Email: weimin-han@uiowa.edu]






Abstract.
This paper is devoted to a study of time-dependent hemivariational inequality. We prove existence and uniqueness of its solution, provide fully discrete scheme and reformulate this scheme as a series of nonsmooth optimization problems. This theory is later applied to a sample quasistatic contact problem describing a viscoelastic body in frictional contact with a foundation. This contact is governed by a nonmonotone friction law with dependence on normal component of displacement and tangential component of velocity. Finally, computational simulations are performed to illustrate obtained results.

Keywords. Nonmonotone friction, optimization problem, error estimate, finite element method, numerical simulations.

AMS Classification. 35Q74, 49J40, 65K10, 65M60, 74S05, 74M15, 74M10, 74G15







§ INTRODUCTION


Currently various engineering applications require scientific analysis of contact phenomena. In order to study such problems a framework in the form of variational and hemivariational inequalities complemented by theory of Finite Element Method emerged. Literature in the field of Contact Mechanics grows rapidly. Main ideas and mathematical tools were introduced in monographs <cit.>. Other comprehensive studies include <cit.>. Theory of nonsmooth optimization used in our approach was presented in <cit.> and applied to a sample static contact problem in <cit.>. Numerical analysis of contact problems can be found for instance in papers <cit.> and for a recent study <cit.>.

This paper is a continuation of our previous work <cit.>, where we reformulated time-independent hemivariational inequality as optimization problem and used it to solve a static mechanical contact problem. Here, we consider time-dependent hemivariational inequality, and use similar idea to solve it numerically. In this case we prove existence and uniqueness of a solution using fixed point argument and apply this abstract framework to solve a quasistatic mechanical contact problem. This approach allows us to substitute Coulomb's law of dry friction by more general, nonmonotone friction law with dependence on tangential component of velocity. Additionally, the friction bound is determined by normal component of displacement, describing penetration of the foundation.

The rest of this paper is organized as follows. In Section <ref> we consider an abstract problem in the form of hemivariational inequality and show that under appropriate assumptions it has a unique solution. Section <ref> contains a discrete scheme which approximates solution to introduced abstract problem. We also prove the theorem concerning numerical error estimate. An application of presented theory to a mechanical contact model is introduced in Section <ref>. Finally, in Section <ref>, we describe computational algorithm used to solve mechanical contact problem. Then we present results of simulations and empirical error estimation for a set of sample data.








§ HEMIVARIATIONAL INEQUALITY


 We 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 →ℝ is locally Lipschitz continuous. The generalized directional derivative of j at x ∈ X in the direction v∈ X is defined by

    j^0(x;v) := lim sup_y → x, λ↘ 0j(y+λ v)  - j(y)/λ.

The generalized subdifferential of j at x is a subset of the dual space X^* given by

    ∂ j(x) := {ξ∈ X^*  |  ⟨ξ, v⟩_X^*× X≤ j^0(x;v)   v ∈ X }.

If j X^n →ℝ is a locally Lipschitz function of n variables, then we denote by ∂_i j and j_i^0 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 γ∈ℒ(V, X) be a linear and continuous operator from V to X, and c_γ:=γ_ℒ(V, X). We denote by  γ^*  X^* → V^* the adjoint operator to γ. Let [0,T] be a time interval with T>0.   Let A, B  V → V^*, J  X × X →ℝ, f  [0,T] → V^* and let u_0 ∈ V. The functional J is assumed to be locally Lipschitz continuous with respect to its second argument. We formulate the considered hemivariational inequality as follows.

Problem P:  Find u  [0,T] → V such that for all t ∈ [0,T],

    ⟨  Au'(t) + B u(t), w⟩_V^*× V + J_2^0(γ u(t),γ u'(t);γ w )
       ≥⟨  f(t),  w⟩_V^*× V  ∀ w∈ V,

and

    u(0) = u_0.



We introduce the following assumptions.

H(A):   Operator A  V → V^* satisfies


  (a)
     A is linear and bounded,
  
  (b)
    A is symmetric, i.e., ⟨ A u, v⟩_V^*× V  = ⟨ A v ,  u  ⟩_V^*× V for all u, v ∈ V,
  
  (c)
     there exists m_A>0 such that ⟨ A u, u⟩_V^*× V≥ m_A u_V^2 for all u ∈ V.




H(B):   Operator B  V → V^* is Lipschitz continuous, i.e., there exists L_B > 0 such that

 a     Bv - Bw_V^*≤ L_B v-w_V  v, w ∈ V. 


H(J):    J  X × X →ℝ satisfies

  
  (a)
     J(v,·) is locally Lipschitz continuous for all v∈ V,
  
  (b)
     there exist c_0, c_1, c_2 ≥ 0 such that ∂_2 J(w, v)_X^*≤ c_0+c_1w_X+c_2v_X

     for all w,  v ∈ X,
  
  (c)
    there exist m_J1, m_J2≥ 0 such that

    J_2^0(w_1, v_1; v_2 - v_1) + J_2^0(w _2, v_2; v_1 - v_2)
    ≤ m_J1 v_1- v_2_X^2+ m_J2 w_1 -  w_2 _X v_1 -  v_2_X 

    for all w_1,  w_2,  v_1,  v_2 ∈ X.


H(f):    f ∈ C([0,T]; V^*) .

H(u_0):    u_0 ∈ V.

(H_s):    m_A - 2m_J1c_γ^2 >0.

We remark that assumption H(J)(c) is equivalent to the following condition

    ⟨∂_2 J(w_1,  v_1) - ∂_2 J(w_2,  v_2), v_1 -  v_2 ⟩_X^*× X≥ - m_J1v_1-v_2_X^2 - m_J2 w_1- w_2_X  v_1- v_2_X

for all w_1,w_2,v_1,v_2∈ X. In the special case when J does not depend on its first variable, we obtain a relaxed monotonicity condition, i.e. for all  v_1, v_2 ∈ X

    ⟨∂ J(v_1) - ∂ J(v_2), v_1 -  v_2 ⟩_X^*× X≥ - m_J1 v_1- v_2_X^2.

So, condition H(J)(c) is a generalization of (<ref>).

Define an operator K  L^2(0,T;V)→ C([0,T];V) by the formula

    Kv(t)= ∫_0^t v(s)   ds + u_0.

We can reformulate Problem P as the following one.

Problem P_hvi:  Find v  [0,T] → V such that for all t ∈ [0,T],

    ⟨  Av(t) + B (K v)(t),w⟩_V^*× V + J_2^0(γ (Kv)(t), γ v(t);γ w)
     ≥⟨  f(t),  w⟩_V^*× V  ∀ w∈ V.



Let us now present the following theorem.


Under assumptions H(A), H(B), H(J), H(f), H(u_0) and (H_s), Problem P_hvi has a unique solution v∈ C([0,T];V).


We use a fixed point argument. Given η∈ C([0,T];V), define

    y_η = Kη.

Then y_η∈ C([0,T];V). Consider the auxiliary problem of finding a function v_η [0,T] → V such that

    ⟨  Av_η(t),w⟩_V^*× V + J_2^0(γ y_η(t),γ v_η(t); γ w)
    ≥⟨ f(t)-B y_η(t),w⟩_V^*× V

for all w∈ V. Applying Theorem 4.2 in <cit.> with φ≡ 0, we know that there exists a unique element v_η(t)∈ V
which solves this inequality for each t∈[0,T]. Let us show that the function v_η∈ C([0,T];V).
For simplicity, with t_1, t_2 ∈ [0,T], denote 
v_η(t_i) = v_i, y_η(t_i) = y_i,
f(t_i) = f_i for i = 1,2. We get that

    ⟨ Av_1,w⟩_V^*× V+J_2^0(γ y_1,γ v_1;γ w)   ≥⟨ f_1-By_1,w⟩_V^*× V,
    
    ⟨ Av_2,w⟩_V^*× V+J_2^0(γ y_2,γ v_2;γ w)   ≥⟨ f_2-By_2,w⟩_V^*× V.

Taking w = v_2-v_1 in (<ref>), w = v_1-v_2 in (<ref>) and adding, we get

    ⟨ Av_1-Av_2,v_1-v_2⟩_V^*× V   ≤ J_2^0(γ y_1,γ v_1;γ v_2-γ v_1) + J_2^0(γ y_2,γ v_2;γ v_1-γ v_2)
         + ⟨  f_1 - B y_1 - f_2 + B y_2, v_1 - v_2⟩_V^*× V.

Using strong monotonicity of the operator A guaranteed by H(A)(c), assumptions  H(B), H(J)(c) and H(f) we obtain

    m_Av_1-v_2_V^2    ≤ m_J1γ v_1-γ v_2_X^2+m_J2γ y_1-γ y_2_Xγ v_1-γ v_2_X
         + (f_1-f_2_V^* + L_By_1 - y_2_V) v_1 - v_2_V
       ≤ m_J1c_γ^2 v_1-v_2_V^2+(f_1-f_2_V^*+(m_J2c_γ^2+L_B) y_1-y_2_V)v_1-v_2_V,

i.e.,

    (m_A-m_J1c_γ^2) v_1-v_2_V≤f_1-f_2_V^* +(m_J2c_γ^2+ L_B)y_1-y_2_V.

By the smallness assumption (H_s), m_A-m_J1c_γ^2>0, and from continuity of y_η and f,
we deduce that [0,T] ∋ t ↦ v_η (t) ∈ V
is a continuous function. This allows us to define an operator Λ C([0,T];V)→ C([0,T];V) via the relation

    Λη = v_η  η∈ C([0,T];V).


Let us prove that the operator Λ has a unique fixed point η^* ∈ C([0,T];V). For two arbitrary functions
η_1,η_2∈ C([0,T];V), let v_i be the solution of (<ref>) for η = η_i, i=1,2.
Similarly to (<ref>), we have

    v_1(t) - v_2(t)_V≤ c y_1(t) - y_2(t)_V.

Since

    y_1(t) - y_2(t)= ∫_0^t (η_1(s) - η_2(s))   ds,

we derive from the previous inequality that

    Λη_1(t) - Λη_2(t)_V ≤ c∫_0^t η_1(s) - η_2(s)_V  ds.

This shows that Λ is a history dependent operator. Applying Theorem 3.20 in <cit.>, we know that the operator Λ has a unique fixed point η^* ∈ C([0,T]; V). Moreover, by the definition of Λ, η^* is a solution to Problem P_hvi. Uniqueness of a solution to Problem P_hvi is a consequence of uniqueness of fixed point.





§ NUMERICAL SCHEME


Let us now fix h > 0 and let V^h ⊂ V be a finite dimensional space with a discretization parameter h>0.
For a given N∈ℕ, we introduce the time step k = T/N and the temporal nodes t_j = jk, 0≤ j≤ N.
Let u_0^h ∈ V^h be an approximation of element u_0 in space V^h. We write
v^hk ={v_j^hk}_j=1^N ⊂ V^h and define

    (K^k v^hk)_0 = u_0^h    
      (K^k v^hk)_j = k ∑_i=1^j v_i^hk + u_0^h,    1≤ j ≤ N.

We present the following discretized version of Problem P_hvi in the form of an operator inclusion problem.
Problem P_hvi^h:  Find v^hk ={v_j^hk}_j=1^N ⊂ V^h such that for every w^h∈ V^h and j ∈{1,…,N}

    ⟨  A v^hk_j + B (K^k v^hk)_j-1,  w^h - v^hk_j  ⟩_V^*× V + J_2^0(γ (K^k v^hk)_j-1, γ v^hk_j;  γ w^h - γ v^hk_j )
         ≥⟨  f_j,  w^h - v^hk_j  ⟩_V^*× V.



Problem P_incl^h: Find v^hk = {v^hk_j}_j=1^N ⊂ V^h such that for all j ∈{1,…,N}

    A v^hk_j + B (K^kv^hk)_j-1 + γ^* ∂_2 J (γ  (K^k v^hk)_j-1, γ v_j^hk) ∋ f_j.


Now we introduce some preliminary material, namely we recall a special case of the Jensen inequality and the discrete version of Gronwall inequality. On several occasions, we will also apply the elementary inequality

    (∑_i=1^l a_i)^2 ≤ l∑_i=1^l a_i^2,   a_1,⋯,a_l ∈ℝ.


  (the Jensen inequality)
Let I ⊂ℝ be a set of positive measure and let f I →ℝ be an integrable function. Then

    (1/|I|∫_I f(s)  ds )^2 ≤1/|I|∫_I (f(s))^2  ds.

In particular, for t_m, t_m+1∈ I, t_m < t_m+1, t_m+1 - t_m = k, we have

    (∫_t_m^t_m+1 f(s)  ds )^2 ≤  k ∫_t_m^t_m+1 (f(s))^2  ds.



  (the Gronwall inequality, <cit.>) Let T be given. For N>0 we define k=T/N. Let {g_n}_n=1^N, {e_n}_n=1^N be two nonnegative sequences satisfying for c>0 and for all n∈{1,…,N}

    e_n≤ cg_n + ck∑_j=1^n-1 e_j.

Then there exists a constant ĉ>0 such that

    max_1≤ n≤ N e_n ≤ĉmax_1≤ n≤ N g_n.



Let us now prove the following lemma.

 
Under assumptions H(A), H(B), H(J), H(f), H(u_0) and (H_s) if Problem P_incl^h has a solution v^hk, then it is unique and satisfies

    v_j^hk_V ≤ c(1+ f_j _V^*)

for all j ∈{1,…,N} with a positive constant c.


To simplify the notation in this proof, we write v instead of v^hk.
Let v be a solution to Problem  P_incl^h and let us fix j ∈{1,…,N}. This means that there exists ζ_j ∈∂_2 J (γ  (K^k v)_j-1, γ v_j ) such that

    A v_j + B (K^kv)_j-1 + γ^* ζ_j = f_j.

From the definition of generalized subdifferential of J(γ (K^k v)_j-1, ·) we have for all w^h ∈ V^h,

    ⟨ f_j - A v_j - B (K^kv)_j-1,  w^h⟩_V^*× V   = ⟨γ^* ζ_j,  w^h⟩_V^*× V
        = ⟨ζ_j , γ w^h⟩_X^*× X
       ≤ J_2^0(γ (K^k v)_j-1, γ v_j;  γ w^h).

After reformulation, we obtain discretized version of hemivariational inequality P_hvi,

    ⟨ A v_j+B(K^kv)_j-1,w^h⟩_V^*× V+J_2^0(γ (K^k v)_j-1,γ v_j;γ w^h)
    ≥⟨ f_j, w^h ⟩_V^*× V.

Let us now assume that Problem P_incl^h has two solutions v^1 and v^2. We will prove inductively
that these solutions are equal.
For j = 1 we get (K^kv^1)_0 = (K^kv^2)_0 = u_0^h. In inequality (<ref>) for a solution v^1 we set
w^h =v^2_1-v^1_1, and for a solution v^2 we set w^h  =  v_1^1-v^2_1.
Then adding these inequalities, we obtain

    ⟨ A v_1^1 - A v_1^2,  v_1^1 - v_1^2 ⟩_V^*× V≤ J_2^0(γ  u_0^h, γ v_1^1; γ v_1^2 - γ v_1^1 ) +J_2^0(γ  u_0^h, γ v_1^2; γ v_1^1 - γ v_1^2 ).

Hence, H(A)(c) and H(J)(c) yield

    (m_A - m_J1c_γ^2)   v_1^1 - v_1^2_V^2 ≤ 0,

and consequently from (H_s) we have v^1_1 = v^2_1. We now show that if v_i^1 = v_i^2 for i=1,…,j-1, then v_j^1 = v_j^2. Similarly, in (<ref>) for a solution v^1 we set w^h  =  v^2_j-v^1_j, and for a solution v^2 we set  w^h  =  v_j^1-v^2_j. Adding these inequalities we obtain

    ⟨ A v_j^1 - A v_j^2,  v_j^1 - v_j^2 ⟩_V^*× V   ≤⟨ B (K^kv^1)_j-1 - B (K^kv^2)_j-1,  v_j^2 - v_j^1 ⟩_V^*× V
         +J_2^0(γ (K^k v^1)_j-1, γ v_j^1; γ v_j^2 - γ v_j^1 )
    +J_2^0(γ (K^k v^2)_j-1, γ v_j^2; γ v_j^1 - γ v_j^2 ).

We observe that (K^kv^1)_j-1 = (K^kv^2)_j-1, hence

    ⟨ A v_j^1 - A v_j^2,  v_j^1 - v_j^2 ⟩_V^*× V≤ J_2^0(γ (K^k v^1)_j-1, γ v_j^1; γ v_j^2 - γ v_j^1 ) +J_2^0(γ (K^k v^2)_j-1, γ v_j^2; γ v_j^1 - γ v_j^2 ).

Again H(A)(c) and H(J)(c) yield

    (m_A - m_J1c_γ^2)   v_j^1 - v_j^2_V^2 ≤ 0.

Under assumption (H_s), we obtain that v^1_j = v^2_j. This equality holds for j=1,…,N, hence a solution of P_incl^h is unique.

Now, in order to prove (<ref>), we set w^h  =-v_j in (<ref>) to obtain

    ⟨ A v_j+B(K^kv)_j-1,-v_j⟩_V^*× V+J_2^0(γ (K^k v)_j-1,γ v_j;-γ v_j
    ≥⟨ f_j,  - v_j ⟩_V^*× V.

Hence,

    ⟨ A v_j,  v_j ⟩_V^*× V≤ J_2^0(γ (K^k v)_j-1, γ v_j; -γ v_j)
    +⟨ B (K^kv)_j-1, - v_j ⟩_V^*× V  + ⟨  f_j , v_j ⟩_V^*× V.





Using H(J)(c), we get

    J_2^0(γ (K^k v)_j-1, γ v_j;  - γ v_j ) + J_2^0(0, 0; γ v_j )
         ≤ m_J1γ v_j_X^2 + m_J2γ (K^k v)_j-1_Xγ v_j_X
       ≤ m_J1c_γ^2 v_j_V^2 + m_J2c_γ^2 (K^k v)_j-1_Vv_j_V.

From Proposition 3.23(ii) in <cit.> and assumption H(J)(b) we have

    J_2^0(0, 0; γ v_j ) ≤ c_0γ v_j_X ≤ c_0 c_γ v_j_V.

Hence,

    J_2^0(γ (K^k v)_j-1, γ v_j;  - γ v_j ) ≤ m_J1c_γ^2 v_j_V^2 + m_J2c_γ^2 (K^k v)_j-1_V  v_j_V + c_0c_γ v_j_V.

From H(B),

    ⟨ B (K^kv)_j-1, v_j ⟩_V^*× V≤B(K^k v)_j-1_V^* v_j_V ≤ (L_B  (K^k v)_j-1_V +c)  v_j_V.

Using H(A)(c), (<ref>) and (<ref>) in (<ref>), we get

    m_A v_j_V^2 ≤ m_J1c_γ^2 v_j_V^2 + (L_B + m_J2c_γ^2)(K^k v)_j-1_Vv_j_V
    + cv_j_V + f_j_V^*v_j_V

Then

    (m_A-m_J1c_γ^2-ε)v_j_V^2≤ c_ε(K^k v)_j-1_V^2+(f_j_V^*+c) v_j_V
    
    .

Taking sufficiently small ε > 0, using (H_s) and inequality (<ref>),
we obtain

    v_j_V^2 ≤ ck^2N∑_i=1^j-1 v_i_V^2 + c u_0_V^2 + c (f_j_V^* + 1)  v_j_V

From discrete version of the Gronwall inequality (Lemma <ref>) with

    e_j = v_j_V^2,     
    g_j = c u_0_V^2 + c(f_j_V^* + 1)  v_j_V
    
    ,

we get

    v_j_V^2 ≤  c(f_j_V^* + 1)  v_j_V,

and then

    v_j_V ≤  c(1+f_j_V^*).

which concludes the proof of (<ref>).


We now consider an optimization problem, which is equivalent to Problem P_incl^h under the stated assumptions.
To this end, let us fix j ∈{1,…,N} and assume v_0,…,v_j-1 to be known, hence also (K^kv)_j-1
is given. Let a functional ℒ_j V →ℝ be defined for all w ∈ V as follows

    ℒ_j(w)=1/2⟨ Aw ,w⟩_V^*× V
    +⟨ B (K^kv)_j-1 -  f_j,  w ⟩_V^*× V + J (γ (K^k v)_j-1, γ w).


The next lemma collects some properties of the functional ℒ_j.

 
Under assumptions  H(A), H(B), H(J), H(f), H(u_0) and (H_s), the functional ℒ_j V →ℝ defined by (<ref>) has the following properties

 
  * ℒ_j is locally Lipschitz continuous,
 
  * ∂ℒ_j( w ) ⊆ A w + B (K^kv)_j-1 - f_j + γ^*  ∂_2 J(γ (K^k v)_j-1, γ w) for all w ∈ V,
 
  * ℒ_j is strictly convex,
 
  * ℒ_j is coercive.



The proof of (i) is immediate since for v ∈ V the functional ℒ_j is locally Lipschitz continuous as a sum of locally Lipschitz continuous functions with respect to w.

For the proof of (ii), we observe that from H(A), H(B) and H(f), the functions

    ψ_1 V ∋  w ↦1/2⟨ A w , w  ⟩_V^*× V∈ℝ,   ψ_2 V ∋ w ↦⟨  B (K^kv)_j-1 - f_j, w  ⟩_V^*× V∈ℝ

are strictly differentiable with

    ψ_1'( w ) = Aw,  ψ_2'( w ) = B (K^kv)_j-1 - f_j.

Now, using the sum and the chain rules for generalized subgradient (cf. Propositions 3.35 and 3.37 in <cit.>), we obtain


    ∂ℒ_j(w )    = ψ_1'( w )+ψ_2'( w ) + ∂_2 (J∘γ)(γ (K^k v)_j-1, w)
    
       ⊆ A w + B (K^kv)_j-1 - f_j + γ^* ∂_2 J(γ (K^k v)_j-1, γ w),

which concludes (ii).

In order to prove (iii), let us fix w^1, w^2 ∈ V. We take η_i ∈∂ℒ_j( w^i) for i=1, 2.
From (ii), there exists ζ_i ∈∂_2J(γ (K^k v)_j-1 , γ w^i) such that

    η_i = A w^i + B (K^kv)_j-1 - f_j + γ^* ζ_i.

From the equivalent condition to H(J)(c), and consequently from (<ref>), we have

    ⟨∂_2J(γ (K^k v)_j-1 , γ w^1)
     - ∂_2J(γ (K^k v)_j-1 , γ w^2), γ w^1 - γ w^2  ⟩_X^*× X
    
      ≥ -m_J1γ w^1 - γ w^2  _X^2.

Hence and from H(A)(c), we obtain

    ⟨η  _1 -  η  _2,  w^1 -  w^2 ⟩_V^*× V   =⟨ A w^1 - A w^2,  w^1-w^2⟩_V^*× V+⟨γ^*ζ_1-γ^*ζ_2,w^1-w^2⟩_V^*× V
       ≥ m_A w^1 -  w^2_V^2 + ⟨ζ_1 - ζ_2, γ w^1 - γ w^2  ⟩_X^*× X
       ≥ m_A w^1 -  w^2_V^2  -  m_J1γ w^1 - γ w^2  _X^2 
       ≥ (m_A - m_J1 c_γ^2) w^1 -  w^2_V^2.

From (H_s) we see that ∂ℒ_j is strongly monotone. This is equivalent to the fact that ℒ_j is strongly convex (see Theorem 3.4 in <cit.>), which implies that it is strictly convex.


For the proof of (iv), let us fix w ∈ V. From H(A)(c) we obtain

    ℒ_j(w)    =1/2⟨  Aw , w  ⟩_V^*× V +⟨  B (K^kv)_j-1, w  ⟩_V^*× V -  ⟨  f_j,  w ⟩_V^*× V+ J (γ (K^k v)_j-1, γ w)
       ≥1/2 m_A w_V^2 - B (K^kv)_j-1_V^*w_V- f_j_V^*w_V +J (γ (K^k v)_j-1,γ w).

Now, using the Lebourg mean value theorem (cf. Proposition 3.36 in <cit.>), we get that there exist λ∈ (0,1) and η∈∂_2(J∘γ)(γ (K^k v)_j-1, λ w) such that

    J(γ  (K^k v)_j-1, γ w) = ⟨η, w  ⟩_V^*× V + J(γ  (K^k v)_j-1, 0).

Since ∂_2(J ∘γ)(γ (K^k v)_j-1, λ w) ⊆γ^* ∂_2 J(γ (K^k v)_j-1 , λγ w) we have η∈γ^* ∂_2 J(γ (K^k v)_j-1, λγ w ). Then there exists z_1 ∈∂_2 J(γ (K^k v)_j-1 , λγ w ) such that η =  γ^* z_1. Taking z_2 ∈∂_2 J(γ (K^k v)_j-1, 0) and by (<ref>), we obtain

    λ⟨γ^* z_1 -  γ^* z_2, w  ⟩_V^*× V = ⟨ z_1 - z_2, λγ w ⟩_X^*× X≥ - m_J1λγ w _X^2 ≥ - m_J1λ^2 c_γ^2  w _V^2,

and this, along with the fact that λ∈ (0,1), leads to

    ⟨η, w⟩_V^*× V≥ - m_J1c_γ^2w_V^2+⟨ z_2,γ w⟩_X^*× X.

Using H(J)(b), we get for given v

    ∂_2 J(γ (K^k v)_j-1, 0)_X^*≤ c_0 + c_1 γ (K^k v)_j-1_X≤ c.

Hence

    ⟨ z_2, γ w ⟩_X^*× X≥ -|⟨ z_2, γ w ⟩_X^*× X| ≥ - ∂_2 J(γ (K^k v)_j-1, 0)_X^*γ w _X ≥ -cw _V.

Combining (<ref>)–(<ref>) and because J(γ (K^k v)_j-1, 0) is bounded from below for fixed first argument,  we get

    ℒ_j(w)    ≥1/2 m_A w_V^2 - B (K^kv)_j-1_V^*w_V - f_j_V^*w_V
    - m_J1 c_γ^2  w _V^2 -cw _V  - c
       ≥(1/2m_A - m_J1 c_γ^2)  w _V^2 -cw _V - c.

From (H_s) we see that ℒ is coercive.


The problem under consideration reads as follows.

Problem P_opt^h: Find v^hk = {v_j^hk}_j=1^N ⊂ V^h such that for j∈{1,…,N},

    0 ∈   ∂ℒ_j( v_j^hk ).


We are now in a position to prove the existence and uniqueness result for the above optimization problem.

 
Assume H(A), H(B), H(J), H(f), H(u_0) and (H_s). Then Problem P_opt^h has a unique solution.


From Lemma <ref> (i), (iv), we see that the functional ℒ_j is proper, lower semicontinuous
and coercive. This implies that it attains a global minimum. Uniqueness of that
minimum is guaranteed by Lemma <ref> (iii).


Let us conclude results from the previous lemma by the following theorem.

 
Assume H(A), H(B), H(J), H(f), H(u_0) and (H_s). Then Problems P_incl^h and P_opt^h are equivalent,
have a unique solution v^h∈ V^h and this solution satisfies

    v_j^hk_V ≤ c(1+ f_j _V^*),   1≤ j≤ N

with a positive constant c.


Lemma <ref> (ii) implies that every solution to Problem P_opt^h solves Problem P_incl^h.
Using this fact, Lemmas <ref> and <ref>, we see that a unique solution to Problem P_opt^h is also a unique solution to Problem P_incl^h. Because of the uniqueness of the solution to Problem P_incl^h we get that Problems P_incl^h and P_opt^h are equivalent. The estimate in the statement of the theorem follows from Lemma <ref>.


Now let us conclude this section by presenting the following main theorem concerning error estimate of the numerical scheme.

 
Assume that H(A), H(B), H(J), H(f) and (H_s) hold. Let v and v^hk be unique solutions to
Problems P_hvi and P_opt^h, respectively. Additionally assume that v∈ W^1,∞(0,T;V).
Then there exists a constant c>0 such that

    max_1≤ j ≤ Nv_j-v ^hk_j_V^2
       ≤ cmax_1≤ j≤ N inf_w^h∈ V^h{k^2+v_j-w^h_V^2+γ v_j-γ w^h_X
    +|R_j(v_j-w^h)| }
          + c u_0 - u_0^h_V^2,

where a residual quantity is given by

    R_j(w)=⟨ A v_j,w⟩_V^*× V+⟨ B(K v)_j, w⟩_V^*× V-⟨ f_j, w⟩_V^*× V.



Let v be a solution to Problem P_hvi and v^hk be a solution to Problem P_opt^h,
and let us fix  j ∈{1,…,N}.  Then we have

    ⟨  Av_j + B (Kv)_j,  w⟩_V^*× V + J_2^0(γ (Kv)_j, γ v_j;  γ w)
    ≥⟨  f_j,  w⟩_V^*× V  ∀  w ∈ V,

and, as in (<ref>),

    ⟨ Av^hk_j+B(K^k v^hk)_j-1,w^h⟩_V^*× V+J_2^0(γ(K^k v^hk)_j-1,γ v^hk_j; γ w^h)≥⟨  f_j,  w^h⟩_V^*× V  ∀ w^h ∈ V^h.

Taking (<ref>) with w=v^hk_j-v_j and adding to (<ref>) with w^h replaced by w^h - v^hk_j,
after some calculation we obtain for all w^h ∈ V^h,

    ⟨ A v^hk_j - A v_j,  v^hk_j - v_j ⟩_V^*× V≤⟨  A v^hk_j - A v_j,  w^h - v_j  ⟩_V^*× V
         +⟨ B (Kv)_j - B (K^k v^hk)_j-1, v^hk_j - v_j ⟩_V^*× V
         + ⟨ B (K^k v^hk)_j-1 - B (K v)_j, w^h - v_j ⟩_V^*× V
         + J_2^0(γ (K v)_j, γ v_j;  γ v^hk_j - γ v_j) + J_2^0(γ (K^k v^hk)_j-1, γ v^hk_j;  γ w^h - γ v^hk_j)
         + ⟨  A v_j,  w^h - v_j  ⟩_V^*× V + ⟨ B (K v)_j, w^h - v_j ⟩_V^*× V + ⟨ f_j, v_j - w^h ⟩_V^*× V

Using H(A) we have

    m_A  v^hk_j - v_j_V^2 ≤⟨ A v^hk_j - A v_j,  v^hk_j - v_j ⟩_V^*× V.

From the Schwartz inequality, the fact that A∈ℒ(V,V^*) and the Young inequality, we get

    ⟨  A v^hk_j - A v_j,  w^h - v_j  ⟩_V^*× V   ≤ A v^hk_j - A v_j_V^*w^h - v_j_V ≤A_ℒ(V,V^*) v^hk_j - v_j_V w^h - v_j_V 
       ≤ε v^hk_j - v_j_V^2 + c_ε A^2_ℒ(V,V^*) w^h - v_j_V^2.

Analogously for the operator B,

    ⟨ B (Kv)_j - B (K^k v^hk)_j-1, v^hk_j - v_j ⟩_V^*× V≤ L_B^2 c_ε (K^k v^hk)_j-1 - (Kv)_j_V^2 + εv^hk_j - v_j_V^2,

and

    ⟨ B (K^k v^hk)_j-1 - B (K v)_j, w^h - v_j ⟩_V^*× V≤L_B^2/2 (K^k v^hk)_j-1 - (Kv)_j_V^2 + 1/2w^h - v_j_V^2.

Now, from (<ref>) there exits M > 0 such that for all i ∈{1,…,N}

    v_i^hk_V ≤ c(1+ f_i _V^*) ≤ c(1+max_1≤ i ≤ N f_i _V^*) ≤ M,

and, since kN = T, we have

    (K^k v^hk)_j-1_V ≤  k ∑_i=1^j-1v_i^hk_V  + u_0^h_V ≤ kN M + u_0^h_V ≤ T M + c.

From these estimations combined with Proposition 3.23(ii) in <cit.> and by assumption H(J)(b) we have

    |J_2^0(γ (K^k v^hk)_j-1, γ v^hk_j;  γ w^h - γ v_j )|
      
       ≤(c_0 +  c_1c_γ (K^k v^hk)_j-1_V + c_2c_γ v^hk_j_V) γ w^h -  γ v_j _X
       ≤ c γ w^h -  γ v_j _X.

 We now use assumption H(J)(c) and the Young inequality to obtain
 
    J_2^0(γ (K^k v^hk)_j-1, γ v^hk_j;  γ v_j - γ v^hk_j )+ J_2^0(γ (K v)_j, γ v_j;  γ v^hk_j - γ v_j )
    
     
           ≤ m_J1c_γ^2   v^hk_j -  v_j _V^2 + m_J2^2  c_ε c_γ^4   (K^k v^hk)_j-1 - (K v)_j _V^2 + ε   v^hk_j - v_j _V^2.

Hence, from subadditivity of generalized directional derivative (Proposition 3.23(i) in <cit.>) we get

    J_2^0(γ(K^k v^hk)_j-1,γ v^hk_j;  γ w^h - γ v^hk_j ) + J_2^0(γ (K v)_j, γ v_j;  γ v^hk_j - γ v_j )
         ≤ J_2^0(γ (K^k v^hk)_j-1, γ v^hk_j;  γ w^h - γ v_j ) + J_2^0(γ (K^k v^hk)_j-1, γ v^hk_j;  γ v_j - γ v^hk_j ) 
           + J_2^0(γ (K v)_j, γ v_j;  γ v^hk_j - γ v_j )
         ≤ c γ w^h -  γ v_j _X + m_J1c_γ^2   v^hk_j -  v_j _V^2
           + m_J2^2  c_ε c_γ^4   (K^k v^hk)_j-1 - (K v)_j _V^2 + ε   v^hk_j - v_j _V^2.

Using inequalities (<ref>)–(<ref>) in (<ref>), we obtain

    (m_A - m_J1c_γ^2 - 3ε)  v^hk_j - v_j_V^2
           ≤  c_ε (K^k v^hk)_j-1 -  (K v)_j _V^2 + c_εv_j - w^h _V^2 + c γ w^h -  γ v_j _X
              + ⟨  A v_j,  w^h - v_j  ⟩_V^*× V + ⟨ B (K v)_j, w^h - v_j ⟩_V^*× V + ⟨  f_j,  v_j - w^h  ⟩_V^*× V.

We now calculate
 
    (K^k v^hk)_j-1 - (K v)_j _V^2 ≤  2( (K^k v)_j-1 - (K v)_j _V^2 +  (K^k v^hk)_j-1 - (K^k v)_j-1_V^2 ) ,

where

    (K^k v)_j = k ∑_i=1^j v_i + u_0^h,     v ∈ C([0,T];V).

Using the triangle inequality, the elementary inequality (<ref>),
and the Jensen inequality (Lemma <ref>) with I = [t_i-1,t_i] (cf. (<ref>)), since kN=T, we obtain
 
    (K v)_j - (K^k v)_j-1_V^2 =  ∫_0^t_j v(s)   ds + u_0 - (k ∑_i=1^j-1 v_i + u_0^h) _V^2
       =  ∑_i=1^j-1∫_t_i-1^t_i  (v(s) - v_i)   ds +  ∫_t_j-1^t_j  v(s)  ds  + u_0 - u_0^h_V^2
       ≤(∑_i=1^j-1∫_t_i-1^t_i v(s) - v_i _V  ds + ∫_t_j-1^t_jv(s)_V  ds + u_0 - u_0^h_V)^2
       ≤ cN ∑_i=1^j-1(∫_t_i-1^t_i v(s) - v_i _V  ds)^2 + c (∫_t_j-1^t_jv(s)_V  ds)^2 + cu_0 - u_0^h_V^2
       ≤ cT ∑_i=1^j-1∫_t_i-1^t_i v(s) - v_i _V^2   ds + c(∫_t_j-1^t_jv(s)_V  ds)^2 + cu_0 - u_0^h_V^2.

Since v ∈ W^1,∞(0,T;V), we obtain for all i ∈{1,…,N}, s ∈ (t_i-1,t_i]

    v(s) -  v_i _V = ∫_t_i^s v'(τ)  dτ_V≤∫_t_i^t_i+1v'(τ)_V dτ≤ kv'_L^∞(0,T;V)≤ kv_W^1,∞(0,T;V),

and

    ∫_t_j-1^t_jv(s)_V ds ≤ kv_L^∞(0,T;V)≤ kv_W^1,∞(0,T;V).

Hence, from the above inequalities, we get

    (K v)_j - (K^k v)_j-1_V^2
       ≤ ck^2 ∑_i=1^j-1∫_t_i-1^t_iv_W^1,∞(0,T;V)^2  ds+ ck^2v_W^1,∞(0,T;V)^2
    + cu_0 - u_0^h_V^2 
       ≤ cTk^2 + ck^2 + cu_0 - u_0^h_V^2.

And now, from the triangle inequality

    (K^k v^hk)_j-1 - (K^k v)_j-1_V^2
       = k∑_i=1^j-1 v_i^hk + u_0^h - (k∑_i=1^j-1 v_i + u_0^h)_V^2
       ≤(k∑_i=1^j-1v_i^hk- v_i_V)^2≤ Tk∑_i=1^j-1 v_i^hk - v_i_V^2.

Combining (<ref>), (<ref>) and (<ref>), we get

    (K^k v^hk)_j-1 - (K v)_j _V^2 ≤ ck^2 + ck∑_i=1^j-1 v_i^hk-v_i_V^2 +cu_0-u_0^h_V^2.

Returning to (<ref>), we obtain

    (m_A - m_J1c_γ^2 - 3ε)  v^hk_j - v_j_V^2
       ≤ c_ε k^2 + c_εv_j - w^h _V^2 + c γ w^h -γ v_j_X+|R_j(v_j-w^h)|
          + c_εu_0 - u_0^h_V^2 + c_ε k∑_i=1^j-1 v_i^hk - v_i_V^2.

Taking sufficiently small ε, from smallness assumption (H_s) we get

    v^hk_j - v_j_V^2    ≤ c k^2 + cv_j - w^h _V^2 + c γ w^h -  γ v_j _X 
         + c |R_j(v_j-w^h)| + cu_0 - u_0^h_V^2 + c k∑_i=1^j-1 v_i^hk - v_i_V^2.

From discrete version of the Gronwall inequality (Lemma <ref>) with

    e_j =  v^hk_j - v_j_V^2, 
        g_j = c k^2 + cv_j - w^h _V^2 + c γ w^h -γ v_j _X + c |R_j(v_j-w^h)| + cu_0 - u_0^h_V^2,

we obtain the required estimation.





§ APPLICATION TO A CONTACT PROBLEM
 

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.

A viscoelastic body occupies domain Ω⊂ℝ^d, where d = 2,3 in application.
We assume that its boundary Γ is divided into three disjoint measurable parts:
Γ_D, Γ_C, Γ_N, where Γ_D has a positive measure.
Additionally Γ is Lipschitz continuous, and therefore the outward 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 inclusion. Forces and contact conditions may be time dependent. We are interested in finding  body displacement in the time interval [0,T], with T>0.

Let us denote by “·” and · the scalar product and the Euclidean norm in ℝ^d or 𝕊^d, respectively, where  𝕊^d = ℝ^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(Ω;ℝ^d) and H^1(Ω)^d = H^1(Ω;ℝ^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 quasistatic mechanical contact problem.

Problem P^M: Find a displacement field uΩ× [0,T] →ℝ^d and a stress field σΩ× [0,T] →𝕊^d such that for all t ∈[0,T],

    σ= 𝒜(ε(u'))+ℬ(ε(u))      in Ω, 
    Div σ + f_0 = 0      in Ω, 
    u = 0        on Γ_D, 
    σν = f_N        on Γ_N, 
    
    -σ_ν = g_ν(u_ν)         on Γ_C, 
    
    -σ_τ∈ g_τ (u_ν) ∂ j_τ(u_τ')         on Γ_C, 
    u(0) = u_0        in Ω.


Here, equation (<ref>) represents an viscoelastic constitutive law, 𝒜 is a viscosity operator and
ℬ is an elasticity operator. Equilibrium equation (<ref>) reflects the fact that problem is quasistatic.
Equation (<ref>) represents clamped boundary condition on Γ_D and (<ref>) reflects the forces acting
on Γ_N. Relation (<ref>) describes the response of the foundation in normal direction, whereas the friction
is modeled by inclusion (<ref>), where j_τ is a given superpotential, and g_τ is a given friction bound.
Finally, equation (<ref>) represents the initial condition with the initial displacement u_0.

We use the following Hilbert spaces

    H = L^2(Ω; ℝ^d),   ℋ = L^2(Ω; 𝕊^d),
       H_1 = {u∈ H | ε(u) ∈ℋ},
      ℋ_1 = {σ∈ℋ | Div σ∈ H},
       V = {v∈ H_1 | v = 0  on Γ_D}

endowed with the inner products

    (u,v)_H=∫_Ω u_iv_i dx,  (σ,τ)_ℋ=∫_Ωσ_ijτ_ij dx, 
       (u,v)_H_1 = (u,v)_H +(ε(u),ε(v))_ℋ,   
    (σ,τ)_ℋ_1=(σ,τ)_ℋ
    +(Div σ,Div τ)_H,
       (u,v)_V = (ε(u),ε(v))_ℋ

and corresponding norms _X with X being H, ℋ, H_1, ℋ_1, V.
The fact that space V equipped with the norm ·_V is complete follows from Korn's inequality, and its application is allowed because we assumed that meas(Γ_D) > 0.
We consider the trace operator γ V → L^2(Γ_C)^d=X.
By the Sobolev trace theorem we know that γ∈ℒ(V, X) with the norm equal to c_γ.

Now we present the hypotheses on data of Problem P^M.

H(𝒜):   𝒜Ω×𝕊^d →𝕊^d satisfies

  
  (a)
     𝒜(x,τ) = (a_ijkh(x)τ_kh)
     for all τ∈𝕊^d, a.e. x∈Ω, a_ijkh∈ L^∞(Ω),
  
  (b)
    𝒜(x,τ_1) ·τ_2 = τ_1 ·𝒜(x,τ_2) for all τ_1, τ_2 ∈𝕊^d, a.e. x∈Ω,
  
  (c)
     there exists m_𝒜>0 such that 𝒜(x,τ) ·τ≥ m_𝒜τ^2 for all τ∈𝕊^d, a.e. x∈Ω.



H(ℬ):   ℬΩ×𝕊^d →𝕊^d satisfies

  
  (a)
     ℬ(·,τ) is measurable on Ω for all
        τ∈𝕊^d, ℬ(·,τ)∈ℋ,
  
  (b)
    there exists L_ B > 0 s.t. B(,τ_1)- B(,τ_2)≤ L_ Bτ_1- τ_2 for all τ_1,τ_2∈𝕊^d, a.e. ∈Ω,
  
  (c)
     there exists m_ℬ>0 such that (ℬ(x,τ_1)-ℬ(x,τ_2))· (τ_1-τ_2) ≥ m_ℬτ_1-τ_2^2 for all τ_1, τ_2 ∈𝕊^d, a.e. x∈Ω.


H(j_τ):   j_τΓ_C ×ℝ^d→ℝ satisfies

  
  (a)
    j_τ(·, ξ) is measurable on Γ_C for all ξ∈ℝ^d and there exists e∈ L^2(Γ_C)^d such that j_τ(·,e(·))∈ L^1(Γ_C),
  
  (b)
    there exists c_τ>0 such that 

    |j_τ(x, ξ_1) - j_τ(x, ξ_2)| ≤ c_τξ_1 - ξ_2  for all ξ_1, ξ_2 ∈ℝ^d, a.e. x∈Γ_C,
  
  (c)
    there exists α_τ≥ 0 such that  

    j_τ^0(x,ξ_1;ξ_2-ξ_1) + j_τ^0(x,ξ_2;ξ_1-ξ_2)≤α_τξ_1-ξ_2^2

     for all ξ_1, ξ_2 ∈ℝ^d, a.e. x∈Γ_C.


H(g):   g_ιΓ_C×ℝ→ℝ, ι∈{ν, τ}, satisfies

  
  (a)
    g_ι(·, r) is measurable on Γ_C for all r ∈ℝ,
  
  (b)
    there exists g_ι > 0 such that 0 ≤ g_ι(x, r) ≤g_ι for all r ∈ℝ, a.e. x∈Γ_C,
  
  (c)
    there exists L_g_ι>0 such that 

     |g_ι(x, r_1) - g_ι(x, r_2)| ≤ L_g_ι |r_1 - r_2|  for all r_1, r_2 ∈ℝ, a.e. x∈Γ_C.


(H_0):   f_0 ∈ C([0,T]; L^2(Ω)^d),   f_N ∈ C([0,T]; L^2(Γ_N)^d).

We remark that condition H(j_τ)(b) is equivalent to the fact that j_τ(x,·) is 
Lipschitz continuous and there exists c_τ≥ 0 such that ∂ j_τ(x,ξ)≤ c_τ for all ξ∈ℝ^d and a.e. x∈Γ_C. Moreover, condition H(g)(b) is enough to obtain presented mathematical results, but from mechanical point of view we should additionally assume that g_ι(r)=0 for r≤ 0. This corresponds to the situation when body is separated from the foundation and normal response of the foundation and friction force vanishes.

Using the standard procedure, the Green formula and the definition of generalized subdifferential, we obtain a weak formulation of Problem P^M in the form of hemivariational inequality.


Problem P^M_hvi:  Find a velocity v [0,T]→ V satisfying
for all w∈ V and  t∈ [0,T]

    ⟨ Av(t)+B(Kv)(t),w⟩_V^*× V+∫_Γ_Cj_2^0
    (γ(Kv)(t),γv(t);γw) da≥⟨f(t),w⟩_V^*× V.


Here, operators A,B  V → V^*, K L^2(0,T;V) → L^2(0,T;V) and f(t) ∈ V^* are defined for all
w,v∈ V, t ∈ [0,T] as follows:

    ⟨ Av,w⟩_V^*× V=(𝒜(ε(v)),ε(w))_ℋ,
       ⟨ Bv,w⟩_V^*× V=(ℬ(ε(v)),ε(w))_ℋ,
       (Kv)(t) = ∫_0^t v(t)   ds + u_0,
       ⟨f(t),w⟩_V^*× V=∫_Ωf_0(t)·w dx
    +∫_Γ_Cf_N(t)·γw  da,

and jΓ_C×ℝ^d ×ℝ^d →ℝ be defined for all ξ, η∈ℝ^d and x∈Γ_C by

    j(x, η, ξ) =  g_ν(x, η_ν) ξ_ν + g_τ(x, η_ν)  j_τ(x, ξ_τ).

It is easy to check that under assumptions H(𝒜), H(ℬ) and (H_0), the operators A and B and the functional f satisfy H(A), H(B) and H(f), respectively.
We also define the functional J  L^2(Γ_C)^d × L^2(Γ_C)^d →ℝ for all w, v∈ L^2(Γ_C)^d by

    J(w, v) =  ∫_Γ_C j(x, w(x), v(x))  da.


Below we present some properties of the functional J.

 
Assumptions H(j_τ) and H(h) imply that functional J defined by (<ref>)–(<ref>) satisfies H(J).



We first observe that from H(j_τ)(a),(b) and H(h)(a),(c) the function j(·,η,ξ) is measurable on Γ_C, there exist e_1, e_2 ∈ L^2(Γ_C)^d such that j(·, e_1(·), e_2(·))∈ L^1(Γ_C), j(x, ·, ξ) is continuous and j(x, η, ·) is locally Lipschitz. Moreover, by H(j_τ)(b) and H(h)(b) we easily conclude

    ∂_2 j(x, η,ξ)   ≤ g_ν(x,η_ν) + g_τ(x,η_ν) ∂ j_τ (x,ξ_τ)≤g_ν+g_τ  c_τξ.

Applying similar procedure to one presented in the proof of Corollary 4.15 in <cit.>, we obtain that functional J is well defined, locally Lipschitz with respect to the second variable and the growth condition H(J)(b) holds with c_0= √(2 meas(Γ_C)) g_ν, c_1=√(2) g_τ  c_τ and  c_2=0.

To prove H(J)(c), we take η_i, ξ_i ∈ℝ^d, i = 1, 2, and by the sum rules
(cf. Proposition 3.35 in <cit.>) and from H(j_τ)(b),(c) and H(h)(b),(c), we obtain

    j_2^0(x,η_1,ξ_1;ξ_2-ξ_1)+j_2^0(x,η_2,ξ_2;ξ_1- ξ_2)
           ≤ g_ν(x,η_1ν)(ξ_2ν- ξ_1ν) - g_ν(x,η_2ν)(ξ_1ν- ξ_2ν)  
              + g_τ(x,η_1ν)(j_τ^0(x,ξ_1τ;
    ξ_2τ-ξ_1τ) + j_τ^0(x,ξ_2τ; ξ_1τ-ξ_2τ)) 
              + (g_τ(x, η_2ν)-g_τ( x,η_1ν)) j_τ^0(x,ξ_2τ; ξ_1τ -  ξ_2τ) 
           ≤g_τ α_τ ξ_1 - ξ_2^2
      + (L_g_ν + L_g_τ c_τ) η_1 - η_2 ξ_1 - ξ_2.

Consequently, since

    J_2^0(w, v; z) ≤∫_Γ_C j_2^0 (x, w(x), v(x); z(x)) da

(cf. Corollary 4.15 (iii) in <cit.>), we have

    J_2^0(w_1, v_1; v_2- v_1) + J_2^0(w_2, v_2; v_1- v_2) 
           ≤∫_Γ_C(g_τα_τv_1(x)-v_2(x)^2
    +(L_g_ν+L_g_τc_τ)w_1(x)-w_2(x) v_1(x)-v_2(x))  da.

Hence, by the Hölder inequality, we obtain H(J)(c) with m_α = g_τα_τ and m_L =L_g_ν + L_g_τc_τ.


With the above properties, we can check that assumptions of Theorems <ref>, <ref> and <ref>
are satisfied. We can employ previously presented abstract framework and conclude that Problem P^M_hvi
has a unique solution.

We now turn to the numerical solution of P^M_hvi.  For simplicity, we assume Ω is a polygonal/polyhedral
domain, and express the three parts of the boundary, Γ_D, Γ_N and Γ_C  as unions of closed
flat components with disjoint interiors:

    Γ_D =∪_i=1^i_DΓ_D,i,  Γ_N =∪_i=1^i_NΓ_N,i,
      Γ_C =∪_i=1^i_CΓ_C,i.

Let { T^h}_h be a regular family of finite element partitions of Ω into
triangular/tetrahedral elements, compatible with the partition of the boundary ∂Ω into
Γ_D,i for 1≤ i≤ i_D, Γ_N,i for 1≤ i≤ i_N, and Γ_C,i for 1≤ i≤ i_C,
i.e. if the intersection of one side/face of an element with one of these sets has a positive measure,
then the side/face lies entirely in that set.  Here h→0 denotes the finite element mesh-size.
Corresponding to the partition T^h, we introduce the linear finite element space

    V^h = {^h∈ C(Ω)^d |^h|_T∈ℙ_1(T)^d
        ∀ T∈ T^h, ^h=  on Γ_D}.

As in the previous section, given a positive integer N, define the step-size k=T/N and the nodes t_n=nk,
0≤ n≤ N.   We will assume

    u_0∈ H^2(Ω)^d.

Let u^h_0∈ V^h be the interpolant of u_0 in V^h.  Then (<cit.>)

    u_0-u^h_0_V≤ c h.


Then we introduce the following discretized version of Problem P_hvi.

Problem P_hvi^h:  Find a velocity v^hk={v^hk_j}_j=1^N⊂ V^h
such that for 1≤ j≤ N,

    ⟨ Av^hk_j+B(K^k v^hk)_j-1,w^h⟩_V^*× V+∫_Γ_Cj_2^0
    (γ(K^k v^hk)_j-1,γv^hk_j;γw^h) da≥⟨f_j,w^h⟩_V^*× V  ∀ w^h∈ V^h.


Similar to Problem P^M_hvi, under the stated assumptions, Problem P_hvi^h has a unique solution.
Assume ∈ W^1,∞(0,T;V).  It is easy to see from the proof of Theorem <ref> that the
inequality (<ref>) remains valid when J^0_2(u,v;w) is replaced by
∫_Γ_C j^0_2(u,v;w) da.  Thus, we have a constant c>0 such that

    max_1≤ j ≤ N_j- ^hk_j_V^2
       ≤ cmax_1≤ j≤ N inf_ w^h∈ V^h{k^2+_j-^h_V^2
    +γ_j-γ^h_L^2(Γ_C)^d+|R_j(_j-^h)| }
          + c _0 - _0^h_V^2,

where

    R_j()=⟨ A_j+ B(K)_j,⟩_V^*× V-⟨_j,⟩_V^*× V.




Assume H(𝒜), H(ℬ), H(j_τ), H(g), (H_0), (H_s), and (<ref>),
and assume the solution regularity

    ∈ W^1,∞(0,T;V),  ∈ C([0,T]; H^2(Ω)^d),  ∈ C([0,T]; L^2(Γ_C)^d),
    
       |_Γ_C,i∈ C([0,T]; H^2(Γ_C,i)^d), 1≤ i≤ i_C.

Then, for the solution v to Problem P^M_hvi and the solution v^hk to Problem P_hvi^h there exists a constant c>0 such that

    max_1≤ j ≤ N_j- ^hk_j_V≤ c(k+h).



We bound the residual term defined by (<ref>) using similar procedure to one described in <cit.>.
Let w∈ C^∞(Ω)^d be arbitrary with w=0 on Γ_D∪Γ_C in
the inequality (<ref>). Let us fix j ∈{1,…,N}. We can derive the equality

    ⟨ Av_j+B(Kv)_j, w⟩_V^*× V  = ⟨f_j, w⟩_V^*× V  ∀ w∈ C^∞(Ω)^d,  w=0  on Γ_D∪Γ_C.

Then it is possible to deduce that

    Div (𝒜(ε(u'_j)) + ℬ(ε(u_j)))
        + (f_0)_j = 0   in the sense of distributions.

Since (f_0)_j ∈ L^2(Ω)^d, we have

    Div (𝒜(ε(u'_j)) + ℬ(ε(u_j)))
        + (f_0)_j = 0  a.e. in Ω.

It is also possible to deduce that
 
    σ_jν = (f_N)_j   a.e. on Γ_N.

For details, see <cit.>.  We multiply equation (<ref>) by w∈ V to obtain

    ∫_Γσ_jν·w da-∫_Ω[𝒜(ε(u'_j))
    + ℬ(ε(u_j))]·ε(w) dx
    +∫_Ω(f_0)_j·w dx=0.

Using the homogeneous Dirichlet boundary condition of w on Γ_D and the traction boundary condition
given by (<ref>) we have

    ∫_Ω[𝒜(ε(u'_j))
    + ℬ(ε(u_j))]·ε(w) dx
    =∫_Γ_C_j· da + ∫_Ω(f_0)_j·w dx+∫_Γ_N(f_N)_j·w da.

Thus, for R_j() defined by (<ref>), we have

    R_j()=  ∫_Γ_Cσ_jν·w da≤ c w_L^2(Γ_C)^d.

We denote by Π^h u∈ V^h the finite element interpolant of u.
From (<ref>), (<ref>) and (<ref>),  we get

    max_1≤ j ≤ N_j- ^hk_j_V^2
    ≤ cmax_1≤ j≤ N{k^2+_j-Π^h _j_V^2
    +_j-Π^h _j_L^2(Γ_C)^d}+ c h^2.

By the standard finite element interpolation error bounds (<cit.>), due to the solution regularity
(<ref>) and (<ref>), we have

    _j-Π^h _j_V ≤ c h,
       _j-Π^h _j_L^2(Γ_C)^d≤ c h^2.

Using these bounds in (<ref>), we obtain the error estimate (<ref>).





§ NUMERICAL RESULTS
 

In this section, we report computer simulation results on a numerical example. We apply numerical scheme P_opt^h to approximate solution of P^M_hvi. The linear finite element space V^h based on uniform triangular partition of Ω and the uniform partition of the time interval [0,1] with the time step size k=1/N for a positive integer N are used. In order to minimize not necessarily differentiable functional ℒ we use Powell's conjugate direction method. This method does not require the assumption that optimized function is differentiable. Other, more refined nonsmooth optimization algorithms described for instance in <cit.>, could also be adapted.

We set d=2 and consider a square-shaped set Ω=(0,1)×(0,1) presented in Figure <ref> with the following partition of the boundary

    Γ_D = {0}× [0,1],   Γ_N = ([0,1] ×{1}) ∪ ({1}× [0,1]),
      Γ_C = [0,1] ×{0}.



We employ the Kelvin-Voigt type short memory viscoelastic law for the isotropic body. The viscosity operator 𝒜 and the elasticity operator ℬ are defined by

    𝒜(τ) = 2ϕτ + ξ(τ)I,    τ∈𝕊^2,
    ℬ(τ) = 2ητ + λ(τ)I,    τ∈𝕊^2.


Here I denotes the identity matrix,  denotes the trace of the matrix, λ and η
are the Lamé  coefficients, whereas ϕ and ξ represent the viscosity coefficients,  λ, η, ϕ, ξ > 0.
In our simulations, we choose ϕ = ξ = 2, λ = η = 4, T=1.


We first demonstrate the effect of some input data on the deformation of the body.
In all cases, we show the shape of the body at final time t=1, as well as the contact interface forces on Γ_C. The numerical solutions correspond
to the time step size 1/32 and the boundary Γ_C of the body divided into 32 equal parts.
Let us take the following data

    u_0(x) = (0,0),   x∈Ω,
       g_ν(x, η)= {[           0, η∈ (-∞,  0),;        30 η, η∈ [0, 0.1),;           3,  η∈ [0.1,∞),;              ]. x∈Γ_C,
       g_τ(x, η)= g_ν(x, η),   η∈ℝ,  x∈Γ_C,
       j_τ(x, ξ) = -0.3 e^-ξ + 0.7 ξ ,   ξ∈ℝ^2, x∈Γ_C,
       f_N(x,t) = (0,0),   x∈Ω, t ∈ [0,T],
       f_0(x,t) = (-2.5,-0.5),   x∈Ω, t ∈ [0,T].



We note that function j_τ, based on Example 7.26 in <cit.>, is nondifferentiable and nonconvex.
Our aim is to investigate reaction of the body to various modifications of input data.


In Figure <ref> we present output obtained without any modifications. We push the body down and to the left
with force f_0. As a result the body penetrates the foundation, but frictional forces restrict its movement
in proximity of Γ_C. Next, we modify the function g_ν to be given by

    g_ν(x, η) = {[          0, η∈ (-∞, 0),;      200 η, η∈ [0,0.1),;         20, η∈ [0.1,∞),;             ]. x∈Γ_C.

In Figure <ref> we observe that this modification models more rigid foundation by increasing its response in normal direction on Γ_C. The result is decreased penetration of the foundation u_ν. The friction also decreases, due to influence of function g_τ which depends on u_ν.
Now we return to original data and only change the direction of force f_0 to the following

    f_0(x,t) = (2.5,  -0.5),   x∈Ω, t ∈ [0,T].

In Figure <ref> we observe that the body displaces in opposite direction to previous examples.
Because of frictional forces it moves to the right more in the higher part than in the lower part.
The penetration of the foundation and friction increase as we get closer to bottom right corner of the body.
In the last experiment we once more return to original data and modify the function g_τ as follows

    g_τ(x, η) = {[                            0,            η∈ (-∞, 0), x∈Γ_C,;                            0,   η∈ [0,∞), x∈ [0.5, 1] ×{0},;                         30 η, η∈ [0,0.1), x∈ [0, 0.5) ×{0},;                            3, η∈ [0.1,∞), x∈ [0, 0.5) ×{0}.;                               ].

In this case Γ_C is divided into two parts, and right part is covered in grease. In Figure <ref> we see that contact of left part of the body with the foundation creates friction, whereas contact of right part is frictionless.














In order to illustrate the error estimate obtained in Section <ref>, we present the empirical convergence orders of the numerical method on the model problem. We take the following data

    u_0(x) = (0,0),   x∈Ω,
       g_ν(x, η)= {[           0, η∈ (-∞,  0),;        60 η, η∈ [0, 0.1),;           6,  η∈ [0.1,∞),;              ]. x∈Γ_C,
       g_τ(x, η)= {[           0, η∈ (-∞,  0),;       120 η, η∈ [0, 0.1),;          12,  η∈ [0.1,∞),;              ]. x∈Γ_C,
       j_τ(x, ξ) = ξ,   ξ∈ℝ^2, x∈Γ_C,
       f_N(x,t) = (-0.2,-0.2),   x∈Ω, t ∈ [0,T],
       f_0(x,t) = (-1,-0.4),   x∈Ω,  t ∈ [0,T].

We have v_V≐ 0.06738.
We present a comparison of numerical errors v - v^hk_V computed for a sequence of solutions to discretized problems.
We use a uniform discretization of the problem domain and time interval according to the spatial
discretization parameter h and the time step size k, respectively. The boundary Γ_C of Ω
is divided into 1/h equal parts. The numerical solution corresponding to h = 1/256 and k = 1/256 is taken as the “exact” solution v.

For the first experiment, we fix h = 1/256 and start with k = 1/2, which is successively halved.
The results are presented in Table <ref> and Figure <ref>, where the dependence of the relative error estimates v  - v^hk_V / v_V  with respect to k are plotted on a log-log scale.  A first order convergence can be observed for the numerical solutions of the displacement.

For the second experiment, we fix k = 1/256 and start with h = 1/2, which is also successively halved.
The results are presented in Table <ref> and Figure <ref>. Again, a first order convergence can be observed.


Acknowledgments

The project has received funding from the European Union's Horizon 2020 Research and Innovation Programme under the Marie Sklo­do­wska-Curie grant agreement no. 823731 CONMECH.
It is supported by the projects financed by the Ministry of Science and Higher Education of Republic of Poland under
Grants Nos. 4004/GGPJII/H2020/2018/0 and 440328/PnH2/2019.



00

BKM
    A. Bagirov, N. Karmitsa, M. M. Mäkelä,
    Introduction to Nonsmooth Optimization Theory, Practice and Software,
    Springer International Publishing, 2014.

BBK
    M. Barboteu, K. Bartosz, P. Kalita,
    An analytical and numerical approach to a bilateral contact problem with nonmonotone friction,
    International Journal of Applied Mathematics and Computer Science, 23(2) (2013), 263–276.

BHM
    M. Barboteu, W. Han, S. Migórski,
    On numerical approximation of a va­ria­tio­nal–he­mi­va­ria­tio­nal inequality modeling contact problems for locking materials,
    Computers and Mathematics with Applications, 77 (2018), 898-1221.

BS2008
    S.C. Brenner,  L.R. Scott,
    The Mathematical Theory of Finite Element Methods, Third Edition,
    Springer-Verlag, New York, 2008.

Ci1978
    P. G. Ciarlet,
    The Finite Element Method for Elliptic Problems,
    North Holland, Amsterdam, 1978.

C
    F. H. Clarke,
    Optimization and Nonsmooth Analysis,
    Wiley Interscience, New York, 1983.

FLG
    L. Fan, S. Liu, S. Gao,
    Generalized monotonicity and convexity of non-differentiable functions,
    Journal of Mathematical Analysis and Applications, 279 (2003), 276–289.

H
    W. Han,
    Numerical analysis of stationary variational-hemivariational inequalities with applications in contact mechanics,
    Mathematics and Mechanics of Solids, 23 (2018), 279–293.

HSM
    W. Han, S. Migórski, M. Sofonea,
    A class of variational-hemivariational inequalities with applications to frictional contact problems,
    SIAM Journal on Mathematical Analysis, 46(6) (2014), 3891–3912.

HS
    W. Han, M. Sofonea,
    Quasistatic Contact Problems in Viscoelasticity and Viscoplasticity,
    American Mathematical Society and International Press, 2002.

HS2
    W. Han, M. Sofonea,
    Numerical analysis of hemivariational inequalities in contact mechanics,
    Acta Numerica, 28 (2019), 175–286.

HSB
    W. Han, M. Sofonea, M. Barboteu,
    Numerical analysis of elliptic hemivariational inequalities,
    SIAM Journal on Numerical Analysis, 55(2) (2017), 640–663.

HSD
    W. Han, M. Sofonea, D. Danan,
    Numerical analysis of stationary va­ria­tio­nal–he­mi­va­ria­tio­nal inequalities,
    Numerische Mathematik, 139(3) (2018), 563–592.

HMP
    J. Haslinger, M. Miettinen, P.D. Panagiotopoulos,
    Finite Element Method for Hemivariational Inequalities. Theory, Methods and Applications,
    Kluwer Academic Publishers, Boston, 1999.

JO2
    M. Jureczka, A. Ochal,
    A nonsmooth optimization approach for hemivariational inequalities with applications in Contact Mechanics,
    Applied Mathematics and Optimization (2019),
    doi.org/10.1007/s00245-019-09593-y.

JO1
    M. Jureczka, A. Ochal,
    Numerical analysis and simulations of contact problem with wear,
    Computers and Mathematics with Applications 77(11), (2018), 2980–2988.

KO
    N. Kikuchi, J. T. Oden,
    Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods,
    SIAM, Philadelphia, 1988.

MH
    M. Miettinen, J. Haslinger,
    Finite element approximation of vector-valued hemivariational problems,
    Journal of Global Optimization, 10(1) (1997), 17–35.

MOS
    S. Migórski, A. Ochal, M. Sofonea,
    Nonlinear Inclusions and Hemivariational Inequalities. Models and Analysis of Contact Problems,
    Advances in Mechanics and Mathematics, vol. 26, Springer, 2013.

P
    P.D. Panagiotopoulos,
    Hemivariational Inequalities, Applications in Mechanics and En­gi­nee­ring,
    Springer-Verlag, 1993.

SM
    M. Sofonea, S. Migórski,
    Variational-Hemivariational Inequalities with Applications,
    Chapman and Hall/CRC Press, Boca Raton-London, 2018.

W
    P. Wriggers,
    Computational Contact Mechanics,
    Wiley, Chichester, 2002.



