Abstract

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.

1. 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, [16]. Among these laws, the normal compliance condition introduced in [7] remains one of the most popular contact models used in the literature; see [4, 811]. 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 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 [1417] 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 [1416], 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.

2. 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 for the set of positive integers and will represent the set of nonnegative real numbers; that is, . Let . Then, we denote by the space of second-order symmetric tensors on . The inner product and norm on and are defined by Here and below the indices , , , and run between and and, unless otherwise stated, the summation convention over repeated indices is used.

Let be a bounded domain with a Lipschitz continuous boundary and let be a measurable part of such that . We use the notation 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 spacesHere and represent the deformation and divergence operators given by The spaces , , , and are real Hilbert spaces endowed with the inner products and the associated norms , , and , respectively. We denote by and the normal and the tangential component of on , respectively, given by , . For a regular function 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 , the directional derivative of at in the direction is defined by whenever this limit exists. The Clarke generalized directional derivative of a locally Lipschitz function at the point in the direction is defined by The Clarke subdifferential of at is a subset of given byA locally Lipschitz function is said to be regular (in the sense of Clarke) at if, for all , the directional derivative exists and . 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 be given. For a positive integer we define . Assume that and are two sequences of nonnegative numbers satisfying for a positive constant independent of or . Then there exists a positive constant , independent of or , such that

