Numerical Analysis of a Sliding frictional contact problem with Normal Compliance and Unilateral Contact

. This paper represents a continuation of [15] and [18]. Here, we consider the numerical analysis of a non trivial frictional contact problen in a form of a system of evolution nonlinear partial diﬀerential equations. The model describes the equilibrium of a viscoelastic body in sliding contact with a moving foundation. The contact is modeled with a multivalued normal compliance condition with memory term restricted by a unilateral constraint, and is associated to a sliding version of Coulomb’s law of dry friction. After a description of the model and some assumptions, we derive a variational formulation of the problem, which consists of a system coupling a variational inequality for the displacement ﬁeld and a nonlinear equation for the stress ﬁeld. Then, we introduce a fully discrete scheme for the numerical approximation of the sliding contact problem. Under certain solution regularity assumptions, we derive an optimal order error estimate and we provide numerical validation of this result by considering some numerical simulations in the study of a two-dimensional problem.


Introduction
The modelling and the analysis of frictional contact problems represent important topics both in Engineering Sciences and Applied Mathematics, see for instance the numerous works in these fields of research [1,2,6,7,8,11,12,13,14,15,16,17,18,19,20] and the references therein. The mathematical study of contact and friction models concerns the variational and the numerical analysis with also significant contributions related to numerical simulations.
In the construction of the models, various contact conditions have been considered. We can cite for instance, the so-called Signorini condition, introduced in [14], the so-called normal compliance contact condition, introduced in [12] and used in a large number of papers, see [2,6,7,13,15,16,17,18]. Recently, a more general contact condition, called the normal compliance condition restricted by unilateral constraint introduced in [8], models the contact with an elastic-rigid foundation. The mathematical analysis of models involving the frictionless contact condition with normal compliance and unilateral constraint can be found in [7,8,16]. When friction is considered, the unique solvability of the variational problems can be proven by considering a smallness assumption of the friction coefficient, see for instance [2,17,18]. Recently, the study of a model of frictional contact problem between a linearly elastic body and a moving rigid foundation, considered in [1], has shown both the obtention of multiple solutions when the coefficient of friction is larger than a critical value, and the existence and uniqueness of the solution at low coefficient of friction.
In this work, we consider the frictional contact model introduced in [18] which describes the contact of deformable body with a moving foundation not perfectly rigid. Therefore, the contact law with normal compliance and unilateral constraint was associated to a sliding version of Coulomb's law of dry friction. Furthermore both the material constitutive law of the body and the frictional contact model are characterized by memory terms in order to take into account physcical relaxation behaviors. Such kind of historydependent problems have been considered in [17,18,19].
The present paper represents a continuation of several papers [2,17,18]. Here, our goal is to provide the numerical analysis of the frictional contact model with normal compliance and unilateral constraint introduced in [18], and to illustrate the error estimate of the discretisation by numerical simulations. The mathematical model is based on a viscoelastic constitutive law with long memory, a contact conditions combining normal compliance, memory term, unilateral constraint and a frictional sliding version of the Coulomb's law. This nonstandard mathematical model can be formulated by a history-dependent quasivariational inequality for the displacement field.
The rest of the paper is structured as follows. In Section 2 we present the mathematical model by describing its equations and boundary conditions. Some notation, assumptions and preliminary material have been introduced in order to derive a variational formulation of the problem. The Section 3 is devoted to the numerical approximation and the numerical analysis of the variational problem considered in the previous section. Under certain solution regularity assumptions, we derive an optimal order error estimate that represents the main result of this work. Finally, in Section 4 we consider some numerical simulations in the study of a two-dimensional problem and we provide a numerical validation of the error estimate.

