Advances in Variational and Hemivariational

Weimin Han Stanislaw Migórski Mircea Sofonea Editors

Advances in Variational and Hemivariational Inequalities with Applications Theory, Numerical Analysis, and Applications November 6, 2014

Springer

Contents

1

A Hyperelastic Dynamic Frictional Contact Model with Energy-Consistent Properties M. Barboteu, D. Danan and M. Sofonea . . . . . . . . . . . . . . . . . . 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Notation and physical setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 A specific frictional contact model . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Mechanical problem and variational formulation . . . . . . . . . . . . . 1.4.1 Mechanical problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.2 Variational formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Discretization and variational approximation . . . . . . . . . . . . . . . . 1.6 Numerical solution with energy-consistent properties . . . . . . . . . 1.6.1 Usual discrete energy-conserving framework . . . . . . . . . . . 1.6.2 An improved energy-consistent approach . . . . . . . . . . . . . 1.6.3 Analysis of the discrete energy evolution . . . . . . . . . . . . . 1.7 Numerical experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7.1 Impact of a linear elastic ball against a foundation . . . . . 1.7.2 Impact of a hyperelastic ring against a foundation . . . . . 1.8 References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1 1 3 6 8 8 9 10 13 13 14 15 17 18 21 25

1 A Hyperelastic Dynamic Frictional Contact Model with Energy-Consistent Properties M. Barboteu, D. Danan and M. Sofonea

Abstract. In this chapter we present an energy-consistent numerical model for the dynamic frictional contact between a hyperlastic body and a foundation. Our contribution has two traits of novelty. The first one arises from the specific frictional contact model we consider, which provides intrinsic energy-consistent properties. The contact is modeled with a normal compliance condition of such a type that the penetration is limited with unilateral constraint and friction is described with a version of Coulomb’s law of dry friction. The second trait of novelty consists in the construction and the analysis of an energy-consistent scheme, based on recent energy-controlling time integration methods for nonlinear elastodynamics. Some numerical results for representative impact problems are provided. They illustrate both the specific properties of the contact model and the energy-consistent properties of the numerical scheme. Keywords. Contact and friction, Normal compliance, Unilateral constraint, Coulomb friction, Nonlinear elastodynamics, Hyperelasticity, Time integration schemes, Energy-conserving algorithms. AMS Classification. 74M15, 74M20, 74M10, 74B20, 74H15, 74S30, 49M15, 90C53.

1.1 Introduction An important topic concerning the modelling of dynamic frictional contact problems is the enforcement of suitable contact interface conditions with M. Barboteu, D. Danan and M. Sofonea Laboratoire de Mathématiques et Physique, Université de Perpignan Via Domitia, 52 Avenue Paul Alduy, 66860 Perpignan, France, email: [email protected], [email protected] and [email protected]

1

2

M. Barboteu, D. Danan and M. Sofonea

energy-consistent properties. This leads to consider acceptable physical models and to elaborate numerical schemes with long term time integration accuracy and stability properties. During the last twenty years, many works have been devoted to the construction of energy-conserving methods for elastodynamic contact problems, see [2, 4, 6, 19, 20, 22, 25, 29, 30], for instance. Consistent energy dissipation extensions can be found in [4, 6, 22, 30] and [5], in the study of models with friction and viscosity, respectively. In this current work we consider new frictional contact conditions between an hyperelastic body and a foundation, which take into account the asperities of the contact surface and some energy-consistent properties, as well. More precisely, we consider a specific frictional contact model which provides intrinsic energy-consistent properties, characterized by a conserving behaviour for frictionless impacts and admissible dissipation for friction phenomena. In our model the contact is described with a normal compliance condition of such a type that the penetration is limited with unilateral constraint. This penetration can be assimilated to the flattening of the asperities on the contact surface of the foundation. The associated model of friction, which represents a generalization of the models introduced in [9, 10], is based on a version of Coulomb’s law in which we assume that the coefficient of friction depends on both the depth of the penetration and the tangential slip rate. In order to provide numerical energy-consistent properties related to the specific frictional contact, a crucial issue is to focus on computational aspects with long term time integration accuracy and stability properties. To this end, we consider an energy-consistent scheme based on recent energycontrolling time integration methods for nonlinear elastodynamics, developed in [2, 4, 6, 19, 20, 22, 25, 30]. In particular, we use a Newton continuation method and augmented Lagrangian arguments, already used in [4]. Next, to obtain additional energy conserving properties, we combine the specific penalized methods considered in [4, 18, 20], with the procedure of equivalent mass matrix introduced in [25]. We also provide some numerical experiments of representative impact problems which illustrate, in terms of conservation of energy, the good properties of both the frictional contact model and the numerical scheme. More precisely, we compare our numerical results with results obtained by using time integration schemes used in [4, 20, 25] and we discuss issues related to energy conservation properties. The rest of the chapter is structured as follows. In Section 1.2 we introduce the notation as well as some preliminary material used for the physical setting of hyperelastic contact problems. In Section 1.3 we describe the specific contact model, including the associated Coulomb’s law of friction. In Section 1.4, we present the classical formulation of the dynamic frictional contact problem and derive its variational formulation. Then, the Section 1.5 we introduce the fully discrete approximation of the problem. Section 1.6 is devoted to the analysis of the energy-consistent approach used to solve nonlinear elastodynamic frictional contact problems. Thus, after presenting the usual energyconserving frameworks used, we focus on the analysis of the discrete energy

