Analysis of a Dynamic Viscoelastic Contact Problem with Normal Compliance , Normal Damped Response , and Nonmonotone Slip Rate Dependent Friction

We consider a mathematical model which describes the dynamic evolution of a viscoelastic body in frictional contact with an obstacle.The contact is modelled with a combination of a normal compliance and a normal damped response law associated with a slip rate-dependent version of Coulomb’s law of dry friction. We derive a variational formulation and an existence and uniqueness result of the weak solution of the problem is presented. Next, we introduce a fully discrete approximation of the variational problem based on a finite element method and on an implicit time integration scheme. We study this fully discrete approximation schemes and bound the errors of the approximate solutions. Under regularity assumptions imposed on the exact solution, optimal order error estimates are derived for the fully discrete solution. Finally, after recalling the solution of the frictional contact problem, some numerical simulations are provided in order to illustrate both the behavior of the solution related to the frictional contact conditions and the theoretical error estimate result.


Introduction
So far, dynamic contact problems abound in industry and everyday life.For this very reason, special care has been given to the modelling, mathematical analysis, and numerical solution of such problems for several decades, be it in the engineering or the mathematical literature.Such problems involve frictional contact phenomena that leads, because of their inherent complexity, to nonlinear, nonsmooth, and nonconvex mathematical problems.
Over the last decades, various contact boundary conditions have been used to model contact phenomena and their modelling is still under investigation; see, for instance, [1][2][3][4][5][6].Among these laws, the normal compliance condition introduced in [7] remains one of the most popular contact models used in the literature; see [4,[8][9][10][11].It represents a regularization of the so-called Signorini contact condition, expressed in terms of unilateral constraints for the displacement field, and described the contact with a deformable foundation.Such a law allows penetrations, while the Signorini condition does not.Note that these contact conditions, used in most of the related works on the subject, are formulated in terms of the displacement field.However, normal compliance conditions expressed in terms of the normal velocity seem to be more appropriate when the contact surfaces are lubricated, as mentioned in [5,12].Such kind of conditions is called normal damped response conditions.In this current paper, we model the contact with new and nonstandard conditions which involve both the normal compliance law and normal damped response law and, as such, describes the contact with a specific particular deformable foundation.Furthermore, one cannot hope to fully grasp the phenomenon of contact without taking into account the friction.In the literature, it is generally modelled by either the so-called Tresca or Coulomb friction laws.Normally, one would expect the friction coefficient to be constant; however, such a classical approach has shown its own limits for friction-induced phenomena such as stick-slip motion.For this very reason, several authors have introduced nonmonotone versions of friction laws [5,13].
Due to all the motivations mentioned above, in this work, we consider the mathematical study of a dynamic 2 Advances in Mathematical Physics frictional contact problem in which the contact is modelled with a combination of normal compliance law and normal damped response law, associated with a nonmonotone friction law involving a slip rate dependent friction.
The current study represents continuation of [14][15][16][17] and its aim is to provide the variational and numerical analysis of a new and nonstandard dynamic frictional contact problem supported by numerical simulation.With respect to these articles, this paper presents several traits of novelty that we describe in what follows.Here, we no longer have finite penetration in contrast to [14][15][16], where a combination of a normal compliance law and unilateral constraint was used, nor finite velocity as in [17], where a combination of a normal damped response law and unilateral constraint in velocity was considered.Note that we consider here a dynamic contact model, since the static case had been studied in [15,16] for a material whose behavior was described by only an elastic constitutive law.We insist on the fact that, in the problem considered in this paper, the behavior of the foundation is described by a Kelvin-Voigt-like foundation modelled by a combination of normal compliance and normal damped response.Furthermore, such considerations lead to nonstandard Coulomb's law of dry friction where the threshold depends on the tangential velocity, through the coefficient of friction, and both the normal displacement and normal velocity, because of the normal compliance and normal damped response, respectively.Another trait of novelty arises from the mathematical analysis with in particular the numerical analysis of such a problem by considering a dynamic process.Indeed, in this study, we derive optimal order error estimates in the fully discrete case under regularity assumptions imposed on the exact solution.Finally, we provide numerical simulations to validate the bound of the error estimate and to illustrate the mechanical behavior of the concrete physical body.
The rest of the paper is structured as follows.In Section 2 we introduce the notation we shall use as well as some preliminary material.In Section 3 we present the classical formulation of the frictional contact problem, we list the assumptions on the data, and we derive the variational formulation of the initial and approximate problem.Then, the existence and uniqueness result of the weak solution of the problem is presented.Section 4 is devoted to the presentation of the error estimate result in the fully discrete numerical case.In Section 5, the numerical strategy used to solve the frictional contact problem is presented, followed by numerical simulations on a two-dimensional example including numerical validation of the optimal error estimate established in Section 4.

Notation and Preliminaries
We present the notation and some preliminary material which will be of use later on.Everywhere in this paper, we use the notation N for the set of positive integers and R + will represent the set of nonnegative real numbers; that is, R + = [0, +∞).Let  ∈ N.Then, we denote by S  the space of second-order symmetric tensors on R  .The inner product and norm on R  and S  are defined by Here and below the indices , , , and  run between 1 and  and, unless otherwise stated, the summation convention over repeated indices is used.
Let Ω be a bounded domain Ω ⊂ R  ( = 1, 2, 3) with a Lipschitz continuous boundary Γ and let Γ 1 be a measurable part of Γ such that meas(Γ 1 ) > 0. We use the notation x = (  ) for a typical point in Ω ∪ Γ and we denote by ^= (]  ) the outward unit normal at Γ. Also, an index that follows a comma represents the partial derivative with respect to the corresponding component of the spatial variable; for example,  , =   /  .We use standard notations for the Lebesgue and Sobolev spaces associated with Ω and Γ and, moreover, we consider the spaces Here  :  1 →  and Div :  1 →  represent the deformation and divergence operators given by The spaces , ,  1 , and  1 are real Hilbert spaces endowed with the inner products and the associated norms ‖ ⋅ ‖  , ‖ ⋅ ‖  , ‖ ⋅ ‖  1 and ‖ ⋅ ‖  1 , respectively.We denote by V ] and v  the normal and the tangential component of v on Γ, respectively, given by For a regular function  : Ω ∪ Γ → S  we denote by  ] and   the normal and the tangential components of the vector ^on Γ, respectively, and we recall that  ] = (^) ⋅ ^and   = ^−  ] ^.Next recall the definitions of classical (one-sided) directional derivative and its generalization in the sense of Clarke.Let  be a Banach space and  * its dual.For a function  :  → R  , the directional derivative of  at  ∈  in the direction V ∈  is defined by whenever this limit exists.The Clarke generalized directional derivative of a locally Lipschitz function  :  → R  at the point  ∈  in the direction V ∈  is defined by The Clarke subdifferential of  at  is a subset of  * given by A locally Lipschitz function  :  → R  is said to be regular (in the sense of Clarke) at  ∈  if, for all V ∈ , the directional derivative   (; V) exists and   (; V) =  0 (; V).The function  is regular (in the sense of Clarke) on  if it is regular at every point  ∈ .Now let us consider some material useful for the numerical analysis of the problem.We will need the following Gronwall inequalities, proved in [18].Lemma 1.Let  > 0 be given.For a positive integer  we define  = /.Assume that {  }  =1 and {  }  =1 are two sequences of nonnegative numbers satisfying for a positive constant c independent of  or .Then there exists a positive constant , independent of  or , such that

