Numerical analysis of doubly-history dependent variational inequalities in contact mechanics

This paper is devoted to numerical analysis of doubly-history dependent variational inequalities in contact mechanics. A fully discrete method is introduced for the variational inequalities, in which the doubly-history dependent operator is approximated by repeated left endpoint rule and the spatial variable is approximated by the linear element method. An optimal order error estimate is derived under appropriate solution regularities, and numerical examples illustrate the convergence orders of the method.


Introduction
Contact phenomena are commonly seen in engineering applications and in daily life. Due to the complicated mechanism involved in contact processes, mathematical models of contact problems are in the form of variational inequalities or hemivariational inequalities. An extensive amount of references exist on modeling, mathematical analysis, numerical solution, and engineering applications of contact problems. A few representative references in this aspect include [5,7,9,10,18,19] on variational inequalities arising in contact mechanics and [6,11] on related hemivariational inequalities.
Mathematical models of certain contact problems contain so-called history-dependent operators that reflect the dependence of the models on the history of physical quantities, usually through an integration of the physical quantities with respect to the temporal variable. The notion of a family of history-dependent (quasi)variational inequalities was first introduced in [15]. For mathematical analysis of history-dependent variational inequalities and hemivariational inequalities in contact problems, see [12,13,17], and for the numerical solution of the problems, see [8,20,21].
In this paper, we consider the numerical solution of a new kind of history dependent variational inequalities. The major novelty of this paper is the numerical treatment of a doubly-history dependent operator in the form of a repeated time integration. The rest of the paper is organized as follows. In Sect. 2, we present a variational inequality in viscoelas- tic contact. Solution existence and uniqueness of the variational inequality are shown in [13] through an application of a result on an abstract sweeping process. We provide an elementary proof of the same result that is more accessible to the reader. In Sect. 3, we introduce and study a fully discrete scheme for the variational inequality; in particular, we derive an optimal order error estimate for the numerical solution under appropriate solution regularity assumptions. In Sect. 4, we report simulation results on numerical examples to illustrate the performance of the method with an emphasis on numerical convergence orders.

A variational inequality in viscoelastic contact
We follow [13] and consider a quasi-static contact problem for viscoelastic material. Let ⊂ R d be a Lipschitz domain, representing the configuration of the viscoelastic body. The boundary = ∂ is split into three measurable parts 1 , 2 , and 3 ; the portion 3 is further split into two parts: 3,1 and 3,2 where different contact conditions will be specified. We assume meas( 1 ) > 0 and meas( 3,1 ) + meas( 3,2 ) > 0. The latter assumption allows the case where one of the two subsets 3,1 and 3,2 is empty, and then the corresponding contact condition below is suppressed from the problem. The body is subject to the action of volume forces of a total density f 0 in and surface tractions of a total density f 2 on 2 , and it is fixed on 1 . We assume a frictionless contact with unilateral constraint in the velocity variable on 3,1 , and a version of Coulomb's law of dry friction on 3,2 . Let [0, T] be the time interval of interest. We use a prime to indicate the derivative with respect to the temporal variable.
For a vector v ∈ R d , we write v = (v i ) with its components v i ∈ R, 1 ≤ i ≤ d. We adopt the summation convention on a repeated index. Over the space R d , we use the canonical inner product and norm defined by We use the symbol S d to denote the space of symmetric matrices of order d.
The canonical inner product and norm on the space S d are The unit outward normal vector exists a.e. on , and it is denoted as ν. For a vector field v defined on , its normal and tangential components Similarly, for a stress tensor σ defined on , the normal and tangential components are σ ν = (σ ν) · ν, σ τ = σ νσ ν ν.
The classical formulation of the contact problem is to find a displacement field u : × [0, T] and a stress field σ : × [0, T] such that, for t ∈ [0, T], and A brief description on the mechanical interpretations of the equations and relations in the above problem is as follows. Formula (2.1) represents a general constitutive law for a viscoelastic material with long memory, where A is a viscosity operator, B is an elasticity operator, the long memory feature is reflected by the integral term in which R is the relaxation tensor (cf. [12]). The process is quasi-static and (2.2) is the equilibrium equation, where f 0 is the density of the volume forces applied to the deformable body. The homogeneous displacement boundary condition (2.3) is specified on 1 . A traction boundary condition (2.4) is specified on 2 , where f 2 denotes the density of the traction; note that 2 is allowed to be an empty set, and then the boundary condition (2.4) is suppressed. Contact conditions are specified on 3,1 and 3,2 . Over 3,1 , (2.5) describes the frictionless contact with unilateral constraints on the velocity (cf. [4]). Over 3,2 , a version of Coulomb's law of dry friction with the simplifying assumption that the normal stress on the contact boundary is known (cf. [16]); in (2.6), F is a nonnegative valued function, μ ≥ 0 is a friction coefficient, and μF represents a friction bound. Since the problem involves the time derivative of the displacement, an initial value condition (2.7) is supplemented, u 0 being a given initial displacement field.
Let us introduce some function spaces and sets. The function space for stress and strain fields at a fixed time is This is a Hilbert space with the inner product and the induced norm τ Q = (τ , τ ) 1/2 Q . The function space for displacement and velocity fields at a fixed time is Since meas( 1 ) > 0, by the Korn inequality (cf. [14, p. 79]), V becomes a Hilbert space with the inner product (u, v) V = (ε(u), ε(v)) Q , and the associated norm v V = ε(v) Q is equivalent to the standard H 1 ( ) d -norm on V . The set of admissible velocity fields is We further introduce the space of bounded, symmetric fourth-order tensor fields which is a real Banach space with the norm We make the following assumptions on the problem data.
for t ∈ [0, T]. By a standard procedure, one can obtain the following weak formulation of the contact problem (2.1)-(2.7). and The following well-posedness result is shown in [13]. The proof of Theorem 2.2 is carried out in [13] through an application of a result on an abstract sweeping process. In the rest of the section, we provide an elementary proof of the result through a fixed-point argument. For this purpose, we introduce the velocity variable (2.10) Note that given w(t) and the initial displacement u 0 , we can recover the displacement u(t) by Then Problem 2.1 can be reformulated in terms of w.
Note that the variational inequality (2.12) contains a doubly-history dependent term For this reason, (2.12) is called a doubly-history dependent variational inequality. The existence and uniqueness result of Problem 2.3 is the following.