Frictional contact problem and variational formulation
First, we present some preliminary material usefull for the setting of the problem. Let Ω a regular domain of R d (d = 2, 3) with its boundary Γ that is partitioned into three disjoint measurable parts Γ 1 , Γ 2 and Γ 3 , such that meas (Γ 1 ) > 0. We use the notation x = (x i ) for a typical point in Ω and we denote by ν = (ν i ) the outward unit normal defined almost everywhere (a.e.) on Γ. We denote by u = (u i ), σ = (σ ij ), and ε(u) = (ε ij (u)) the displacement vector, the stress tensor, and the linearized strain tensor, respectively. Here and below the indices i, j, k, l run between 1 and d and, unless stated otherwise, the summation convention over repeated indices is used. We note that sometimes we do not indicate the dependence of various vectors and tensors on the spatial variable x and we recall that the components of the linearized strain tensor ε(u) are given by ε ij (u) = 1 2 (u i,j + u j,i ) where the index that follows a coma indicates a partial derivative with the corresponding component of the spatial variable x, e.g. u i,j = ∂u i /∂x j . We denote by t the time variable and a dot superscript represents the time derivative with respect to the time variable t, e.g.u = ∂u/∂t.Furthermore, we use the notation N for the set of positive integers and R + will represent the set of nonnegative real numbers, e.g. R + = [0, +∞). Then, we denote by S d the space of second order symmetric tensors on R d . The inner product and norm on R d and S d are defined respectively by For all vector v, we denote by v ν and v τ the normal and tangential components of v on Γ given by We also recall that, if σ is a regular function, then the normal and tangential components of the stress field σ on Γ are defined by In this work, we consider an viscoelastic body that occupies the bounded domain Ω with Γ, its boundary. The body is clamped on Γ 1 and, therefore, the displacement field vanishes there. A volume force of density f 0 acts in Ω, surface tractions of density f 2 act on Γ 2 . On Γ 3 , the body is in frictional contact with a moving obstacle, the so-called foundation. We suppose that the foundation is plane and moves steadily, i.e. its velocity v * = v * τ * is assumed to be larger than the tangential velocityu τ on the surface contact Γ 3 (i.e. v * ≫ u τ ), where τ * = v * v * denotes a given unitary vector in the tangential plane and the value v * > 0 is also given. The contact between the body and the foundation is modeled by a multivalued normal compliance and a unilateral constraint. The associated friction is based on a version of Coulomb's law of dry friction, in which contact surfaces are assumed to be in relative slip status. The problem is studied in the interval of time R + . Thereby, let us consider the formulation of our quasi-static frictional contact problem defined as follows. Problem P. Find a displacement field u : Ω × R + → R d and a stress field σ : Ω × R + → S d such that, for all t ∈ R + ,

7)
and there exists a normal reaction π : Now, we shortly describe the physical meaning of relations (2.3)-(2.8). Equation (2.3) represents the viscoelastic constitutive law with long memory in which E is the elasticity operator and M is a relaxation tensor. Details and mechanical interpretation concerning such kind of laws can be found in [5], for instance. Equation (2.4) represents the equation of equilibrium in which Div denotes the divergence operator for tensor valued functions (Div σ = (σ ij,j )). Conditions (2.5) and (2.6) are the displacement boundary condition and the traction boundary condition, respectively. Finally, (2.7) and (2.8) represents the friction Coulomb's law and the multivalued normal compliance contact condition with unilateral constraint and memory term, respectively. The friction condition (2.7) represents a regularized form of a version of Coulomb's law in slip status where µ represents the coefficient of friction and R is a regularization operator. Here, σ ν and σ τ represent the normal contact stress and the tangential friction stress on the contact surface Γ 3 , respectively. The conditions (2.8) represents a version of the contact boundary conditions with normal compliance and unilateral constraint, in which the memory effects of the foundation are taken into account. Note that the condition (2.8) models the contact with a foundation made of a rigid material and covered by a layer of soft material (asperities) of thickness g with a thin crust. Let us describe the different terms. g is the maximal penetration allowed in the foundation and represents the size of the soft material. −σ ν (t) = p(u ν (t)) + π(t) is the normal compliance condition with memory term, where the normal compliance function p is a Lipschitz continuous increasing function which vanishes for a negative argument and π describes the memory properties of the foundation. Here, F and N are positive functions related to the history of the soft material.
We now turn to the variational formulation of Problem P and, to this end, we consider standard notation for the Lebesgue and the Sobolev spaces associated to Ω and Γ. Also, we introduce the spaces The spaces Q and Q 1 are Hilbert spaces with the canonical inner product given by and the associated norms · Q and · Q1 . Since meas(Γ 1 ) > 0, V is a real Hilbert space with the inner product and the associated norm · V . By using the Sobolev trace theorem, a positive constant c 0 exists such that In addition, for σ ∈ Q 1 we denote by σ ν ∈ H − 1 2 its normal component, in the sense of traces. Let R : H − 1 2 (Γ) → L 2 (Γ) be a linear continuous operator. Then, there exists a positive constant c R > 0 such that Since σ is regular, the following Green's formula holds: Let us denote by Q ∞ the space of fourth order tensor fields given by Let us recall that Q ∞ is a real Banach space with the norm E Q∞ = max 0≤i,j,k,l≤d E ijkl L ∞ (Ω) and, the following inequality holds, Based on definitions of the norms · Q and · Q∞ , a proof of the inequality (2.12) is obtained by using simple calculations. Furthermore, we need to consider the sets of admissible displacement and admissible constraints defined by (2.14) Let us consider the assumptions on the operators, tensors and data. To this end we assume that the elasticity operator E and the relaxation tensor M satisfy the following conditions. (2.16) The densities of body forces and surface traction have the following regular- The normal compliance function p and the surface yield function F satisfy for any r ∈ R.  Finally, the surface memory function N and the coefficient of friction µ verify In what follows we assume that (u, σ) are sufficiently regular functions and satisfy (2.3)-(2.8). Here, let v ∈ U and t > 0 be given. First, by using the Green's formula (2.11) and the equilibrium equation (2.4), we obtain that For more details on getting inequality (2.22), we refer to [18], pp 2874-2875.
To simplify the notation, we defined the operators E : (2.29) Using the previous notation and substituting (2.3) in (2.22), we obtain the following variational formulation of Problem P.
Problem P V . Find a displacement field u : R + → U and a stress field σ : R + → Σ, for all t ∈ R + , such that, Under the assumptions (2.15)-(2.21) and an additional smallness assumption on the friction coefficient µ, a result of existence and uniqueness for the variational problem P V was provided in [18]. Based on this previous variational formulation, our goal in the next section is to provide the numerical analysis of this specific and non trivial contact problem.

