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

Authors: Mikaël Barboteu; David Danan
Journal: Advances in Mathematical Physics
 (2016)

---

## 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.


---

## Body

## 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, [1–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–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 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–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–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 Section2 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 notationN for the set of positive integers and R+ will represent the set of nonnegative real numbers; that is, R+=[0,+∞). Let d∈N. Then, we denote by Sd the space of second-order symmetric tensors on Rd. The inner product and norm on Rd and Sd are defined by (1)u·v=uivi,v=v·v1/2,∀u,v∈Rd,σ·τ=σijτij,τ=τ·τ1/2,∀σ,τ∈Sd. Here and below the indices i, j, k, and l run between 1 and d and, unless otherwise stated, the summation convention over repeated indices is used.LetΩ be a bounded domain Ω⊂Rd(d=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=(xi) for a typical point in Ω∪Γ and we denote by ν=(νi) 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, ui,j=∂ui/∂xj. We use standard notations for the Lebesgue and Sobolev spaces associated with Ω and Γ and, moreover, we consider the spaces(2)H=L2Ωd=u=ui∣ui∈L2Ω,Q=τ=τij∈L2Ωd×d∣τij=τji,H1=u=ui∣εu∈Q,Q1=τ∈Q∣Div⁡τ∈H.Here ε:H1→Q and Div:Q1→H represent the deformation and divergence operators given by(3)εv=εijv,εijv=12vi,j+vj,i,Div⁡σ=σij,j. The spaces H, Q, H1, and Q1 are real Hilbert spaces endowed with the inner products(4)u,vH=∫Ωuividx,σ,τQ=∫Ωσijτijdx,u,vH1=u,vH+εu,εvQ,σ,τQ1=σ,τQ+Div⁡σ,Div⁡τH, and the associated norms ·H, ·Q, ·H1 and ·Q1, respectively. We denote by vν and vτ the normal and the tangential component of v on Γ, respectively, given by vν=v·ν, vτ=v-vνν. For a regular function σ:Ω∪Γ→Sd 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. LetX be a Banach space and X∗ its dual. For a function φ:X→Rd, thedirectional derivative of φ at x∈X in the direction v∈X is defined by(5)φ′x;v=limλ↓0⁡φx+λv-φxλ whenever this limit exists. The Clarke generalized directional derivative of a locally Lipschitz function φ:X→Rd at the point x∈X in the direction v∈X is defined by(6)φ0x;v=limsupy→x,λ↓0⁡φy+λv-φyλ. The Clarke subdifferential of φ at x is a subset of X∗ given by(7)∂φx=ζ∈X∗:φ0x;v≥ζ,vX∗×X∀v∈X.A locally Lipschitz function φ:X→Rd is said to beregular (in the sense of Clarke) at x∈X if, for all v∈X, the directional derivative φ′(x;v) exists and φ′(x;v)=φ0(x;v). The function φ is regular (in the sense of Clarke) on X if it is regular at every point x∈X.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.
LetT>0 be given. For a positive integer N we define k=T/N. Assume that {gn}n=1N and {en}n=1N are two sequences of nonnegative numbers satisfying(8)en≤c~gn+c~∑j=1nkej,n=1,…,N for a positive constant c~ independent of N or k. Then there exists a positive constant c, independent of N or k, such that(9)max1≤n≤Nen≤cmax1≤n≤Ngn.