Theorem 2.4 Under assumptions H(A), H(B), H 1 (R), H(f), and H(C), for any
Proof We rewrite (2.12) as Let z ∈ C([0, T]; K). Then at any time t ∈ [0, T], z(t) ∈ K and the variational inequality has a unique solution w(t) ∈ K (cf. [1, Theorem 11.3.1]). Now we show that the function w : [0, T] → K is continuous. We take t 1 , t 2 ∈ [0, T] in (2.14) respectively to get and We take v = w(t 2 ) in (2.15) and v = w(t 1 ) in (2.16) then add the inequalities to get Using assumption H(A), we get Together with (2.18), we apply the Cauchy-Schwarz inequality on (2.17) to derive By assumption H(B) and (2.11), we obtain Furthermore, where Using (2.11), we obtain  Let z 1 , z 2 ∈ C([0, T]; K) and w 1 , w 2 be the solution to (2.14) with z = z 1 and z = z 2 , respectively. Similar to the arguments in the derivation of (2.19), we get i.e.,

Numerical analysis of the variational inequality
In this section, we introduce a fully discrete method for Problem 2.3 and derive an optimal order error bound. Assume that is a polygonal/polyhedral domain, and let {T h } be a regular family of partitions of into triangles/tetrahedrons that is consistent with the splitting of the boundary into four subsets 1 , 2 , 3,1 , and 3,2 , i.e., if a side/face of an element has a nontrivial intersection with one of the four subsets, then the side/face lies entirely on the closure of that subset. Let be the linear element space for {T h } and K h be a subset of V h which is defined by In this section, C denotes a generic constant which takes on different values in different places. Denote w(t n ) = w n , f(t n ) = f n . Let u h 0 ∈ K h be an approximation of u 0 . Then the fully discrete scheme for Problem 2.3 is introduced as follows.

Problem 3.1 Find a discrete velocity field w kh
We have the following result for Problem 3.1.

Theorem 3.2 Under assumptions H(A), H(B), H 1 (R), H(f), and H(C)
, for any u 0 ∈ V , Problem 3.1 has a unique solution w kh ⊂ K h .
Proof For n = 0, (3.1) is to find w kh 0 ∈ K h such that From the assumptions on the data, we know this elliptic variational inequality has a unique solution w kh 0 ∈ K h [1, Theorem 13.3.1]. For an integer 1 ≤ n ≤ N , suppose that w kh j , 0 ≤ j ≤ n -1, are known. We rewrite (3.1) as Fv h ν da.
Again, from the assumptions on the data, this elliptic variational inequality has a unique solution w kh n ∈ K h . Thus, by an induction argument, we know that there is a unique solution w kh ⊂ K h to Problem 3.1.
To derive an optimal order error estimate for Problem 3.1, we make further assumptions: Now we provide an optimal order error estimate for the numerical method defined by Problem 3.1.