Variational Approximation and error analysis
This section is devoted to the numerical discrete approximation of the Problem P V . In particular, an optimal error estimate is provided and represents the main result of the paper. More precisely, we are interested in solving the Problem P V over a finite time interval [0, T ], with T > 0 arbitrary but fixed. Thus, let N be a positive integer; we define the size of the time step k = T /N and we consider the uniform temporal discretization characterized by the time instants t n = nk for 0 ≤ n ≤ N . We use the notations f n = f (t n ), u n = u(t n ). For simplicity, we assume that Ω is a polygonal domain. Let consider {T h } a regular family of partitions of Ω into triangles that are compatible with the partition of the boundary Γ into Γ 1 , Γ 2 , and Γ 3 , in the sense that if the intersection of one side of an element with one of the three sets has a positive surface measure, then the side lies entirely in that set. Then we consider linear element spaces corresponding to T h , where P 1 (T ) stands for the space of polynomials of a degree less than or equal to 1 on T . We define U h ⊂ V h a non-empty, convex, and closed set, approximating U by where h > 0 denotes the spatial discretization parameter. For the discretization of the integral term of Problem P V , we use a variant of the trapezoid method where a prime indicates that the first and last terms in the summation are to be halved. Then, we introduce the approximations of operators S, M and N ,

4)
From now on, c will represent various positive constants which do not depend on h and k and whose value may change from line to line. Then the fully discrete scheme of Problem P V is the following.
Problem P hk V . Find u hk := u hk n N n=0 ∈ U h such that n = 0, ...., N , Note that the existence and uniqueness of the solution of the Problem P hk V can be provided by using the same arguments as for the solvability of the variational problem P V (see [15] and [18]). We now focus on the error analysis between the solutions to Problems P V and P hk V . Taking t = t n and v = u hk n in (2.31), we have Then, we consider the term Furthermore, by using (2.15(c)), we obtain We combine the inequalities (3.6), (3.7) and (3.8) to obtain

10)
K 1 = j(S hk n u hk , u hk n , u n ) − j(S hk n u hk , u hk n , u hk n ) +j((Su) n , u n , u hk n ) − j((Su) n , u n , u n ), (3.11) (3.12) Now, we proceed to the estimation of each terms of (3.9). First, we use the definition (2.23) and the inequalities (2.10), (2.15(b)) and (2.18(c)) to obtain For the term K 1 , we use arguments similar to them used in [18, p. 2877], to obtain where α and β are constants, for more details about this constants, we refer to [18].
Moreover, we use the triangular inequality as well as the trapezoid method introduced in (3.2) to obtain Then we combine the inequalities (3.14), (3.15), (3.18), (3.19) and the trace inequality (2.10) to obtain By a similar reasoning applied for the term K 1 , we obtain that Finally, we gather together the inequalities (3.9), (3.13), (3.20) and (3.21) to deduce that Under the assumption of smallness for the uniqueness of the Problem P V (see, [18, p. 2874]) and by using the elementary inequality ab ≤ δa 2 Now let us consider the following Gronwall inequalities.
To proceed further, we make the following solution regularity assuptions: Then it can be shown that the relations (2.3)-(2.8) are satisfied pointwise a.e. and we have Let Π h u ∈ V h be the linear finite element interpolant of the solution u. As the solution u ∈ U , i.e., u ν ≤ g, then Π h u ∈ U h . The standard finite element interpolation theory yields (cf. [4]) and we have In conclusion, we have shown the following result.