## 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 Γ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 u¨ 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,T] be the time interval of interest with T>0. The body is clamped on Γ1×(0,T) and, therefore, the displacement field vanishes there. A volume force of density f0 acts in Ω×(0,T), and surface tractions of density f2 act on Γ2×(0,T). In this study, we consider the dynamic contact with friction between a viscoelastic body and a foundation on Γ3×(0,T). The material’s behavior is modelled by a viscoelastic constitutive law of the form(10)σ=Cεu˙+Eεu.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,(11)-σν=puν+qu˙ν.Here p represents the positive normal compliance function such that p(r)=0 for r≤0 and q represents the positive normal damped response function such that q(r)=0 for r≤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(12)στ≤μu˙τpuν+qu˙ν,-στ=μu˙τpuν+qu˙νu˙τu˙τifu˙τ≠0.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,T)→Rd and a stress field σ:Ω×(0,T)→Sd such that(13)σt=Cεu˙t+EεutinΩ×0,T,(14)Div⁡σt+f0t=ρu¨tinΩ×0,T,(15)ut=0onΓ1×0,T,(16)σtν=f2tonΓ2×0,T,(17)-σνt=puνt+qu˙νtonΓ3×0,T,(18)στt≤μu˙τtpuνt+qu˙νt,-στt=μu˙τtpuνt+qu˙νtu˙τtu˙τtifu˙τt≠0onΓ3×0,T,for all t∈R+ and, moreover,(19)u0=u0u˙0=u1inΩ.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 u0 and u1 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 ProblemP, we consider a closed subspace of H1(20)V=v=vi∈H1∣v=0onΓ1.Completeness of the space(V,·V) follows from the assumption meas(Γ1)>0, which allows the use of Korn’s inequality, for some constant CK, depending only on Ω and Γ1,(21)εvQ≥CKvH1∀v∈V.On V, we use the inner product given by(22)u,vV=εu,εvQ∀u,v∈V,and let ·V be the associated norm; that is,(23)vV=εvQ∀v∈V.It follows from (21) and (23) that ·H1 and ·V are equivalent norms on V and therefore (V,·V) is a real Hilbert space. Also, note that since(24)vH1=vH2+vV21/2∀v∈H1, we obtain that(25)vV≥CKvH∀v∈V.The duality pairing betweenV and V∗ is denoted by 〈·,·〉. Identifying H with its dual, we have an evolution triple V⊂H⊂V∗ with dense, continuous, and compact embedding. By the Sobolev trace theorem and by (25), there exists a constant c0 depending only on the domains Ω, Γ1, and Γ3 such that(26)vL2Γ3;Rd≤c0vV∀v∈V.By (26) there exists a continuous trace operator γ:V→L2(Γ3;Rd) and for the function v∈V we still denote by v its trace γv. In what follows, we need the spaces V=L2(0,T;V), H=L2(0,T;H), 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 vW=(vV2+vV∗2)1/2, the space W becomes a separable Hilbert space. We also have W⊂V⊂H⊂V∗ and the following inequality holds, where ce is a constant,(27)vV∗≤cevL20,T,L2Γ3;Rd∀v∈V.It is well known that the embedding W⊂C([0,T];H) and embedding {w∈V∣w˙∈W}⊂C([0,T];V) are continuous. Moreover, choosing η∈0,1/2 we introduce the space Z=V∩H1-η(Ω;Rd) equipped with H1-η norm. Note that V⊂Z and Z⊂H with both embeddings being compact, and hence H⊂Z∗⊂V∗ where both embeddings are again compact.We introduce now Ehrling’s lemma, proved in [21].Lemma 2.
LetX0, X, and X1 be three Banach spaces such that (28)X0⊂X⊂X1, with the injection of X into X1 being continuous and the injection of X0 into X being compact. Then, for each ε>0, there exists some constant C(ε) such that(29)vX≤εvX0+CεvX1,∀v∈X0.Using this lemma, we get the following inequality:(30)vL2Γ3,Rd2≤εvV2+CεvH2,∀v∈V.In the study of the contact problem, we assume the following properties on the data.H(C)_.
The viscosity operatorC:Ω×Sd→Sd satisfies(31)aC·,εismeasurableonΩ∀ε∈Sd,bCx,ε1-Cx,ε2:ε1-ε2≥c2ε1-ε22∀ε1,ε2∈Sd,a.e.x∈Ωwithc2>0,cCx,ε1-Cx,ε2≤LCε1-ε2∀ε1,ε2∈Sd,a.e.x∈ΩwithLC>0.H(E)_.
The elasticity operatorE:Ω×Sd→Sd satisfies(32)aE·,εismeasurableonΩ∀ε∈Sd,bEx,ε1-Ex,ε2:ε1-ε2≥0∀ε1,ε2∈Sd,a.e.x∈Ω,cEx,ε1-Ex,ε2≤LEε1-ε2∀ε1,ε2∈Sd,a.e.x∈ΩwithLE>0.H(p)_.
The normal compliance functionp:Γ3×R→R satisfies(33)aThereexistsLp>0suchthatpx,u-px,v≤Lpu-v∀x∈Γ3,bThereexistsc3>0suchthatpx,u≤c3∀x,u∈Γ3×R,cpx,u=0∀u≤0a.ex∈Γ3.H(q)_.
The normal damped response functionq:Γ3×R→R satisfies(34)aThereexistsLq>0suchthatqx,u-qx,v≤Lqu-v∀x∈Γ3,bThereexistsc4>0suchthatqx,u≤c4∀x,u∈Γ3×R,cqx,u=0∀u≤0a.ex∈Γ3.H(μ)_.
The friction boundμ:[0,∞)→R satisfies(35)aμiscontinuous,bμs≤c,c>0,cμs1-μs2≥-λs1-s2∀s1>s2≥0withλ>0.The assumption (35)(c) represents a one-side Lipschitz condition, which allows the function μ(·) to decrease with a rate not faster than λ.H(f)_.
The force and the traction densities satisfy(36)f0∈L20,T;H,f2∈L20,T;L2Γ2;Rd.Also, we assume that the initial values satisfy the following.H0_.
u0∈V and u1∈H.Using the Clarke subdifferential (7), the friction condition (18) can be expressed in another form. Indeed, define a function j:Rd→R by(37)jη=∫0ημsds∀η∈Rd.Then, assuming H(μ)_(a)–(c), the condition (18) is equivalent to the following subdifferential inclusion: (38)-στ∈puν+qu˙ν∂ju˙τonΓ3×0,T,where ∂j(η) denotes the Clarke subdifferential of j at the point η∈Rd.Properties of the functionj are summarized in the next lemma.Lemma 3.
If the assumptionsH(μ)_(a)-(b) hold, then the function j defined by (37) is locally Lipschitz, and(39)ξ·u≥0∀u∈Rd,ξ∈∂ju,ξ≤c∀u∈Rd,ξ∈∂ju. If furthermore the assumption H(μ)_(c) holds, then we have(40)ξ1-ξ2·u1-u2≥-λu1-u22∀u1,u2∈Rd,ξi∈∂jui,i=1,2.Proof.
First, we have to proof thatj is locally Lipschitz. In order to do so, let η∈Rd and r>0. For η1,η2∈B(η,r), we get(41)jη1-jη2≤maxz∈Bη,r⁡μzη1-η2. With H(μ)_(a) and Proposition 5.6.28(ii) in [22], we have the following characterization of the Clarke subdifferential ∂j:(42)∂jη=B0,μ0forη=0,μηηηforη≠0. Then, the other properties follow straightforwardly.We define the operatorsA:V→V∗, B:V→V∗, functional J:L2(Γ3;Rd)→R, and f:(0,T)→V∗, defined by(43)Au,vV∗×V=∫ΩCεuεvdxforu,v∈V,Bu,vV∗×V=∫ΩEεuεvdxforu,v∈V,Jv=∫Γ3jvdΓforv∈L2Γ3;Rd,ft,vV∗×V=f0t,vH+f2t,vL2Γ2;Rdforv∈V.IfH(C)_ holds, the operator A has the following properties:(44)aA∈LV;V∗issymmetric,bThereexistsLC>0suchthatAv,vV∗×V≤LCvV2∀v∈V,cThereexistsc2>0suchthatAv,vV∗×V≥c2vV2∀v∈V.Under the assumptions H(E)_, the operator B satisfies the following properties:(45)aB∈LV;V∗issymmetric,bThereexistsLE>0suchthatBv,vV∗×V≤LEvV2∀v∈V,cBv,vV∗×V≥0∀v∈V.Assuming H(μ)_, the functional J:L2(Γ3;Rd)→R is locally Lipschitz, and we have the following inequalities on ξ∈∂J(v):(46)aξL2Γ3;Rd≤c∀ξ∈∂Jv,bξ·v≥0∀ξ∈∂Jv,cξ1-ξ2,v1-v2L2Γ3;Rd≥-λv1-v2L2Γ3;Rd2∀ξi∈∂Jvi,i=1,2,∀v1,v2∈L2Γ3;Rd.Now, multiplying the equation of motion by the test functionv, integrating over Ω×(0,T), using the Green formula and the definition of the operators, we obtain the following variational formulation of Problem P.Problem P V. Find a displacement field u∈V with u˙∈W and a friction stress field ξ∈L2(0,T,L2(Γ3;Rd)) such that, for a.e. t∈(0,T) and for all v∈V,(47)ρu¨t+Au˙t+But-ft,vV∗×V+puνt,vνL2Γ3;R+qu˙νt,vνL2Γ3;R+puνt+qu˙νtξt,vτL2Γ3;Rd=0,with(48)ξt∈S∂ju˙τt2onΓ3×0,T, where S∂j(u˙τ(t))2 is the set of all L2 selections of ∂j(u˙τ(t)) and(49)u0=u0,u˙0=u1.Finally, we formulate the result concerning the existence and uniqueness of solution to ProblemPV. In order to keep the paper with a reasonable length, we skip the proof of this result.Theorem 4.
Assume thatH(C)_, H(E)_, H(p)_, H(q)_, H(μ)_, H(f)_, and H0_ hold. Then, Problem PV 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 ProblemPV in order to bound the error of the fully discrete solutions. To this end, on a finite time interval [0,T] with T>0 given, we consider a positive integer N and we define the time step-size k=T/N and the time nodal points tn=nk and 0≤n≤N.We assume thatΩ is a polyhedral domain and we consider a regular family of finite element partition of Ω¯. Let Vh⊂V be the associate finite element space of continuous piecewise linear functions which vanish on Γ1. Here h>0 denotes the spatial discretization parameter.In what follows, we denotefn=f(tn) and un=u(tn). For a sequence uhk={unhk}n=0N, we use the notation δnunhk=(unhk-un-1hk)/k, n=1,…,N, for the backward divided differences, as well as the additional notation δuhk={δnuhk}n=1N. Let u0h,u1h∈Vh be suitable approximations of u0, u1 characterized by(50)u0h-u0,vhV=0,u1h-u1,vhH=0for all vh∈Vh. It is easy to observe that(51)u0hV≤u0V,u1hH≤u1H. The fully discrete approximation of Problem PV is the following.Problem P V h k. Find a velocity field {wnhk}n=0N⊂Vh and a friction stress field {ξnhk}n=0N⊂L2(Γ3;Rd) such that(52)ρδwnhk+Awnhk+Bunhk-fn,vhV∗×V+punνhk,vνhL2Γ3;R+qwnνhk,vνhL2Γ3;R+punνhk+qwnνhkξnhk,vτhL2Γ3;Rd=0∀vh∈Vh,with(53)ξnhk∈S∂jwnτhk2a.e.onΓ3,n=1,…,N,T=Nk,u0hk=u0h,w0hk=u1h. We consider the following discrete displacement(54)unhk=u0hk+∑j=1nkwjhk.Remark 5.
From now on, the constant denoted byC may differ from line to line.In what follows, our aim is to analyze the error between the solution of ProblemPV in displacement and its fully discrete approximation, solution of Problem PVhk. 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.
Assumek is sufficiently small. Then, under the assumptions H(C)_, H(E)_, H(p)_, H(q)_, H(μ)_, H(f)_, and H0_ and the regularity u∈C1(0,T;[H2(Ω)]d)∩H3(0,T;H) and u˙τ∈C(0,T;H2(Γ3)d) we have the optimal order error estimate(55)max1≤n≤N⁡un-unhkV+wn-wnhkH≤Ch+k, with C a constant independent of h and k.Proof.
We takev=vh in (47) and we combine it with (52) in order to get(56)ρw˙n-δwnhk+Awn-wnhk+Bun-unhk,vhV∗×V+punν-punνhk,vνhL2Γ3;R+qwnν-qwnνhk,vνhL2Γ3;R+punν+qwnνξn-punνhk+qwnνhkξnhk,vτhL2Γ3;Rd=0.Now, taking vh=wnhk in (56) and combining these 2 inequalities, we obtain(57)ρw˙n-δwnhk+Awn-wnhk+Bun-unhk,vhV∗×V+punν-punνhk,vνhL2Γ3;R+qwnν-qwnνhk,vνhL2Γ3;R+punν+qwnνξn-punνhk+qwnνhkξnhk,vτhL2Γ3;Rd=ρw˙n-δwnhk+Awn-wnhk+Bun-unhk,wnhkV∗×V+punν-punνhk,wnνhkL2Γ3;R+qwnν-qwnνhk,wνhL2Γ3;R+punν+qwnνξn-punνhk+qwnνhkξnhk,wnτhkL2Γ3;Rd.After reformulation of (57) under the form (·,wn-wnhk)=(·,wn-vh), we have(58)ρw˙n-δwnhk+Awn-wnhk+Bun-unhk,wn-wnhkV∗×V+punν-punνhk,wnν-wnνhkL2Γ3;R+qwnν-qwnνhk,wnν-wnνhkL2Γ3;R+punν+qwnνξn-punνhk+qwnνhkξnhk,wnτ-wnτhkL2Γ3;Rd=ρw˙n-δwnhk+Awn-wnhk+Bun-unhk,wn-vhV∗×V+punν-punνhk,wnν-vνhL2Γ3;R+qwnν-qwnνhk,wnν-vνhL2Γ3;R+punν+qwnνξn-punνhk+qwnνhkξnhk,wnτ-vτhL2Γ3;Rd.
From that, we deduce that(59)ρw˙n-δwnhk+Awn-wnhk+Bun-unhk,wn-wnhkV∗×V+punν-punνhk,wnν-wnνhkL2Γ3;R+qwnν-qwnνhk,wnν-wnνhkL2Γ3;R+punν+qwnνξn-punνhk+qwnνhkξn,wnτ-wnτhkL2Γ3;Rd+punνhk+qwnνhkξn-punνhk+qwnνhkξnhk,wnτ-wnτhkL2Γ3;Rd=ρw˙n-δwnhk+Awn-wnhk+Bun-unhk,wn-vhV∗×V+punν-punνhk,wnν-vνhL2Γ3;R+qwnν-qwnνhk,wnν-vνhL2Γ3;R+punν+qwnνξn-punνhk+qwnνhkξnhk,wnτ-vτhL2Γ3;Rd.
Therefore, after some elementary manipulations, it is easy to see that(60)ρδwn-δwnhk,wn-wnhkH+Awn-wnhk,wn-wnhkV∗×V+punνhk+qwnνhkξn-punνhk+qwnνhkξnhk,wnτ-wnτhkL2Γ3;Rd=ρw˙n-δwn,wnhk-wnH+ρδwn-δwnhk,wn-vhHρw˙n-δwn,wn-vhH+Awn-wnhk,wn-vhV∗×V+Bun-unhk,wnhk-wn+wn-vhV∗×V+punν-punνhk,wnνhk-wnν+wnν-vνhL2Γ3;R+qwnν-qwnνhk,wnνhk-wnν+wnν-vνhL2Γ3;R+punν+qwnνξn-punνhk+qwnνhkξnhk,wnτ-vτhL2Γ3;Rd+punνhk+qwnνhkξn-punν+qwnνξn,wnτ-wnτhkL2Γ3;Rd.
Remark  7. From now on, the constant denoted by C may differ from line to line.
It is easy to see that(61)2a-b,a=a-b2+a2-b2. Then, we take a=wn-wnhk and b=wn-1-wn-1hk(62)ρ2kwn-wnhkH2-wn-1-wn-1hkH2≤ρδwn-δwnhk,wn-wnhkH.From the coercivity of the operator A in (44), we get(63)c2wn-wnhkV2≤Awn-wnhk,wn-wnhkV∗×V.The operator A is Lipschitz continuous, from (44):(64)Awn-wnhk,wn-vhV∗×V≤LCwn-wnhkVwn-vhV.The operator B is also Lipschitz continuous, from (45):(65)Bun-unhk,wnhk-vhV∗×V≤LEun-unhkVwnhk-vhV.With (26), (33)(a), and (34)(a), we have(66)punν-punνhk+qwnν-qwnνhk,wnνhk-vνhL2Γ3;R≤c0Lpunν-unνhkV+Lqwnν-wnνhkL2Γ3;Rwnνhk-vνhL2Γ3;R≤c0Lpun-unhkV+Lqwn-wnhkL2Γ3;Rwnhk-wnL2Γ3;Rd+c0wn-vhV.From (30), (33)(b), (34)(b), and (46)(c), we get(67)punνhk+qwnνhkξn-punνhk+qwnνhkξnhk,wnτ-wnτhkL2Γ3;Rd≥-λc3+c4wnτ-wnτhkL2Γ3;Rd2≥-λc3+c4εwn-wnhkV2+Cεwn-wnhkH2.Using now (26), (33)(a), and (46)(a), we have(68)punνξn-punνhkξnhk≤punνξn-punνhkξn+punνhkξn-punνhkξnhk≤c0Lpcun-unhkV+2c3c.Then, with (26), (34)(a), and (46)(a), we have(69)qwnνξn-qwnνhkξnhk≤qwnνξn-qwnνhkξn+qwnνhkξn-qwnνhkξnhk≤c0Lqcwn-wnhkV+2c4c.Therefore, using this time only (26) and (33)(a) we deduce that(70)punνhkξn-punνξn≤c0Lpcun-unhkV,and from (26), (30), and (34)(a), we obtain(71)qwnνhkξn-qwnνξn≤Lqcεwn-wnhkV+Cεwn-wnhkH.
Thus, from (62)–(71), we get(72)ρ2kwn-wnhkH2-wn-1-wn-1hkH2+c2-λc3+c4+Lqcεwn-wnhkV2≤ρδwn-δwnhk,wn-vhH+ρw˙n-δwnHwn-wnhkH+wn-vhH+LCwn-wnhkVwn-vhV+LEun-unhkVwn-vhV+wn-wnhkV+c0Lpun-unhkV+Lqwn-wnhkL2Γ3;Rdwnhk-wnL2Γ3;Rd+c0wn-vhV+2c3+c4cwnτ-vτhL2Γ3;Rd+c0Lqcwn-wnhkVwnτ-vτhL2Γ3;Rd+c0Lpcun-unhkVwnτ-vτhL2Γ3;Rd+c0wnτ-wnτhkV+LqcCεwn-wnhkH2.
Then, from (25), (26), (30), the inequalities ab≤2a2+2b2 and ab≤εa2+b2/4ε, and (72), we obtain(73)ρ2kwn-wnhkH2-wn-1-wn-1hkH2+c2-λc3+c4+Lqcεwn-wnhkV2≤ρδwn-δwnhk,wn-vhH+Cw˙n-δwnH2+un-unhkV2+wn-vhV2+εwn-wnhkV2+Cwnτ-vτhL2Γ3;Rd+Cwnτ-vτhL2Γ3;Rd2+Cεwn-wnhkH2,with C a positive constant that may change from line to line.
Remark 8. Note that if ε is small enough, c2-(λ(c3+c4)+Lqc+1)ε>0.
Now, we replacen by j and sum over j from 1 to n to obtain(74)wn-wnhkH2-w0-w0hkH2+2ρc2-λc3+c4+Lqc+1ε∑j=1nkwj-wjhkV2≤2∑j=1nkδwj-δwjhk,wj-vjhH+C∑j=1nkw˙j-δwjH2+uj-ujhkV2+wj-vjhV2+C∑j=1nkwjτ-vjτhL2Γ3;Rd+C∑j=1nkwjτ-vjτhL2Γ3;Rd2+Cε∑j=1nkwj-wjhkH2.
We consider the term∑j=1nk((δwj-δwjhk),wj-vjh)H; we have(75)∑j=1nkδwj-δwjhk,wj-vjhH=∑j=1nwj-wjhk-wj-1-wj-1hk,wj-vjhH=∑j=1nwj-wjhk,wj-vjhH-∑j=0n-1wj-wjhk,wj+1-vj+1hH≤ε′wn-wnhkH2+Cwn-vnhH2+Cw0-w0hkH2+Cw1-v1hH2+∑j=1n-1wj-wjhkHwj-vjh-wj+1-vj+1hH≤ε′wn-wnhkH2+Cwn-vnhH2+Cw0-w0hkH2+Cw1-v1hH2+∑j=1n-1k4wj-wjhkH2+∑j=1n-11kwj-vjh-wj+1-vj+1hH2,with ε′<1/2. We recall the following classical inequality:(76)uj-ujhkV≤u0-u0hV+∑p=1jkwp-wphkV+Ij,where Ij is defined by the following relation(77)Ij=∫0tjwsds-∑p=1jkwpV≤kuH20,T;V.Thus, since 2ab≤a2+b2,(78)
u
j
-
u
j
h
k
V
2
≤
u
0
-
u
0
h
V
2
+
∑
p
=
1
j
k
w
p
-
w
p
h
k
V
2
+
k
2
u
H
2
0
,
T
;
V
2
+
2
u
0
-
u
0
h
V
∑
p
=
1
j
k
w
p
-
w
p
h
k
V
+
2
u
0
-
u
0
h
V
k
u
H
2
0
,
T
;
V
+
2
∑
p
=
1
j
k
w
p
-
w
p
h
k
V
k
u
H
2
0
,
T
;
V
≤
3
u
0
-
u
0
h
V
2
+
3
∑
p
=
1
j
k
w
p
-
w
p
h
k
V
2
+
3
k
2
u
H
2
0
,
T
;
V
2
≤
C
u
0
-
u
0
h
V
2
+
j
∑
p
=
1
j
k
2
w
p
-
w
p
h
k
V
2
+
k
2
u
H
2
0
,
T
;
V
2
.Therefore,(79)∑j=1nkuj-ujhkV2≤CTu0-u0hV2+∑j=1nk2∑p=1jwp-wphkV2+k2uH20,T;V2.
By considering the constantc5=min⁡(1,(2/ρ)(c2-(λ(c3+c4)+Lqc+1)ε)), it is possible to bound from below wn-wnhkH2+(2/ρ)(c2-(λ(c3+c4)+Lqc+1)ε)∑j=1nkwj-wjhkV2 by c5(wn-wnhkH2+∑j=1nkwj-wjhkV2). Note that c5 is positive due to Remark 8.
Denote(80)en=wn-wnhkH2+∑j=1nkwj-wjhkV2,gn=∑j=1nkw˙j-δwjH2+wj-vjhV2+∑j=1nkwjτ-vjτhL2Γ3;Rd+u0-u0hV2+k2uH20,T;V2+wn-vnhH2+w0-w0hkH2+w1-v1hH2+∑j=1n-11kwj-vjh-wj+1-vj+1hH2+∑j=1nkwjτ-vjτhL2Γ3;Rd2.
Then, we can reformulate (74) as follows:(81)en≤Cgn+∑j=1nkej,forn=1,…,N.
From (81), we apply the discrete form of the Gronwall Lemma 1 to obtain that(82)max1≤n≤N⁡wn-wnhkH2+∑j=1nkwj-wjhkV2≤C∑j=1Nkw˙j-δwjH2+wj-vjhV2+max1≤n≤Nwnτ-vnτhL2Γ3;Rd+max1≤n≤Nwnτ-vnτhL2Γ3;Rd2+u0-u0hV2+k2uH20,T;V2+max1≤n≤Nwn-vnhH2+w0-w0hkH2+w1-v1hH2+∑j=1N-11kwj-vjh-wj+1-vj+1hH2.
Now, letujh∈Vh be the finite element interpolations of uj. It follows from [24] that (83)k∑j=1Nw˙j-δwjH2≤ck2wH20,T;H2,1k∑j=1N-1wj-vjh-wj+1-vj+1hH2≤ch2uH20,T;E2,u0-u0hE≤chu0H2Ω;Rd,u1-u1hH≤chu1H1Ω;Rd. 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 ProblemP, 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 m, the value of the friction coefficient μ(wτRd) is fixed to a given value depending on the tangential velocity solution wτ 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 Γ3. In the case presented below, the discretization is based on “node-to-rigid” contact element, which is composed by one node of Γ3 and one Lagrange multiplier node. For more details on the discretization step and Computational Contact Mechanics, we refer to [2, 6, 25–27]. The numerical solution of the nonsmooth nonconvex variational Problem P is based on the iterative scheme given in the following algorithmic lines characterized by the time stepping loop of index n and the convexification loop of index m. The convexification loop is ended by the stopping criterion given by (85).Letϵconv>0, and let w0 be given.Forn=0,1,…,N, do time stepping loop.Form=0,1,…, do convexification loop.
ProblemPconv. Find a velocity field {wnhk,m+1}n=0N⊂Vh and a friction stress field {ξnhk,m+1}n=0N⊂L2(Γ3;Rd) such that (84)δwnhk,m+1+Anwnhk,m+1+Bunhk,m+1-fn,vhV∗×V+punνhk,m+1,vνhL2Γ3;R+qwnνhk,m+1,vνhL2Γ3;R+punνhk,m+1+qwnνhk,m+1ξnhk,m+1,vτhL2Γ3;Rd=0,∀vh∈Vh,ξnhk,m+1∈S∂jwτhk,m2onΓ3,T=Nk,untilunhk,m+1-unhk,mV≤ϵconvunhk,mV, and (85)ξnhk,m+1-ξnhk,mL2Γ3;Rd≤ϵconvξnhk,mL2Γ3;Rd.We recall that the discrete displacement is given by unhk=u0hk+∑j=1n‍kwjhk.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, Ω=(0,L)×(0,l)⊂R2 with L>0, l>0 and(86)Γ1=0×0,l∪L×0,l,Γ2=0,L×l,Γ3=0,L×0. 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.Figure 1
Reference configuration of the two-dimensional body.On the partΓ1=({0}×[0,l])∪({L}×[0,l]), the body is clamped and, therefore, the displacement field vanishes there. Vertical tractions act on the part [0,L]×{l} 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 Γ3=[0,L]×{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 tensorE given by (87)Eταβ=Eκ1+κ1-2κτ11+τ22δαβ+E1+κταβ,1≤α,β≤2,∀τ∈S2. Here, E 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 (88)pr=min⁡cνr+,c3,qr=min⁡dνr+,c4, where c3 and c4 are very large and r+=max⁡(0,r). Therefore, the normal compliance and normal damped response used here are defined by the following functions (89)pr=cνr+,qr=dνr+. For the coefficient of friction, we choose a function μ:[0,∞)→R+ of the form(90)μu˙τ=a-b·e-αu˙τ+b,with a,b,α>0 and a≥b. 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 a to the limit value b. And, for this reason, the corresponding friction law can be characterized as being nonmonotone. Since the function μ(·) is a contraction on [0,∞), we have α(a-b)<λ, 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 f2, with x, the abscissa of a point on the boundary (91)f2x,t=0,-gx,t,gx,t=Fe-ψx-ct2,∀x∈0,L∀t∈0,T, where F,ψ,c>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:(92)
T
=
1
s
,
L
=
2
m
,
l
=
0.1
m
,
u
0
=
0
m
,
u
1
=
0
m
/
s
,
ρ
=
1000
k
g
/
m
3
,
E
=
100
G
P
a
,
κ
=
0.3
,
f
0
=
0,0
G
P
a
,
F
=
10
G
P
a
/
m
,
Ψ
=
50
m
-
1
,
c
m
a
s
s
e
=
2
m
/
s
,
c
3
=
10
8
G
P
a
/
m
,
c
4
=
10
8
G
P
a
/
m
,
c
ν
=
5
G
P
a
,
d
ν
=
0.1
G
P
a
·
s
,
a
=
1
,
b
=
0.5
,
α
=
2000
.In Figure2, we plot the deformed configuration as well as the interface forces on Γ3 at the moments t=0.2s,0.4s,0.6s, and 0.8s. Note that at a given instant t=t0, 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 cmasse. This is why, with the choice cmasse=L/T, the mass proceeds from left to right through the whole boundary Γ2.Figure 2
Time evolution of the deformed meshes and contact forces onΓ3.Next, we study the influence of the parameterψ on the normal contact stresses on Γ3. At a given instant t=t0, since the mass is located at x0=cmasset0, the forces decrease exponentially with the distance between the node considered on the contact boundary and the node at the abscissa x0. Here, ψ describes the decay rate of g. 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 ψ.Figure 3
Deformed meshes and contact forces onΓ3 for ψ=5 and ψ=20, without friction.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 (93)vE≔ρv¨t+Bvt,vV∗×V1/2which is equivalent to the norm vV. Since the true solution u is not available, we use instead the numerical solution corresponding to fine discretization of Ω as the “reference” solution uref in computing the solution errors. Here, the numerical solution with h=1/128 and k=1/128 is taken to be the “reference” solution uref. 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 4
Numerical errors.



---

*Source: 1562509-2016-10-17.xml*