3. 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 , , and such that the (. We use the notations and for the displacement field and the stress field, respectively; therefore, , , 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 be the time interval of interest with . The body is clamped on and, therefore, the displacement field vanishes there. A volume force of density acts in , and surface tractions of density act on . In this study, we consider the dynamic contact with friction between a viscoelastic body and a foundation on . The material’s behavior is modelled by a viscoelastic constitutive law of the formHere is the linearized strain tensor, is a nonlinear operator which describes the viscous properties of the material, and 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 for and represents the positive normal damped response function such that for . 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 . Therefore, when , the friction condition is described byDetails 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 . Find a displacement field and a stress field such thatfor all 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 and 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 , we consider a closed subspace of

Completeness of the space follows from the assumption , which allows the use of Korn’s inequality, for some constant , depending only on and ,On , we use the inner product given byand let be the associated norm; that is,

It follows from (21) and (23) that 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 depending only on the domains , , and such thatBy (26) there exists a continuous trace operator and for the function we still denote by its trace . In what follows, we need the spaces , , and where the time derivative involved in the definition of is understood in the sense of vector valued distributions. Equipped with the norm , the space becomes a separable Hilbert space. We also have and the following inequality holds, where is a constant,It is well known that the embedding and embedding are continuous. Moreover, choosing we introduce the space equipped with norm. Note that and with both embeddings being compact, and hence where both embeddings are again compact.

We introduce now Ehrling’s lemma, proved in [21].

Lemma 2. Let , , and be three Banach spaces such that with the injection of into being continuous and the injection of into being compact. Then, for each , 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..The viscosity operator satisfies.The elasticity operator satisfies.The normal compliance function satisfies.The normal damped response function satisfies.The friction bound satisfiesThe 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 satisfyAlso, we assume that the initial values satisfy the following.. and .

Using the Clarke subdifferential (7), the friction condition (18) can be expressed in another form. Indeed, define a function byThen, assuming , the condition (18) is equivalent to the following subdifferential inclusion: where denotes the Clarke subdifferential of at the point .

Properties of the function are summarized in the next lemma.

Lemma 3. If the assumptions hold, then the function defined by (37) is locally Lipschitz, and If furthermore the assumption holds, then we have

Proof. First, we have to proof that is locally Lipschitz. In order to do so, let and . For , we get With 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 , and , defined by

If holds, the operator has the following properties:Under the assumptions , the operator satisfies the following properties:Assuming , the functional is locally Lipschitz, and we have the following inequalities on :

Now, multiplying the equation of motion by the test function , integrating over , using the Green formula and the definition of the operators, we obtain the following variational formulation of Problem .

Problem . Find a displacement field with and a friction stress field such that, for a.e. and for all ,with where is the set of all selections of and

Finally, we formulate the result concerning the existence and uniqueness of solution to Problem . In order to keep the paper with a reasonable length, we skip the proof of this result.

Theorem 4. Assume that , , , , , , and hold. Then, Problem has a unique solution.

Note that the proofs of this theorem are based on similar arguments to those used in [14, 17, 23].

4. Fully Discrete Error Estimates

We introduce a fully discrete approximation of Problem in order to bound the error of the fully discrete solutions. To this end, on a finite time interval with given, we consider a positive integer and we define the time step-size and the time nodal points and .

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 . Here denotes the spatial discretization parameter.

In what follows, we denote and . For a sequence , we use the notation , , for the backward divided differences, as well as the additional notation . Let be suitable approximations of , characterized byfor all . It is easy to observe that The fully discrete approximation of Problem is the following.

Problem . Find a velocity field and a friction stress field such thatwith We consider the following discrete displacement

Remark 5. From now on, the constant denoted by may differ from line to line.

In what follows, our aim is to analyze the error between the solution of Problem in displacement and its fully discrete approximation, solution of Problem . 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 , , , , , , and and the regularity and we have the optimal order error estimate with a constant independent of and .

Proof. We take in (47) and we combine it with (52) in order to getNow, taking in (56) and combining these 2 inequalities, we obtainAfter reformulation of (57) under the form , we haveFrom that, we deduce thatTherefore, after some elementary manipulations, it is easy to see thatRemark  7. From now on, the constant denoted by may differ from line to line.
It is easy to see that Then, we take and From the coercivity of the operator in (44), we getThe operator is Lipschitz continuous, from (44):The operator is also Lipschitz continuous, from (45):With (26), (33)(a), and (34)(a), we haveFrom (30), (33)(b), (34)(b), and (46)(c), we getUsing now (26), (33)(a), and (46)(a), we haveThen, with (26), (34)(a), and (46)(a), we haveTherefore, using this time only (26) and (33)(a) we deduce thatand from (26), (30), and (34)(a), we obtainThus, from (62)–(71), we getThen, from (25), (26), (30), the inequalities and , and (72), we obtainwith a positive constant that may change from line to line.
Remark 8. Note that if is small enough, .
Now, we replace by and sum over from to to obtainWe consider the term ; we havewith . We recall the following classical inequality:where is defined by the following relationThus, since ,Therefore,By considering the constant , it is possible to bound from below by . Note that is positive due to Remark 8.
DenoteThen, we can reformulate (74) as follows:From (81), we apply the discrete form of the Gronwall Lemma 1 to obtain thatNow, let be the finite element interpolations of . It follows from [24] that Then, the error bound follows from Theorem 6 and (50).

5. Numerical Simulations

The aim of this section is to present the numerical strategy used to solve the frictional contact Problem , to provide numerical simulations and also to get a numerical evidence of the convergence of the discrete scheme established in Section 4.

Numerical Solution. The numerical solution is based on an iterative procedure which leads to a sequence of convex programming problems already used in [14, 23]. For each “convexification” iteration of index , the value of the friction coefficient is fixed to a given value depending on the tangential velocity solution found in the previous iteration. Then, the resulting nonsmooth convex iterative problems are solved by classical numerical methods. Furthermore, the frictional contact conditions are treated by using a numerical approach based on the augmented Lagrangian method (cf. [25]). To this end, we consider additional fictitious nodes for the Lagrange multiplier in the initial mesh. The construction of these nodes depends on the contact element used for the geometrical discretization of the interface . In the case presented below, the discretization is based on “node-to-rigid” contact element, which is composed by one node of and one Lagrange multiplier node. For more details on the discretization step and Computational Contact Mechanics, we refer to [2, 6, 2527]. The numerical solution of the nonsmooth nonconvex variational Problem is based on the iterative scheme given in the following algorithmic lines characterized by the time stepping loop of index and the convexification loop of index . The convexification loop is ended by the stopping criterion given by (85).

Let , and let be given.

For , do time stepping loop.For , do convexification loop.Problem Find a velocity field and a friction stress field such that until , and We recall that the discrete displacement is given by .

Note that the different numerical methods used for the numerical algorithm (85) have been implemented in a Fortran computer code which is based on a MODULar Finite Element library (MODULEF) developed by INRIA (Institut National de Recherche en Informatique et en Automatique, Rocquencourt, France). For more details, we refer to https://www.rocq.inria.fr/modulef/english.html.

Numerical Example. We consider the physical setting depicted in Figure 1. There, with , and The domain represents the cross section of a three-dimensional linearly viscoelastic body subjected to the action of tractions in such a way that the plane stress hypothesis is assumed.

On the part , the body is clamped and, therefore, the displacement field vanishes there. Vertical tractions act on the part of the boundary. 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 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 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 and are very large and . 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 of the formwith 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 , 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 , we consider the following function for , with , the abscissa of a point on the boundary where .

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 at the moments , and . Note that at a given instant , significant normal contact forces only exist on a portion of the boundary and those forces will eventually vanish over time. Such a behavior is consistent with the movement of a mass along the boundary . Besides, the velocity of the mass along the boundary can be assimilated to the coefficient . This is why, with the choice , the mass proceeds from left to right through the whole boundary .

Next, we study the influence of the parameter on the normal contact stresses on . At a given instant , since the mass is located at , the forces decrease exponentially with the distance between the node considered on the contact boundary and the node at the abscissa . 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 numerical solution errors in the energy norm defined by which is equivalent to the norm . Since the true solution is not available, we use instead the numerical solution corresponding to fine discretization of as the “reference” solution in computing the solution errors. Here, the numerical solution with and is taken to be the “reference” solution . 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.

Competing Interests

The authors declare that they have no competing interests.