Numerical simulations
This section provides computer simulation results on the contact Problem P hk V , including numerical evidence of the theoretical error estimates obtained in Section 3 for the discrete approximation of the variational problem P V . The solution of Problem P hk V is based on numerical methods described in [2,15]. For a large overview on numerical methods to solve contact problems, we can refer for instance to [10,11,20].
Numerical example. The physical setting used for Problem P hk V is depicted in Figure 1. Here, we consider the frictional contact between a deformable body and a moving foundation. This specific foundation is composed by a rigid material covered by a thin crust and a deformable layer of asperities of depth g.
Here g represents the maximum value of the allowed penetration in the foundation. When this value of penetration is reached, the contact follows a unilateral condition without any additional penetration. This kind of foundation is characterized by the contact conditions (2.8). Since the foundation is moving the friction condition is in a slip status within the Coulomb's form (2.7). The deformable body is a rectangle, Ω = (0, 2) × (0, 1) ⊂ R 2 , and its boundary Γ is split as follows: The domain Ω represents the cross section of a three-dimensional linearly viscoelastic body subjected to the action of tractions in such a way that a plane stress hypothesis is assumed. On the part Γ 1 the body is clamped and, therefore, the displacement field vanishes there. Vertical compressions act on the part (0, 2) × {1} of the boundary Γ 2 and the part ({2} × [0, 1)) is traction free. Constant vertical body forces are assumed to act on the viscoelastic body. The body is in frictional contact with an obstacle on the part Γ 3 of the boundary. The compressible material's behaviour of the domain Ω is governed by a viscoelastic constitutive law of the form (2.3). In addition, we assume that the material is homogeneous and isotropic; then, the elasticity Ω deformable body where the coefficients E and κ are the Young's modulus and the Poisson's ratio of the material, respectively, and α is a viscosity parameter. δ ij denotes the Kronecker symbol. For the numerical simulation of Problem P hk V , the data concerning the material are the following: For the contact boundary conditions on Γ 3 , we recall that the friction is in a sliding status and the normal contact response follows a multivalued normal compliance condition with respect to the normal displacement u ν and for which the maximal penetration is restricted by a unilateral constraint. Then the data concerning this specific frictional contact model (2.7) and (2.8) are the following: p(r) = c ν max{0, r}, c ν = 100 N/m 2 , F (s) = s, N (s) = βs, β = 10 and g = −0.04 m.
Numerical solution of Problem P hk V . In Figure 2, we present the deformed configuration as well as the interface forces on Γ 3 . On the right extremity of the boundary Γ 3 , we can see that a large number of contact nodes of Γ 3 are in several contact status with the foundation. We note that some of these nodes have first broken the crust (crust contact), other have crushed the asperities Errors and numerical convergence orders. The aim of this part is to illustrate the convergence of the discrete scheme and to provide a numerical evidence of the optimal error estimate obtained in Section 3. To this end, we computed a sequence of numerical solutions by using uniform discretization of Problem P hk V according to the spatial discretization parameter h and the time step k, respectively. For instance, the deformed configuration and the interface forces plotted in Figure 2 are obtained for for h = 1/128 and k = 1/128 and corresponds to a problem with 17028 degrees of freedom.
The numerical estimations of u − u hk V computed for several discretization parameters of h and k, have been estimated by using the energy norm · E defined by v hk E := 1 √ 2 (E(ε(v hk )), ε(v hk )) 1/2 H . Since the exact solution u can not be calculated analytically, we consider as "reference" solution u ref a numerical solution which corresponds to a fine approximation of Problem P hk V . For this procedure, the boundary Γ of Ω is divided into 1/h equal parts and the time interval [0, T ] is divided into 1/k time steps. We start with h = 1/2 and k = 1/2 which are successively halved. The numerical solution u ref corresponding to h = 1/256 and k = 1/256 was taken as the "reference" solution. This fine discretization corresponds to a problem with 66820 degrees of freedom; the simulation runs in 129521 (expressed in seconds) CPU time on a IBM computer equipped with Intel Dual core processors (Model 5148, 2.33 GHz). The numerical results presented in Figure 3 provide a good numerical evidence of the theoretically predicted order convergence of the numerical solution measured in the energy norm.