Mechanical Problem and Variational Formulation
We consider a viscoelastic body that occupies the bounded domain Ω with Γ, its boundary.We assume that Γ is a Lipschitz continuous boundary divided into three measurable parts Γ 1 , Γ 2 , and Γ 3 such that the meas(Γ 1 ) > 0. We use the notations u and  for the displacement field and the stress field, respectively; therefore, (u), u , and ü represent the linearized strain tensor, the velocity field, and the acceleration field, respectively.We denote by ] the unit outward normal, defined almost everywhere on Γ.Let [0, ] be the time interval of interest with  > 0. The body is clamped on Γ 1 × (0, ) and, therefore, the displacement field vanishes there.A volume force of density f 0 acts in Ω × (0, ), and surface tractions of density f 2 act on Γ 2 × (0, ).In this study, we consider the dynamic contact with friction between a viscoelastic body and a foundation on Γ 3 × (0, ).The material's behavior is modelled by a viscoelastic constitutive law of the form Here (u) is the linearized strain tensor, C is a nonlinear operator which describes the viscous properties of the material, and E is a nonlinear operator which describes its elastic behavior.Various examples and mechanical interpretations in the study of viscoelastic materials of the form (10) can be found in [19] and the references therein.Such kinds of constitutive laws were used in the literature in order to model the behavior of real materials like rubbers, rocks, metals, pastes, and polymers.In addition, the foundation is made of a soft material, where the contact is modelled with normal compliance and normal damped response.The friction is based on nonmonotone Coulomb's law of dry friction during the contact with the soft material.The boundary conditions on the contact surface are described in what follows.
When the body moves towards the obstacle, the contact is described by both a normal compliance condition and a normal damped response condition.In this case, the foundation is characterized by a Kelvin-Voigt like model.Therefore, Here  represents the positive normal compliance function such that () = 0 for  ≤ 0 and  represents the positive normal damped response function such that () = 0 for  ≤ 0. A similar rheological model was used in [20] to describe the behavior of soft clay with application on rheological consolidation.Furthermore, the contact is associated with Coulomb's law of dry friction in which the coefficient of friction  is assumed to depend on the slip rate ‖ u  ‖.Therefore, when  ] < 0, the friction condition is described by Details on the normal compliance conditions associated with Coulomb's law of dry friction can be found in [4,5,18], for instance.To simplify the notation and when it is sensible, we omit the dependence in position of the variables.With these preliminaries, the classical formulation of our dynamic frictional contact problem is the following.
Problem P. Find a displacement field u : Ω × (0, ) → R  and a stress field  : Ω × (0, ) → S  such that Advances in Mathematical Physics Div for all  ∈ R + and, moreover, Here, ( 13) represents the viscoelastic constitutive law of the material.Equation ( 14) is the classic equilibrium equation which includes the acceleration term;  denotes the density of the material and is assumed to be a constant for the sake of simplicity.Conditions (15) and ( 16) are the homogeneous displacement and traction boundary conditions, respectively, and (19) represents the initial conditions, with u 0 and u 1 being the initial displacement and velocity, respectively.Our interest lies in the contact conditions (17) and (18).In this case, the contact condition (17) allows nonlimited penetration.It could be assimilated to a Kelvin-Voigt law since the normal stress on the contact boundary is described by a combination of a normal compliance law and normal damped response law.The friction condition, described in (18), represents Coulomb's law of dry friction where  is the slip rate dependent coefficient of friction.
In order to introduce a weak formulation of the mechanical Problem P, we consider a closed subspace of  1 Completeness of the space (, ‖ ⋅ ‖  ) follows from the assumption meas(Γ 1 ) > 0, which allows the use of Korn's inequality, for some constant   , depending only on Ω and Γ 1 , On , we use the inner product given by and let ‖ ⋅ ‖  be the associated norm; that is, It follows from ( 21) and ( 23) that ‖ ⋅ ‖  1 and ‖ ⋅ ‖  are equivalent norms on  and therefore (, ‖⋅‖  ) is a real Hilbert space.Also, note that since we obtain that The duality pairing between  and  * is denoted by ⟨⋅, ⋅⟩.Identifying  with its dual, we have an evolution triple  ⊂  ⊂  * with dense, continuous, and compact embedding.By the Sobolev trace theorem and by (25), there exists a constant  0 depending only on the domains Ω, Γ 1 , and Γ 3 such that By (26) there exists a continuous trace operator  :  →  2 (Γ 3 ; R  ) and for the function v ∈  we still denote by v its trace v.In what follows, we need the spaces V =  2 (0, ; ), H =  2 (0, ; ), and W = {v ∈ V | v ∈ V * } where the time derivative involved in the definition of W is understood in the sense of vector valued distributions.Equipped with the norm , the space W becomes a separable Hilbert space.We also have W ⊂ V ⊂ H ⊂ V * and the following inequality holds, where   is a constant, It is well known that the embedding W ⊂ ([0, ]; ) and embedding {w ∈ V | ẇ ∈ W} ⊂ ([0, ]; ) are continuous.Moreover, choosing  ∈ (0, 1/2) we introduce the space  =  ∩  1− (Ω; R  ) equipped with  1− norm.Note that  ⊂  and  ⊂  with both embeddings being compact, and hence  ⊂  * ⊂  * where both embeddings are again compact.
Lemma 2. Let  0 , , and  1 be three Banach spaces such that with the injection of  into  1 being continuous and the injection of  0 into  being compact.Then, for each  > 0, there exists some constant () such that Using this lemma, we get the following inequality: In the study of the contact problem, we assume the following properties on the data.