A Frictional Contact Model with Energy Consistent Properties

3

evolution of the method. Finally, in Section 1.7 we present numerical results in the study of two representative two-dimensional examples with linear elastic and hyperelastic materials, respectively.

1.2 Notation and physical setting In this section we present the notation we shall use and some preliminary material. Everywhere in this chapter we use the notation N for the set of positive integers and R− will represent the set of non positive real numbers, i.e. R− = (−∞, 0]. We denote by Md the space of second order tensors on Rd or, equivalently, the space of square matrices of order d. The inner product and norm on Rd and Md are defined by u · v = u i vi , Π : τ = Πji τij ,

1

∀ u, v ∈ Rd ,

kvk = (v · v) 2 kτ k = (τ : τ )

1 2

∀ Π, τ ∈ Md .

Let Ω ⊂ Rd (d = 1, 2, 3) be an open bounded connected set with a Lipschitz boundary Γ . We use the notation x = (xi ) for a typical point in Ω ∪ Γ and we denote by ν = (νi ) the outward unit normal at Γ . 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. An index that follows a comma represents the partial derivative with respect to the corresponding component of the spatial variable, e.g. ui,j = ∂ui /∂xj . We consider the spaces V = {v ∈ H 1 (Ω; Rd ) : v = 0 on Γ1 }, H = L2 (Ω; Rd ). These are real Hilbert spaces endowed with their standard inner products (u, v)V and (Π, τ )H and their associated norms k·kV and k·kH , respectively. Note that V ⊂ H ⊂ V ∗ is an evolution triple, with all embeddings being continuous, compact and dense. The duality pairing between V ∗ and V will be denoted by hu, viV ∗ ×V . For an element v ∈ V we still write v for its trace and we denote by vν and v τ the normal and tangential components of v on Γ , given by vν = v · ν and v τ = v − vν ν, respectively. Also, for a regular stress function Π we use the notation Πν and Π τ for its normal and the tangential components, i.e. Πν = (Πν) · ν and Π τ = Πν − Πν ν. Moreover, we recall that the divergence operator is defined by the equality Div Π = (Πij,j ) and, finally, the following Green’s formula holds: Z Z Z Π : ∇v dx + Div Π · v dx = Πν · v dΓ ∀ v ∈ V. (1.1) Ω

Ω

Γ

In the rest of the chapter we consider the time interval of interest [0, T ] with T > 0. We denote by t ∈ [0, T ] the time variable and, as already mentioned, x ∈ Ω ∪Γ will represent the spatial variable. In order to simplify the notation,

4


Fig. 1.1. A deformable body in dynamic contact with a foundation.

we do not indicate the dependence of the functions on x and t. Moreover, we use the dots above to represent the derivatives with respect to the time. We also use u for the displacement field and Π for the first Piola-Kirchoff stress tensor. We consider a dynamic contact problem between a deformable body and a foundation in the framework of finite deformations theory. The material’s behavior is described with a hyperelastic constitutive law. We recall that hyperelastic constitutive laws are characterized by the first Piola-Kirchoff tensor Π which derives from an internal hyperelastic energy density W (F), i.e. Π = ∂F W (F). Here F is the deformation gradient defined by F = I + ∇u and ∂F represents the differential with respect to the variable F, see [12] for details. The deformable body occupies the domain Ω with the boundary partitioned into three disjoint parts Γ 1 , Γ 2 and Γ 3 with Γ1 , Γ2 and Γ3 being relatively open. The part Γ1 is the part in which the displacement field is prescribed. A volume force of density f 0 acts in Ω × (0, T ), and we assume that a density of traction forces, denoted by f 2 , acts on the part Γ2 , i.e. σν = f 2

on

Γ2 × (0, T ).

On the part Γ3 the body can arrive in frictional contact with an obstacle, the so-called foundation, as shown in Figure 1.1. Various contact boundary conditions have been used to model contact phenomena, both in engineering and mathematical literature, see for instance [1, 14, 27, 17, 28, 33, 35, 37, 39, 40, 41, 42] and the references therein. One of the most popular is the Signorini condition, introduced in [38] which describes the contact with a perfectly rigid foundation. Expressed in terms of unilateral constraints for the displacement field, this condition leads to highly nonlinear and nonsmooth mathematical problems. The unilateral contact conditions


5

with a gap between a deformable body and a rigid foundation are given by uν − G ≤ 0,

Πν ≤ 0,

Πν (uν − G) = 0

on Γ3 × (0, T ),

(1.2)

where the gap function G measures the distance between a point on Γ3 and its projection onto the rigid obstacle. In the following, for simplicity, we consider that the gap function vanishes, i.e. we use condition (1.2) with G = 0. The Signorini contact condition is idealistic since a fondation is never perfectly rigid, due to the presence of microscopic asperities on its surface. Furthermore, it induces non negligible dissipation of energy during impacts, as explained in [2, 29]. For this reason, an unilateral contact condition expressed in terms of velocity field has been considered in the literature. Its form is given by u˙ ν ≤ 0,

Πν ≤ 0,

Πν u˙ ν = 0 on