Theorem 3.3 Assume H(A), H(B), H 2 (R), H(f), H(C), and H(u 0 ). Moreover, assume solution regularities H(w) and H(σ ). Let u h
0 ∈ K h be the finite element interpolation or the L 2 ( )-projection of u 0 . The following error bound holds for the numerical solution of Problem 3.1: Proof We establish the framework for error estimation. Using v V = ε(v) Q and assumption H(A), we have Now we introduce a residual term R n (v h , w n ) which plays an important role in the proof. Taking t = t n and v = v h in (2.12), we get Then define R n (v h , w n ) by Note that v h ∈ K h and K h ⊂ K ; therefore R n (v h , w n ) ≥ 0. Taking v h = w kh n in (3.4), we have In addition, (3.1) can be rewritten as follows: By (3.3), (3.6), and (3.7), we derive an error bound of w n -w kh n V as follows: where In the following we bound each of the terms I i , 1 ≤ i ≤ 4. Firstly, we focus on the estimation of I 1 . By Taylor's formula, for s ∈ [t j , t j+1 ], there exists ξ j ∈ [t j , t j+1 ] such that where w (ξ j ) ∈ V is the derivative of w at t = ξ j . Thus the following error bound holds for the left endpoint rule: Using the Cauchy-Schwarz inequality and (3.9), we obtain where Then we bound the term I 2 . First, where By the norm triangle inequality, where We take the square root of both sides of (3.30) to get Recall discrete Gronwall's inequality [5]: Next, we derive an upper bound for the residual term I 4 . The weak formulation (2.12) is used to get point-wise equations. Define a function space (3.33) Note that v = 0 on 3,2 , we get Recall Green's formula Since v is arbitrary, we can follow the technique in [5,Sect. 8.1] to obtain the following point-wise relations: Div σ (t) + f 0 (t) = 0 a.e. in (3.37) and σ (t)ν = f 2 (t) a.e. on 2 . (3.38) Denote σ (t n ) = σ n , f 2 (t n ) = f 2n , and f 0 (t n ) = f 0n . We take t = t n in equation (3.37) and multiply (v h -w n ) to both sides of the equation. Integrating over , we get Applying Green's formula again, we have We use (3.40), (3.38), and (3.39) to derive Combining (3.5), (3.41), H(f), and H(σ ), we get Substituting (3.42) into (3.32), we obtain Finally, we draw the conclusion of the numerical analysis for the fully discrete method. Using the regularity assumptions u 0 ∈ H 2 ( ) d , w ∈ C([0, T]; H 2 ( ) d ) and finite element interpolation error estimates [2,3], we derive the result of Theorem 3.3 that max 0≤n≤N w n -w kh n V ≤ C(h + k), (3.44) i.e., the fully discrete scheme is of first order both in the time step-size and the spatial meshsize.

Numerical examples
In this section we report numerical examples for the contact problem. A threedimensional viscoelastic body is in contact with a foundation. The body is subjected to the action of traction. Let domain represent the cross section of the body and assume that the plane stress hypothesis is valid. As depicted in Fig. 1, = (0, L 1 ) × (0, L 2 ) is a rectangle whose boundary is divided as follows:  The viscosity tensor A satisfies (Aτ ) ij = 2θτ ij + ζ (τ 11 + τ 22 )δ ij , 1≤ i, j ≤ 2, (4.1) where θ and ζ represent the viscosity coefficients satisfying θ > 0 and ζ ≥ 0, and δ ij denotes the Kronecker symbol. The elasticity tensor B satisfies (Bτ ) ij = E 1 + κ τ ij + Eκ 1κ 2 (τ 11 + τ 22 )δ ij , 1≤ i, j ≤ 2, (4.2) where E denotes the Young modulus, and κ denotes the Poisson ratio of the material. The relaxation tensor R(s) = e -10s I, where I is an identity matrix. The fully discrete method is used to solve the contact problem. The intervals [0, L 1 ] and [0, L 2 ] are divided into 1/h equal parts, and uniform triangular finite element partitions are introduced; see   The configuration of the deformable body with h = k = 1 32 at t = 1s is illustrated in Fig. 3. The numerical solution with h = k = 1 128 is used as the "reference" solution, and the convergence orders in H 1 norm are reported in Table 1, Table 2, and Table 3 which illustrate the performance of the fully discrete method.