𝐻(C). The viscosity operator
(E).The elasticity operator E : Ω×S  → S  satisfies (33) ().The normal damped response function  : The assumption (35)(c) represents a one-side Lipschitz condition, which allows the function (⋅) to decrease with a rate not faster than .().The force and the traction densities satisfy (36) Also, we assume that the initial values satisfy the following. 0 .u 0 ∈  and u 1 ∈ .
Using the Clarke subdifferential (7), the friction condition (18) can be expressed in another form.Indeed, define a function  : R  → R by Then, assuming ()(a)-(c), the condition ( 18) is equivalent to the following subdifferential inclusion: where () denotes the Clarke subdifferential of  at the point  ∈ R  .Properties of the function  are summarized in the next lemma.
Lemma 3. If the assumptions ()(a)-(b) hold, then the function  defined by ( 37) is locally Lipschitz, and If furthermore the assumption ()(c) holds, then we have Proof.First, we have to proof that  is locally Lipschitz.In order to do so, let  ∈ R  and  > 0. For  1 ,  2 ∈ (, ), we get With ()(a) and Proposition 5.6.28(ii) in [22], we have the following characterization of the Clarke subdifferential : Then, the other properties follow straightforwardly.
We define the operators  :  →  * ,  :  →  * , functional  :  2 (Γ 3 ; R  ) → R, and f : (0, ) →  * , defined by Under the assumptions (E), the operator  satisfies the following properties: Assuming (), the functional  :  2 (Γ 3 ; R  ) → R is locally Lipschitz, and we have the following inequalities on  ∈ (v): Advances in Mathematical Physics 7 Now, multiplying the equation of motion by the test function v, integrating over Ω × (0, ), using the Green formula and the definition of the operators, we obtain the following variational formulation of Problem P.
Problem P  .Find a displacement field u ∈ V with u ∈ W and a friction stress field  ∈  2 (0, ,  2 (Γ 3 ; R  )) such that, for a.e. ∈ (0, ) and for all v ∈ , with where  2 ( u  ()) is the set of all  2 selections of ( u  ()) and Finally, we formulate the result concerning the existence and uniqueness of solution to Problem P  .In order to keep the paper with a reasonable length, we skip the proof of this result.
Note that the proofs of this theorem are based on similar arguments to those used in [14,17,23].