Γ3 × (0, T ).

(1.3)

Even if it induces good properties of conservation of energy, condition (1.3) is not realistic, since it could lead to infinite penetrations. The contact with a deformable foundation is modelled by the so-called normal compliance contact condition. It assigns a reactive normal pressure that depends on the interpenetration of the asperities on the body’s surface and those on the foundation. The normal compliance contact condition was first introduced in [33, 37] in the study of dynamic contact problems with elastic and viscoelastic materials. A general expression for the normal compliance condition with a zero gap function is −Πν = p(uν )

on Γ3 × (0, T )

(1.4)

where p(·) is a nonnegative prescribed function which vanishes for negative argument. This condition can be viewed as a regularization of the Sigorini unilateral condition. It is obvious to see that the normal compliance condition is characterized by a non limited penetration. In the case of the frictional contact, the contact condition is usually associated to Coulomb’s law of dry friction, given by ( kΠ τ k ≤ −µΠν if u˙ τ = 0, on Γ3 × (0, T ). (1.5) ˙τ −Π τ = −µΠν ku u˙ k if u˙ τ 6= 0, τ

Here µ is a positive variable, the coefficient of friction. As noted in [36], the tribological laws (1.2) and (1.5) can be written in the form of subdifferential inclusions which derive from non-differentiable convex potentials: Πν ∈ ∂IR− (uν ) and

− Π τ ∈ −µΠν ∂ku˙ τ k on Γ3 × (0, T ).

Here ∂IR− and ∂ku˙ τ k denote the subdifferential of the indicator function IR− of the negative half-line of R and the subdifferential of the norm of the slip rate, ku˙ τ k, respectively. Furthermore, in the case of large deformations, we can refer to [14] to the definition of objective quantities related to contact mechanics.

6


1.3 A specific frictional contact model The purpose of this section is to present a specific frictional contact model which provide intrinsic energy-consistent properties. The contact is modeled with a normal compliance condition of such a type that the penetration is limited with unilateral constraint and the friction is modeled with a specific version of Coulomb’s law. In order to obtain energy conservation properties, the work of normal R contact reactions, denoted by Wc = Γ3 Πν u˙ ν dΓ , has to vanish. Therefore, for energy conservation purposes, as explained in [2, 29], the following persistency condition has to be added: u˙ ν Πν = 0 on

Γ3 × (0, T ).

(1.6)

This condition means that normal contact reactions can appear only during persistent contact. Note that the unilateral contact condition (1.3), expressed in terms of velocity, could lead to displacements which do not satisfy the non penetration condition. Furthermore, as explained in [4, 6], it is impossible to enforce, at a given moment, both the complementarity condition in displacement and in velocity. To overcome this drawback, in order to take into account the deformability of the foundation (arising from the existence of micro asperities on its surface), we consider the following normal compliance condition restricted by unilateral constraint:   Πν + p(uν ) ≤ 0, uν − g ≤ 0, on Γ3 × (0, T ). (1.7)  (Πν + p(uν ))(uν − g) = 0,

This condition was introduced for an elastic-visco-plastic problem in [23]. In this model, the contact follows a normal compliance condition with penetration but up to the limit g and then, when this limit is reached, the contact follows a Signorini-type unilateral condition with the gap g and without any aditionnal penetration in the foundation. We conclude from above that condition (1.7) models the contact with a foundation which is composed by a thin deformable layer of asperities of thickness g which covers a perfect rigid material. This contact model has two intrinsic advantages: the adequation with energy conservation properties during penetrations for the impact phase (−Πν = p(uν ) for 0 ≤ uν < g), on one hand, and the limitation of penetration into the foundation (Πν ≤ 0, uν −g ≤ 0, Πν (uν −g) = 0), on the other hand. Note that the energy conservation property during penetrations comes from the specific form of the normal compliance function p. Indeed, let us consider a normal compliance function p defined by p(uν ) = ru+ ν

with

u+ ν = max(0, uν ) = distR− (uν ),

(1.8)

where r is a positive constant which represents the deformability of the foundation. Then, with (1.8), the persistency condition is no longer necessary to


7

obtain energy conservation properties and the work Wc of contact reactions takes the form Z Z Z r d + 2 ru+ u ˙ dΓ = − Πν u˙ ν dΓ = − Wc = {uν } dΓ. (1.9) ν ν Γ3 2 dt Γ3 Γ3 If we consider the Gonzalez approach, [15], then the conservation of the energy for a contact system with a normal compliance condition of the form (1.8), is provided by the following energy assessment on [0, t]: Z tZ Z tZ g.u˙ dΓ f .u˙ dΩ + E(t) − E(0) = 0

−

r 2

Z

Γ3

Ω

0

Γ2

+ 2 (uν (t))2 − (u+ dΓ. ν (0))

(1.10)

Here E(t) represents the internal energy of the body Ω at time t and is defined by Z Z 1 ρu˙ 2 dΩ + W (F) dΩ. (1.11) E(t) = 2 Ω Ω

We refer to Section 1.6.2 or [20] for more details about this energy assessment properties for a normal compliance condition of the form (1.8). The same statement can be established for the angular and the linear momentum, as shown in [39]. We turn now to the friction conditions. Our goal is to consider a friction model suited to the previous contact conditions. To this end, we introduce a specific version of Coulomb’s law of dry friction in which the friction bound depends both on the depth of the penetration for 0 ≤ uν ≤ g and on the normal contact stress Πν . Therefore, we consider the following friction condition: ( kΠ τ k ≤ −µ(uν )Πν if u˙ τ = 0, on Γ3 × (0, T ). (1.12) ˙τ −Π τ = −µ(uν )Πν ku u˙ k if u˙ τ 6= 0, τ

Here µ denotes the coefficient of friction and is assumed to depend on the penetration uν as long as uν < g. When there is penetration, as far as the normal displacement does not reach the bound g (i.e. 0 ≤ uν < g), the contact is described with a normal compliance condition associated to the classical Coulomb’s law of dry friction with the friction bound µ(uν )p(uν ). Details on the normal compliance contact condition associated to Coulomb’s law of dry friction can be found in [17, 40, 41], for instance. After the complete flatenning of the asperities, i.e. when the normal displacement reaches the bound g, the magnitude of the normal stress is larger than p(g) and, moreover, friction follows a Coulomb’s law associated to unilateral contact, with the friction bound −µ(g)Πν . Note that the friction bound is characterized by the friction coefficient µ which depends on the depth of the penetration uν and on the size of the asperities g. In what follows, we consider two examples, namely

8


  0 for η ≤ 0, µ1 (η) = ηg µ0 for η ∈ (0, g),   µ0 for η ≥ g.

  0 for η ≤ 0, µ2 (η) = (2 − ηg )µ0 for η ∈ (0, g). (1.13)   µ0 for η ≥ g,

Here µ0 = µ(g) denotes a given coefficient of friction associated to the unilateral contact (uν = g). In the case of function µ1 , we remark that the friction bound increases with respect to the flattening of the asperities. In contrast, in the case of the function µ2 , the friction bound decreases with respect to the flattening of the asperities. Several experimental studies have demonstrated the dependence of the friction coefficient with respect to the normal compression load and the flatenning of asperities. This behavior is generated by the wear of asperities on the contact surfaces. References on this matter include [16, 31, 32], among others.

1.4 Mechanical problem and variational formulation 1.4.1 Mechanical problem With these preliminaries the classical formulation of hyperelastodynamic frictional contact problem is the following. Problem P M . Find a displacement field u : Ω × [0, T ] → IRd and a stress field Π : Ω × [0, T ] → Md such that

(

Π = ∂F W (F)

in Ω × (0, T ),

(1.14)

ρ¨ u − Div Π − f 0 = 0 u=0

in Ω × (0, T ), on Γ1 × (0, T ),

(1.15) (1.16)

Πν = f 2 uν ≤ g, Πν + p(uν ) ≤ 0, (uν − g)(Πν + p(uν )) = 0

on Γ2 × (0, T ),

(1.17)

on Γ3 × (0, T ),

(1.18)

on Γ3 × (0, T ),

(1.19)

in Ω.

(1.20)

kΠ τ k ≤ −µ(uν )Πν u˙ τ −Π τ = −µ(uν )Πν ku ˙ k τ

if if

u˙ τ = 0, u˙ τ 6= 0,

˙ u(0) = u0 , u(0) = u1

We recall that equation (1.14) is the hyperelastic constitutive law. Equation (1.15) represents the equation of motion in which ρ denotes the density of the material and is assumed to be constant, for the sake of simplicity. Conditions (1.16) and (1.17) are the displacement and traction boundary conditions, respectively. Finally, (1.20) represent the initial conditions in which u0 and u1 are the initial displacement and velocity, respectively. Next, conditions (1.18) and (1.19) represent the frictional contact conditions, already introduced in Section 1.3. We recall that (1.18) represents a contact condition with normal compliance contact and unilateral constraint,


9

in which the penetration is limited to the value g. Condition (1.19) represents a version of Coulomb’s law of dry friction in relation with the contact conditions (1.18). Note that the condition (1.18) is equivalent with −Πν ∈ p(uν ) + ∂I(−∞,g] (uν ) on Γ3 × (0, T ),

(1.21)

where ∂ represents the subdifferential operator in the sense of the convex analysis and IA denotes the indicator function of the set A ⊂ R. In the same way, we observe that the condition (1.19) is equivalent with −Π τ ∈ −µ(uν )Πν ∂ku˙ τ k on Γ3 × (0, T ).

(1.22)

In the rest of the Chapter, we will consider the frictional contact conditions in their subdifferential form (1.21), (1.22). 1.4.2 Variational formulation We now introduce a hybrid variational formulation of Problem P M in which the dual variables corresponding to Lagrange multipliers are related to the contact stress and the friction force. In this case, the Lagrange multipliers verify extended subdifferential inclusions derived from the pointwise subdifferential inclusions defined in (1.21) and (1.22). To this end we consider the trace spaces Xν = { vν |Γ3 : v ∈ V } and Xτ = { vτ |Γ3 : v ∈ V }, equipped with their usual norms. We denote by Xν′ and Xτ′ the duals of the spaces Xν and Xτ , respectively. Moreover, we denote by h·, ·iXν′ ,Xν and h·, ·iXτ′ ,Xτ the corresponding duality pairing mappings. To establish the variational formulation, we need additional notations. Thus, we consider the function f : (0, T ) → V ∗ and the operator B : V → V ∗ defined by hf (t), viV ∗ ×V = (f 0 (t), v)H + (f 2 (t), v)L2 (Γ2 ;Rd ) , Z hBu, viV ∗ ×V = Π(u) : ∇v dx

(1.23) (1.24)

Ω

for all t ∈ (0, T ), u, v ∈ V . For the contact conditions, we introduce a function ϕν : Xν → (−∞, +∞] and an operator L : Xν → Xν′ defined by Z ϕν (uν ) = I(−∞,g] (uν ) dΓ ∀ uν ∈ Xν , Γ3

L : Xν →

Xν′ ,

hLuν , vν iXν′ ,Xν =

Z

p(uν )vν dΓ

∀ uν , vν ∈ Xν .

Γ3

We note that, for all t ∈ (0, T ), condition (1.21) leads to the subdifferential inclusion −Πν (t) ∈ ∂ϕν (uν (t)) + Luν (t) in Xν′ . (1.25)

10


To reformulate the friction law, we introduce the function ϕτ : Xτ → (−∞, +∞] defined by Z ϕτ (uτ ) = ku˙ τ k dΓ ∀ uτ ∈ Xτ . Γ3

We note that for all t ∈ (0, T ), condition (1.22) leads to the subdifferential inclusion −Π τ (t) ∈ −µ(uν (t))Πν ∂ϕτ (uτ (t)) in Xτ′ . (1.26) Inclusions (1.25) and (1.26) suggest to introduce two new unknowns, the Lagrange multipliers, which represent the normal and tangential stresses on the contact surface, and which will be denoted in what follows by ξν and ξ τ , respectively. Therefore, multiplying the equation of motion (1.15) by the test function v − u, integrating the result over Ω × (0, T ) and using the Green formula (1.1) and the inclusions (1.25)–(1.26), we obtain the following hybrid variational formulation of the Problem P M , in terms of three unknown fields. Problem PV . Find a displacement field u : [0, T ] → V , a normal stress field ξν : [0, T ] → Xν′ and a tangential stress field ξ τ : [0, T ] → Xτ′ such that hρ¨ u(t) + Bu(t), viV ∗ ×V = hf (t), viV ∗ ×V +hξν (t), vν iXν′ ,Xν + hξ τ (t), vτ iXτ′ ,Xτ

(1.27) ∀ v ∈ V,

−ξν (t) ∈ ∂ϕν (uν (t)) + Luν (t),

(1.28)

−ξ τ (t) ∈ −µ(uν (t))ξν ∂ϕτ ((uτ (t)),

(1.29)

for all t ∈ [0, T ] and, moreover, u(0) = u0 ,

˙ u(0) = u1 .

(1.30)

We note that the existence and uniqueness of weak solution of Problem PV represents, at the best of our knowledge, an open mathematical problem. Nevertheless, the solvability of a dynamic viscoelastic frictional contact problem with a regularized version of the frictional contact conditions (1.28)– (1.29) has been established in [8], in the case of small deformations theory. We also recall that the question of weak solvability for several contact problems is discussed in detail in the books [17, 35, 40].

1.5 Discretization and variational approximation This section is devoted to the discretization of the variational problem PV , based on arguments similar to those used in [4, 5, 7, 8, 9, 10, 11]. First, we recall some preliminary material concerning the time discretizaT be the time step and define tion step. Let N be an integer, let k = N


tn = n k,

11

0 ≤ n ≤ N.

Below, for a continuous function f (t) with values in a function space, we use the notation fj = f (tj ), for 0 ≤ j ≤ N . In what follows, we consider a collection of discrete times {tn }N n=0 which define a uniform partition of the SN time interval [0, T ] = n=1 [tn−1 , tn ] with t0 = 0, tn = tn−1 + k and tN = T . Finally, for a sequence {wn }N n=1 , we denote the midpoint divided differences by 1 δwn− 21 = (wn − wn−1 )/k = (δwn + δwn−1 ), (1.31) 2 and, equivalently, we have δwn = −δwn−1 + k2 (wn − wn−1 ). In the rest of the paper, we use the notation 2n− 21 = 21 (2n + 2n−1 ), where 2n represents the approximation of 2(tn ). Note that the time integration scheme we use is based on the implicit second order midpoint rule given in (1.31). We now present some material concerning the spatial discretization step. Let Ω be a polyhedral domain. Moreover, let {T h } be a regular family of triangular finite element partitions of Ω that are compatible with the boundary decomposition Γ = Γ1 ∪ Γ2 ∪ Γ3 , i.e., if one side of an element T r ∈ T h has more than one point on Γ , then the side lies entirely on Γ1 , Γ2 or Γ3 . The space V is approximated by the finite dimensional space V h ⊂ V of continuous and piecewise affine functions, that is, V h = { v h ∈ [C(Ω)]d : v h |T ∈ [P1 (T r)]d ∀ T r ∈ T h , v h = 0 at the nodes on Γ1 }, where P1 (T r) represents the space of polynomials of degree less or equal to one in T r. For the discretization of the normal contact terms, we consider the space h Xνh = { vν| : vh ∈ V h } Γ 3

equipped with its usual norm. Let us consider the discrete space of piecewise constants Yνh ⊂ L2 (Γ3 ) related to the discretization of the normal stress ξν . Then, we note that the contact condition (1.28) leads to the following discrete subdifferential inclusion at time tn− 21 : hk hk −ξν hk n− 1 ∈ ∂ϕν (uν n− 1 ) + Luν n− 1 2

2

2

in

Yνh .

For the discretization of the tangential friction terms, let us consider the space Xτh = { v hτ|Γ

3

: vh ∈ V h }

equipped with its usual norm. We also consider the discrete space of piecewise constants Yτh ⊂ L2 (Γ3 )d related to the discretization of the friction density ξ τ . In a similar way, we note that the friction condition (1.29) leads to the following discrete subdifferential inclusion at time tn− 21

12

M. Barboteu, D. Danan and M. Sofonea hk h hk −ξ τ hk n− 1 ∈ −µ(uν )ξν n− 1 ∂ϕτ (uτ n− 1 ) in 2

2

2

Yτh .

More details about the discretization step can be found in [24, 26, 42]. Let uh0 ∈ V h and uh1 ∈ V h be finite element approximations of u0 and u1 , respectively. Then, using the previous notations and the midpoint scheme (1.31), the fully discrete approximation of the Problem PV at the time tn− 21 is the following. N h Problem PVhk . Find a discrete displacement field uhk = {uhk n }n=0 ⊂ V , a hk N h hk discrete normal stress field ξν = {ξν n }n=0 ⊂ Yν and a discrete tangential hk N h stress field ξ hk τ = {ξ τ n }n=0 ⊂ Yτ such that, for all n = 1, . . . , N , ρ hk h hk hk h h δuhk n − δun−1 + Bun− 21 , v iV ∗ ×V = hf n− 12 , v iV ∗ ×V k hk h h +hξν hk n− 1 , vν iYνh ,Xνh + hξ τ n− 1 , vτ iYτh ,Xτh 2

2

∀ v ∈ V h,

(1.32)

hk hk −ξν hk n− 1 ∈ ∂ϕν (uν n− 1 ) + Luν n− 1 , 2

2

(1.33)

2

hk hk hk −ξ τ hk n− 1 ∈ −µ(uν n− 1 )ξν n− 1 ∂ϕτ (uτ n− 1 ),

(1.34)

h uhk 0 = u0 ,

(1.35)

2

2

2

2

h δuhk 0 = u1 .

Note that the discrete frictional contact conditions (1.33) and (1.34) are considered at the time tn− 21 . As a consequence the solution uhk n does not verify the contact conditions at the desired time tn and the scheme does not control penetration. In order to overcome this drawback, several authors impose the contact conditions at the time tn , see [4, 5, 15, 20, 24, 30], for instance. Here we use an implicit backward divided difference for the discretization of the tangential velocity u˙ τ (t) given by δuτ n = (uτ n − uτ n−1 )/k, which leads to the following discrete problem. N Problem P¯ hk . Find a discrete displacement field uhk = {uhk ⊂ V h, a n } n=0

V

h N discrete normal stress field ξνhk = {ξν hk n }n=0 ⊂ Yν and a discrete tangential hk N hk h stress field ξ τ = {ξ τ n }n=0 ⊂ Yτ such that, for all n = 1, . . . , N , ρ hk h h hk hk h δuhk n − δun−1 + Bun− 21 , v iV ∗ ×V = hf n− 12 , v iV ∗ ×V k hk h h +hξν hk n , vν iYνh ,Xνh + hξ τ n , vτ iYτh ,Xτh

∀ v ∈ V h,

(1.36)

hk hk −ξν hk n ∈ ∂ϕν (uν n ) + Luν n ,

(1.37)

hk hk hk −ξ τ hk n ∈ −µ(uν n )ξν n ∂ϕτ (uτ n ),

(1.38)

h uhk 0 = u0

(1.39)

h in δuhk 0 = u1 .

Note that the specific discretization used in Problem P¯Vhk represents the starting point to develop improved energy-conserving algorithms for the solution of elastodynamic contact problems with long term time integration


13

accuracy and stability. Some details on the energy-conserving framework, can be found in the next section.

1.6 Numerical solution with energy-consistent properties 1.6.1 Usual discrete energy-conserving framework We start by recaling some preliminaries concerning the usual discrete energyconserving framework in the case without contact. In the rest of the section, to simplify the notation and the readability, we do not indicate the dependence of various variables with respect to the discretization parameters k and h, i.e., for example, we write u instead of uhk . In order to solve a nonlinear elastodynamic problem, we have to use adapted time integration schemes. When nonlinear dynamic problems are considered, the standard implicit schemes (θ-method, Newmark schemes, midpoint or HHT methods) lose their unconditional stability, as explained in [21, 28]. Therefore, there is a need to use implicit energy conservative schemes as those used in [3, 15, 18, 28, 39], which are appropriate, due to their long term time integration accuracy and stability. In all these methods the corresponding discrete mechanical conservation properties are satisfied. To establish these discrete energy conservative properties, one of the most used implicit time integration scheme is the second order midpoint scheme given by 2 δun = −δun−1 + (un − un−1 ). (1.40) k Moreover, according the time integration scheme of Gonzalez [15], the variational inequality (1.36) is characterized by the operator B defined by Z hBun− 12 , viV ∗ ×V = Π algo : ∇v dx for v ∈ V h , (1.41) Ω

algo

in which the discrete tensor Π is introduced in order to satisfied exact discrete energy properties. This tensor defined by Gonzalez in [15] takes the form  algo  Π = Fn− 12 Σalgo ,     ˜ W (1.42) ˜ (Cn ) − W ˜ (Cn−1 ) Σalgo = 2 ∂∂C (Cn− 12 ) + 2[W     ˜ ∆Cn−1  W − ∂∂C (Cn− 12 ) : ∆Cn−1 ] ∆Cn−1 :∆Cn−1 ,

where ∆Cn−1 = Cn − Cn−1 and Cn−1 = t Fn−1 Fn−1 . As shown in [12], ˜ (C). Then, using the axiom of frame indifference implies that W (F) = W the arguments in [15], it follows that (1.42) verifies exact energy conservation characterized by condition

14


˜ (Cn ) − W ˜ (Cn−1 ). Π algo : (∇un − ∇un−1 ) = W

(1.43)

For more details on standard energy-conserving framework, we refer the reader to [3, 15, 18, 28, 39]. 1.6.2 An improved energy-consistent approach Many works have been devoted to extend the previous conservative properties to frictionless impact; more precisely, Laursen and Chawla [29] and Armero and Petocz [2] have shown the benefit of the persistency condition to conserve the energy in the discrete framework. Neverthelerss, in all these works the numerical method shows that the interpenetration vanishes only when the time step tends towards zero. In order to overcome this drawback, Laursen and Love [30] have developed an efficient method, by introducing a discrete jump in velocity; however, this method requires the solution of an auxiliary system in order to compute the velocity update results. Furthermore, Hauret and Le Tallec [20] have considered a specific penalized enforcement of the contact conditions which allows to provide energy conservation properties. Then, Khenous, Laborde and Renard [25] have introduced the Equivalent Mass Matrix method (EMM), based on a procedure of redistribution of the mass matrix. Interpretations and extensions of this method can be found in [19]. The resulting problem exhibits Lipschitz regularity in time and achieves good energy evolution properties, due to the fact that the persitency condition is automatically satisfied. This equivalent mass matrix approach was studied and used in many works; for instance, theoretical and computational aspects related to this model can be found in [18, 22, 25]. In what follows, we present an improved energy-conserving method for hyperelastodynamic contact problems with its extension to frictional dissipation phenomena. This method permits to enforce the normal compliance with unilateral constraint during each time step with controlled contact penetrations and with energy-consistent properties. The strategy developed is based on the solution of the system (1.36) by taking into account only the normal compliance condition with friction, in the first step, then the normal compliance restricted by unilateral constraint with friction, in the second step. This strategy is employed successively when passing from the time moment tn−1 to the time moment tn . To this end, we developed an adapted continuation Newton method, decomposed in two steps, which could be summarized as follows:


step (a): Newton scheme to solve the nonlinear system (1.36) ( −ξν n = Luν n on Γ3 , with on Γ3 . −ξ τ n ∈ −µ(uν n )ξν n ϕτ (uτ n )

15

(1.44)

step (b): Continuation of the Newton scheme to solve (1.36)  (a)  on Γ3 , if uν n < g − ξν n = Luν n (a) with if uν n ≥ g − ξν n ∈ ∂ϕν (uν n ) + Lg on Γ3 ,   −ξ τ n ∈ −µ0 ξν n ϕτ (uτ n ) on Γ3 .

(1.45)

According the work of Hauret and Le Tallec [20], we reproduce in the discrete framework the conservation properties described in (1.10) by taking into account a specific form of the normal contact reaction ξν n defined by 2

−ξν n = Luν n = rp(uν n ) with p(uν n ) =

2

[(uν n )+ ] − [(uν n−1 )+ ] . (1.46) 2(uν n − uν n−1 )

Note that p(uν n ) represents a specific form of the normal compliance function at tn . Here r is a penalization parameter interpreted as the stiffness coefficient of the asperities of the foundation. In the following, the continuation Newton method with the normal compliance form (1.46) will be called the Improved Penalized Method (IPM). To keep this paper in reasonable length, we skip the details of the solution of the nonlinear system (1.36) with conditions (1.44) and (1.45), and we restrict ourselves to recall that the presentation of the algorithms together with their numerical implementation can be found in [4, 8]. Details on the discretization step and Computational Contact Mechanics, including algorithms similar to that used here, can be found in [1, 4, 5, 24, 26, 28, 42]. 1.6.3 Analysis of the discrete energy evolution This part is devoted to establish energy-consistent properties induced by the improved penalized method described in the previous Section 1.6.2. Below we use the notation En and En−1 for the energy E of the hyperelastic frictional contact system at times tn and tn−1 , respectively. For instance, the discrete energy En can be written as follows: Z Z 1 2 ˜ (Cn ) dΓ. W (1.47) ρ [δun ] d x + En = 2 Ω Ω (a)

(b)

(a)

(b)

The notation En , En , En−1 and En−1 have similar significance, being related to the steps (a) and (b) of the numerical method introduced in Section 1.6.2. The general assessment of the discrete energy of the frictional contact Problem P¯Vhk between times tn−1 and tn is based on the following proposition.

16


Proposition 1.1. The following equality holds: En − En−1 = k hf n− 21 , un− 12 iV ∗ ×V Z ξν n δuν n + ξ τ n · δuτ n dΓ. +k

(1.48)

Γ3

Here En and En−1 denote the internal energy E of the hyperelastic frictional contact system at times tn and tn−1 , respectively. Proof. We use the variational formulation (1.36) with v = δun− 12 =

δun + δun−1 un − un−1 = . k 2

Then, we use (1.42) to get the hyperelastic energy conservation. As a consequence we obtain the equality Z Z 1 1 ρ(δun − δun−1 ).(δun + δun−1 )d x + Πalgo : ∇(un − un−1 )d x 2k Ω k Ω Z = hf n− 21 , un− 12 iV ∗ ×V + [ξν n δuν n + ξ τ n .δuτ n ]d Γ. Γ3

Furthermore, using the identity 2

2

(δun − δun−1 )(δun + δun−1 ) = [δun ] − [δun−1 ]

and the conservation property of the Gonzalez scheme given in (1.43), we obtain that Z Z 1 1 2 2 ˜ (Cn ) − W ˜ (Cn−1 )d x W ρ([δun ] − [δun−1 ] )d x + 2k Ω k Ω Z = hf n− 12 , un− 21 iV ∗ ×V + [ξν n δuν n + ξ τ n .δuτ n ]d Γ. Γ3

Finally, using the definition (1.47) of the discrete energy we obtain the identity (1.48). 2 Remark 1.2. Similar results for the discrete angular and linear momenta can also be established, see for instance [18, 20]. Based on the Proposition 1.1, we can state, at the end of the step (a), the assessment of the discrete energy for the specific normal compliance contact. Proposition 1.3. The following equality holds: (a)

En(a) − En−1 = hf n− 12 , un− 12 iV ∗ ×V Z Z r 2 2 − [(uν n )+ ] − [(uν n−1 )+ ] d Γ + k ξ τ n · δuτ n d Γ. (1.49) Γ3 2 Γ3


17

Proof. We use similar arguments as those used in the proof of the Proposition 1.1; in particular we use equality (1.46) combined with equality δuν n = uν n −uν n−1 . 2 k We remark that the form (1.46) of the normal contact reaction allows, in the frictionless case, an energy assessment which is agreement with the formula (1.10). In addition, when the external forces vanishes, the energy statements in Propositions 1.1 and 1.3 during the steps (a) and (b), respectively, allow us to obtain the following situations. Case without friction (a)

(a)

step (a): ξν n δuν n ≈ 0 ⇒ En ≈ En−1 , (b) (b) step (b): ξν n δuν n ≤ 0 ⇒ En ≤ En−1 . Case with friction (a)

(a)

step (a): ξν n δuν n ≈ 0, ξ τ n · δuτ n ≤ 0 ⇒ En ≤ En−1 , (b) (b) step (b): ξν n δuν n ≤ 0, ξ τ n · δuτ n ≤ 0 ⇒ En ≤ En−1 . In the case without friction (ξ τ n = 0), we remark that during step (a), the 2 energy of the system is almost conserved. Indeed, the difference [(uν n−1 )+ ] − 2 [(uν n )+ ] is very small since the penetrations (uν n−1 )+ and (uν n )+ are small. During the step (b), the enforcement of the unilateral constraint allows to limit the penetrations obtained in step (a) by the value g, which represents a small given value. On the other hand, we can easily prove that the product ξν n δuν n is always negative and this represents an unnacceptable physical behavior, since it generates energy dissipation. However, this energy dissipation is low because the impact has occurred during the step (a). Furthermore, when the friction case is considered, in both steps (a) and (b) we observe an admissible dissipation of energy. Indeed, in this case the inner product ξ τ n · δuτ n is always negative. In other words, this strategy limits the energy dissipation between times tn and tn−1 , in the frictionless case, and allows the energy dissipation, in the frictional one. To resume, the advantages of the method arise in the fact that both the dissipation of energy and the penetrations are limited during the impact.

1.7 Numerical experiments In order to recover the theoretical numerical behaviour of the fully discrete scheme discussed in Section 1.6.3, we carried out some numerical simulations based on two representative dynamic contact problems: the impact without friction of a linear elastic ball against a foundation (Subsection 1.7.1) and the impact with friction of a hyperelastic ring against a foundation (Subsection 1.7.2).

18


1.7.1 Impact of a linear elastic ball against a foundation This representative benchmark problem describes the frictionless impact of an linear elastic ball against a foundation (see [24]). The elastic ball is thrown with an initial velocity (u1 = (0, −10)m/s) toward the foundation (x1 , x2 ) ∈ R2 : x2 ≤ 0 . The material’s behavior is described by an elastic linear constitutive law defined by the energy function W (ǫ) =

Eκ E (trǫ)2 + tr(ǫ2 ), 2(1 + κ)(1 − 2κ) 2(1 + κ)

∀ǫ ∈ Mn .

Here E and κ are Young’s modulus and Poisson’s ratio of the material and tr(·) denotes the trace function, respectively. Note that ǫ = 12 (∇uT + ∇u) represents the linearized strain tensor in the framework of the of small deformations theory (kuk