Fully Discrete Error Estimates
We introduce a fully discrete approximation of Problem P  in order to bound the error of the fully discrete solutions.To this end, on a finite time interval [0, ] with  > 0 given, we consider a positive integer  and we define the time stepsize  = / and the time nodal points   =  and 0 ≤  ≤ .
We assume that Ω is a polyhedral domain and we consider a regular family of finite element partition of Ω.Let  ℎ ⊂  be the associate finite element space of continuous piecewise linear functions which vanish on Γ 1 .Here ℎ > 0 denotes the spatial discretization parameter.
In what follows, our aim is to analyze the error between the solution of Problem P  in displacement and its fully discrete approximation, solution of Problem P ℎ  .Note that the analysis of the error estimates in terms of tangential constraint remains an open problem.The main result of this section lies in the following theorem.Theorem 6. Assume  is sufficiently small.Then, under the assumptions (C), (E), (), (), (), (), and  0 and the regularity u ∈  1 (0, ; [ 2 (Ω)]  ) ∩  3 (0, ; ) and u  ∈ (0, ; with  a constant independent of ℎ and . Proof.We take v = v ℎ in (47) and we combine it with (52) in order to get Now, taking k ℎ = w ℎ  in (56) and combining these 2 inequalities, we obtain After reformulation of (57) under the form (⋅, w  − w ℎ  ) = (⋅, w  − k ℎ ), we have From that, we deduce that Therefore, after some elementary manipulations, it is easy to see that Remark 7. From now on, the constant denoted by  may differ from line to line.
It is easy to see that ( From the coercivity of the operator  in (44), we get The operator  is Lipschitz continuous, from (44): The operator  is also Lipschitz continuous, from (45): With ( 26), (33)(a), and (34)(a), we have From ( 30), (33)(b), (34)(b), and (46)(c), we get Using now ( 26), (33)(a), and (46)(a), we have Then, with ( 26), (34)(a), and (46)(a), we have Therefore, using this time only ( 26) and (33)(a) we deduce that Then, from ( 25), ( 26), (30), the inequalities  ≤ 2  No body forces are assumed to act on the viscoelastic body during the dynamic process.The body is in frictional contact with an obstacle on the part Γ 3 = [0, ] × {0} of the boundary.Note that we used the viscosity in (13) for mathematical reasons, so far, but from the practical point of view one may take the viscosity as small as one wishes.Here, since it is not our main interest, it is reasonable to neglect it.Indeed, this choice was motivated by the fact that our aim was to focus only on the behavior of the very specific frictional contact conditions.Therefore, the material response is governed by an elastic linear constitutive law defined by the elasticity tensor E given by Here,  and  are Young's modulus and Poisson's ratio of the material and   denotes the Kronecker delta.Note that, in order to insure the boundedness condition for both functions expressed in (33)(b) and (34)(b), one may consider the functions where  3 and  4 are very large and  + = max(0, ).Therefore, the normal compliance and normal damped response used here are defined by the following functions For the coefficient of friction, we choose a function  : with , ,  > 0 and  ≥ .Such a slip weakening phenomenon appears in the study of geophysical problems; see [28] for details.Indeed, in this case the coefficient of friction decreases with the slip from the value  to the limit value .And, for this reason, the corresponding friction law can be characterized as being nonmonotone.Since the function (⋅) is a contraction on [0, ∞), we have ( − ) < , and, as a consequence, the condition (35)(c) is verified.
Also, in order to describe the motion of a mass along the boundary Γ 2 , we consider the following function for f 2 , with , the abscissa of a point on the boundary where , ,  > 0.
The choice of this particular function is motivated by several points which will be highlighted through the numerical simulations presented below.For the computation, we use the following data:   In Figure 2, we plot the deformed configuration as well as the interface forces on Γ 3 at the moments  = 0.2 s, 0.4 s, 0.6 s, and 0.8 s.Note that at a given instant  =  0 , significant normal contact forces only exist on a portion of the boundary Γ 3 and those forces will eventually vanish over time.Such a behavior is consistent with the movement of a mass along the boundary Γ 2 .Besides, the velocity of the mass along the boundary can be assimilated to the coefficient  masse .This is why, with the choice c masse = /, the mass proceeds from left to right through the whole boundary Γ 2 .
Next, we study the influence of the parameter  on the normal contact stresses on Γ 3 .At a given instant  =  0 , since the mass is located at  0 =  masse  0 , the forces decrease exponentially with the distance between the node considered on the contact boundary and the node at the abscissa  0 .Here,  describes the decay rate of .Regarding the mass, it could be understood as its length: the longer the mass is, the lower  should be.Such an assumption is confirmed by Figure 3, where the deformed configurations and the interface forces are plotted for 2 values of .
This kind of simulation is supposed to model coarsely the action of a train wheel (mass) on a rail (deformable body) placed in contact on a ballast (deformable obstacle).
Error Estimate.In order to check the convergence of the discrete scheme and to illustrate the optimal error estimate obtained in Section 4, we report in Figure 4  solution u ref in computing the solution errors.Here, the numerical solution with ℎ = 1/128 and  = 1/128 is taken to be the "reference" solution u ref .This fine discretization corresponds to a problem with 33538 degrees of freedom and 32768 elements and was computed in 5001 CPU time (expressed in seconds) on an IBM computer equipped with Intel Dual core processors (Model 5148, 2.33 GHz).The curve of the numerical error estimate is asymptotically linear, which is consistent with the theoretically predicted optimal linear convergence of the numerical solution established in Section 4.

Figure 1 :
Figure 1: Reference configuration of the two-dimensional body.

Figure 2 :
Figure 2: Time evolution of the deformed meshes and contact forces on Γ 3 .