A quasistatic contact problem with normal compliance and damage involving viscoelastic materials with long memory

Applied Numerical Mathematics 58 (2008) 1274–1290 www.elsevier.com/locate/apnum

A quasistatic contact problem with normal compliance and damage involving viscoelastic materials with long memory M. Campo a , J.R. Fernández a,∗ , Á. Rodríguez-Arós b a Departamento de Matemática Aplicada, Universidade de Santiago de Compostela Facultade de Matemáticas, Campus Sur s/n,

15782 Santiago de Compostela, Spain b Departamento de Métodos Matemáticos e de Representación, Universidade de A Coruña, E.T.S.E. Camiños, Canais e Portos,

Campus de Elviña s/n, 15071 A Coruña, Spain Available online 24 July 2007

Abstract In this work, a model for the quasistatic frictionless contact between a viscoelastic body with long memory and a foundation is studied. The material constitutive relation is assumed to be nonlinear and the contact is modelled with the normal compliance condition, i.e., the obstacle is assumed deformable. The mechanical damage of the material, caused by excessive stress or strain, is described by the damage function, which is modelled by a nonlinear partial differential equation. We derive a variational formulation for the problem and prove the existence of its unique weak solution. Then, we introduce a fully discrete scheme for the numerical solutions of the problem, based on the finite element method to approximate the spatial variable and an Euler scheme to discretize the time derivatives, and we obtain error estimates on the approximate solutions. Finally, some numerical results are presented in the simulation of two-dimensional test problems. © 2007 IMACS. Published by Elsevier B.V. All rights reserved. MSC: 49J40; 65M15; 74M15; 74S05 Keywords: Viscoelastic material with long memory; Damage; Frictionless contact; Weak solutions; Finite elements; Error estimates; Numerical simulations

1. Introduction In this paper, a contact problem between a viscoelastic body with long memory and a deformable obstacle is studied. The effect due to the damage of the material, caused by the opening and growth of microcracks which lead to the decrease in the load carrying capacity of the body, is also included. The effective functioning and safety of a mechanical system may be deteriorated by this decrease as the material undergoes damaged. Because of the importance of this topic, an increasing number of engineering publications dealing with damage models appeared in the last ten years (see, for instance, [2,26,28,29]). Here, following the early works of Frémond (see, e.g., [17,18], and also the monograph [16] for details), the damage is modelled using an internal variable ζ taking values between 0 and 1, in such a way that when ζ = 1 the material is damage-free (that is, a material with long memory), when ζ = 0 the material is completely damaged, and when * Corresponding author. Tel.: +34981563100x 23244; fax: +34982285926.

E-mail addresses: [email protected] (J.R. Fernández), [email protected] (Á. Rodríguez-Arós). 0168-9274/$30.00 © 2007 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2007.07.003

M. Campo et al. / Applied Numerical Mathematics 58 (2008) 1274–1290

1275

Fig. 1. Contact problem of a viscoelastic body with long memory.

0 < ζ < 1 there is partial damage. Mathematical problems including damage have been studied recently for elastic materials [4,5,9], elastic-viscoplastic materials [7,10,11] or viscoelastic materials [6,20]. Moreover, one-dimensional problems have been also considered (see, e.g., [1]). Contact problems involving viscoelastic materials with long memory have been studied in [15,31,33]. This paper extends the results provided in [31] to the case with damage and nonlinear constitutive functions, including a numerical algorithm to solve the fully discrete problem and performing some numerical simulations, and it is parallel to [8,9]. Here, a frictionless contact problem involving a viscoelastic material with long memory and damage is considered. The obstacle is assumed deformable and thus, a normal compliance contact condition is employed (see, for example, [22]). We perform the variational analysis of this problem, including the existence of a unique solution to the variational problem. Our interest concerns mainly the numerical analysis of the model including numerical simulations. An algorithm is developed, by using a penalty-duality algorithm, for the numerical resolution of the fully discrete problem. The paper is organized as follows. In Section 2 we state the mechanical problem, list the assumptions on the data and derive a variational formulation of the model. An existence and uniqueness result is proved, based on results of parabolic variational equations and Banach’s fixed point theorem. In Section 3, we analyze a fully discrete scheme for the problem. We use the finite element method to approximate the domain and a hybrid combination of the forward and the backward Euler schemes to discretize the time derivatives, and derive related error estimates. Finally, in Section 4 we present numerical simulations in the study of two-dimensional test problems. 2. Problem statement and variational formulation We consider a viscoelastic body which occupies a domain Ω ⊂ Rd for d = 1, 2, 3. The boundary Γ of the domain Ω is assumed to be Lipschitz, and it is divided into three disjoint measurable parts ΓD , ΓF and ΓC such that meas (ΓD ) > 0. The body is clamped on ΓD × (0, T ) and so the displacement field vanishes there; a volume force of density f 0 acts in Ω × (0, T ), and surface tractions of density f 2 act on ΓF × (0, T ). Here, T > 0 and [0, T ] is the time interval of interest. The body may come in contact with a reactive foundation over the part ΓC . A gap g exists between the potential contact surface ΓC and the foundation, and it is measured along the outward normal vector ν (see Fig. 1). We denote by u the displacement field, σ the stress tensor and ε(u) the linearized strain tensor. Moreover, let ζ be the damage field which measures the density of the microcracks in the material, and will be described below. The material is assumed viscoelastic with long memory, and it obeys the following constitutive law,   σ (t) = Aε u(t) +

t

    G t − s, ε u(s) , ζ (s) ds,

0

where A represents the elasticity tensor and G is a prescribed nonlinear constitutive function. This law is a generalization of that used for example in [14,25,30,34] to model the mechanical behavior of a variety of solid materials, as organic polymers or some kinds of wood and it has been used in the context of contact mechanics in, e.g., [15,31,33]. To simplify the notation, we do not indicate explicitly the dependence of various functions on the independent variables x ∈ Ω ∪ Γ and t ∈ [0, T ]. Also, dots above a variable represent time derivatives.

1276

M. Campo et al. / Applied Numerical Mathematics 58 (2008) 1274–1290

We assume that the normal stress on ΓC , σν = σ ν · ν, is given by the following normal compliance contact condition (see [22]), −σν = p(uν − g). Here, uν = u · ν denotes the normal displacement. When uν > g, the difference uν − g represents the interpenetration of the body’s asperities into those of the foundation. The normal compliance function p is prescribed and satisfies p(r) = 0 for r  0, since then there is no contact. Moreover, we assume that the contact is frictionless, i.e. the tangential component of the stress σ τ = σ − σν ν vanishes on the contact surface. We turn to describe the damage process. As a result of the tensile or compressive stresses in the body, microcracks open and grow and this, in turn, causes the load bearing of the material to decrease. This reduction in the strength of the material is modelled by introducing the damage field ζ = ζ (x, t). Following the derivation in [17,18], the evolution of the microscopic cracks responsible for the damage is described by the partial differential equation (see [23,24]),   ζ˙ − κ ζ = φ ε(u), ζ . Here, κ > 0 is a constant and φ is a given constitutive function which describes damage sources in the system. Remark 1. We notice that here, according to [23,24], the damage is assumed to satisfy a partial differential equation. In [10,11,20] (see the monograph [32] for details), this is replaced by a parabolic subdifferential inclusion because a constraint condition must be imposed. In this paper, this condition is obtained from the assumptions required to the damage source function φ. Following [20], an homogeneous Neumann boundary condition for the damage field is used, so the normal derivative of ζ , denoted by ∂∂ζν , vanishes on Γ . Let Sd be the space of second order symmetric tensors on Rd , or equivalently, the space of symmetric matrices of order d, and denote by · the usual inner product defined in these spaces and by  ·  their norm. The classical form of the mechanical problem of the quasistatic frictionless contact with damage of a viscoelastic body with long memory may be written as follows. Problem P . Find a displacement field u : Ω × [0, T ] → Rd , a stress field σ : Ω × [0, T ] → Sd , and a damage field ζ : Ω × [0, T ] → R such that ζ (0) = ζ0 in Ω and, − Div σ = f 0

in Ω × (0, T ), t       σ (t) = Aε u(t) + G t − s, ε u(s) , ζ (s) ds

(1) in Ω, a.e. t ∈ (0, T ),

(2)

0

  ζ˙ − κ ζ = φ ε(u), ζ

in Ω × (0, T ),

(3)

∂ζ = 0 on Γ × (0, T ), ∂ν u = 0 on ΓD × (0, T ),

(5)

σν = f 2

on ΓF × (0, T ),

(6)

σ τ = 0,

−σν = p(uν − g)

(4)

on ΓC × (0, T ).

(7)

We now proceed to obtain a variational formulation of Problem P . For this purpose, we introduce additional notation and assumptions on the problem data. Here and in what follows the indices i and j run between 1 and d, the summation convention over repeated indices is adopted and the index that follows a comma indicates a partial derivative with respect to the corresponding component of the independent variable. Let Y = L2 (Ω), E = H 1 (Ω), H = [L2 (Ω)]d and define the following variational spaces:

M. Campo et al. / Applied Numerical Mathematics 58 (2008) 1274–1290

1277

   d V = v ∈ H 1 (Ω) ; v = 0 on ΓD ,   Q = σ = (σij ); σij = σj i ∈ L2 (Ω) . For every element v ∈ V , we denote by vν and v τ the normal and the tangential components of v on the boundary Γ given by vν = v · ν and v τ = v − vν ν, respectively. Similarly, for a regular (say C 1 ) tensor field σ : Ω → Sd we define its normal and tangential components by σν = (σ ν) · ν and σ τ = σ ν − σν ν. On V we consider the inner product given by (u, v)V = (ε(u), ε(v))Q , and the associated norm vV = ε(v)Q . Since meas (ΓD ) > 0, from Korn’s inequality it follows that (V ,  · V ) is a real Hilbert space. Moreover, by the Sobolev’s trace theorem we have a constant c0 , depending only on the domain Ω, ΓD and ΓC , such that v[L2 (ΓC )]d  c0 vV ∀v ∈ V . For every real Hilbert space X, let  · X be the associated norm on X, we use the classical notation for the spaces Lp (0, T ; X) and W k,p (0, T ; X), 1  p  +∞, k  1, and we denote by C([0, T ]; X) and C 1 ([0, T ]; X) the spaces of continuous and continuously differentiable functions from [0, T ] to X, respectively. In the study of the mechanical problem (1)–(7), we assume that the elasticity tensor A(x) = (aij kl (x))di,j,k,l=1 : τ ∈ Sd → A(x)τ ∈ Sd satisfies ⎫ (a) aij kl = aklij = aj ikl for i, j, k, l = 1, . . . , d. ⎪ ⎪ ⎬ (b) aij kl ∈ L∞ (Ω) for i, j, k, l = 1, . . . , d. (8) 2 (c) There exists mA > 0 such that A(x)τ · τ  mA τ  ⎪ ⎪ ⎭ d ∀τ ∈ S , a.e. x ∈ Ω. The viscoelastic long memory operator G : Ω × [0, T ] × Sd × R → Sd satisfies ⎫ (a) There exists MG > 0 such that  ⎪ ⎪ ⎪

G(x, t, ξ 1 , β1 ) − G(x, t, ξ 2 , β2 )  MG ξ 1 − ξ 2  + |β1 − β2 | ⎪ ⎪ ⎪ ⎪ d ⎪ ∀t ∈ [0, T ], ξ 1 , ξ 2 ∈ S , β1 , β2 ∈ R, a.e. x ∈ Ω. ⎪ ⎪ ⎬ (b) The mapping x → G(x, t, ξ , β) is Lebesgue measurable in Ω, d ∀t ∈ [0, T ], ξ ∈ S , β ∈ R. ⎪ ⎪ ⎪ ⎪ (c) The mapping t → G(x, t, ξ , β) is continuous in [0, T ], ⎪ ⎪ ⎪ d ⎪ ∀ξ ∈ S , β ∈ R, a.e. x ∈ Ω. ⎪ ⎪ ⎭ (d) The mapping x → G(x, t, 0, 0) ∈ Q ∀t ∈ [0, T ].

(9)

The damage source function φ : Ω × Sd × R → R satisfies

⎫ (a) There exists Lφ > 0 such that    ⎪ ⎪ φ(x, ε 1 , β1 ) − φ(x, ε2 , β2 )  Lφ ε1 − ε 2  + |β1 − β2 | ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ∀ε1 , ε 2 ∈ Sd , β1 , β2 ∈ R, a.e. x ∈ Ω. ⎪ ⎪ ⎬ (b) The mapping x → φ(x, ε, β) is Lebesgue measurable on Ω, ∀ε ∈ Sd , β ∈ R. ⎪ ⎪ ⎪ ⎪ (c) The mapping x → φ(x, 0, 0) ∈ L2 (Ω). ⎪ ⎪ ⎪ d ⎪ ⎪ (d) φ(x, ε, β) is bounded ∀ε ∈ S , β ∈ R, a.e. x ∈ Ω. ⎪ ⎭ (e) φ(x, ε, β)  0 if β  1, φ(x, ε, β)  0 if β  ζ∗ .

(10)

Remark 2. We notice that condition (10)(e) allows us to prove that the damage field ζ satisfies the constraint ζ (x)  ζ∗ , ∀x ∈ Ω, where ζ∗ > 0 is assumed small enough. We restrict the damage to ζ∗  ζ since, when the effective modulus is small, the material is full of microcracks and it no longer makes sense to model it by a viscoelastic material law. We suppose that the body forces and surface tractions are given as,     d  f 0 ∈ C [0, T ]; H , f 2 ∈ C [0, T ]; L2 (ΓF ) ,

(11)

and the initial damage field satisfies ζ0 ∈ E,

1  ζ0 (x)  ζ∗

a.e. x ∈ Ω.

(12)

1278

M. Campo et al. / Applied Numerical Mathematics 58 (2008) 1274–1290

The normal compliance function p : ΓC × R → R+ verifies

⎫ (a) There exists Lp > 0 such that ⎪ ⎪ ⎪ p(x, r1 ) − p(x, r2 )  Lp |r1 − r2 | ∀r1 , r2 ∈ R, a.e. x ∈ ΓC . ⎪ ⎬ , ∀r ∈ R. (b) The mapping x → p(x, r) is Lebesgue measurable on Γ C  ⎪ ⎪ (c) p(x, r1 ) − p(x, r2 ) · (r1 − r2 )  0 ∀r1 , r2 ∈ R, a.e. x ∈ ΓC . ⎪ ⎪ ⎭ (d) The mapping x → p(x, r) = 0 for all r  0.

(13)

Next, we define the function f : [0, T ] → V , the functional j : V × V → R and the bilinear form a : E × E → R by the following equations:       f (t), v V = f 0 (t), v H + f 2 (t), v [L2 (Γ )]d ∀v ∈ V , t ∈ [0, T ], F  (14) j (u, v) = p(uν − g)vν da ∀u, v ∈ V , ΓC

 ∇ξ · ∇ψ dx

a(ξ, ψ) = κ

∀ξ, ψ ∈ E.

Ω

Notice that the integral (14) is well defined by (13), and conditions (11) imply the regularity   f ∈ C [0, T ]; V . We turn now to derive a variational formulation of the mechanical problem P . To that end, we assume that (u, σ , ζ ) are regular functions satisfying (1)–(7). Applying a Green’s formula and using the boundary conditions (4)–(7), we obtain the following variational formulation of the mechanical problem P . Problem VP. Find a displacement field u : [0, T ] → V , a stress field σ : [0, T ] → Q, and a damage field ζ : [0, T ] → E such that ζ (0) = ζ0 , and for a.e. t ∈ [0, T ],   σ (t) = Aε u(t) +

t

    G t − s, ε u(s) , ζ (s) ds,

(15)

0

      σ (t), ε(w) Q + j u(t), w = f (t), w V ∀w ∈ V ,           ζ˙ (t), ξ L2 (Ω) + a ζ (t), ξ = φ ε u(t) , ζ (t) , ξ L2 (Ω)

(16) ∀ξ ∈ E.

(17)

The existence of a unique solution to Problem VP is proved using some results for parabolic variational equations and Banach fixed point arguments, which we summarize in the following. Theorem 3. Assume that (8)–(13) hold. Then, there exists a unique solution (u, σ , ζ ) to Problem VP with the following regularity:     u ∈ C [0, T ]; V , σ ∈ C [0, T ]; Q , ζ ∈ W 1,2 (0, T ; Y ) ∩ L2 (0, T ; E). Moreover, the damage field ζ satisfies ζ (t, x) ∈ [ζ∗ , 1] a.e. t ∈ [0, T ] and x ∈ Ω. Proof. As a first step towards an existence and uniqueness result, we introduce the operator A : V → V defined as follows,   (Av, w)V = Aε(v), ε(w) Q + j (v, w) ∀v, w ∈ V . From (8) and (13), we claim that A is a strongly monotone and Lipschitz continuous operator. Now, let θ ∈ C([0, T ]; Y ) be a given auxiliary function and consider the following variational problem. Problem VPθ . Find uθ : [0, T ] → V and ζθ : [0, T ] → E such that ζθ (0) = ζ0 , and for a.e. t ∈ [0, T ],

M. Campo et al. / Applied Numerical Mathematics 58 (2008) 1274–1290

  Auθ (t), w V +

 t

       G t − s, ε uθ (s) , ζθ (s) ds, ε(w) = f (t), w V

  + a ζθ (t), ξ = (θ, ξ )L2 (Ω) L2 (Ω)

(18)

Q

0

  ζ˙θ (t), ξ

∀w ∈ V ,

1279

∀ξ ∈ E.

(19)

We observe that the expression (19) is uncoupled. By using some arguments of evolutionary variational equations (see, e.g., [27]), it follows that there exists a unique solution to (19) satisfying ζθ (0) = ζ0 and the regularity, ζθ ∈ W 1,2 (0, T ; Y ) ∩ L2 (0, T ; E). Now, let μθ ∈ C([0, T ]; Q) be a given auxiliary function. There exists a unique η θ ∈ C([0, T ]; V ) such that     ηθ (t), w V = μθ (t), ε(w) Q ∀w ∈ V , t ∈ [0, T ].

(20)

We consider the following variational problem. Problem VPθ,η . Find uθ,η : [0, T ] → V such that       Auθ,η (t), v V + ηθ (t), v V = f (t), v V

∀v ∈ V , t ∈ [0, T ].

(21)

By using standard arguments of variational equations, we deduce that there exists a unique solution to Problem VPθ,η satisfying   uθ,η ∈ C [0, T ]; V . We introduce the operator Λθ : C([0, T ]; Q) → C([0, T ]; Q) given by t Λθ (μθ )(t) =

    G t − s, ε uθ,η (s) , ζθ (s) ds,

∀t ∈ [0, T ],

0

where η θ is defined as in (20), and uθ,η and ζθ are the solutions to problems VPθ,η and (19) with ζθ (0) = ζ0 , respectively. Now, given μ1 , μ2 ∈ C([0, T ]; Q), we have

Λθ (μ1 )(t) − Λθ (μ2 )(t)  Q

t

       

G t − s, ε uθ,η (s) , ζθ (s) − G t − s, ε uθ,η (s) , ζθ (s) ds 1 2 Q

0

t  MG



uθ,η (s) − uθ,η (s) ds 1 2 V

∀t ∈ [0, T ].

(22)

0

By taking in (21) η θ = η 1 and v = uθ,η2 − uθ,η1 and then η θ = η 2 and v = uθ,η1 − uθ,η2 and doing some algebra, we find that





uθ,η (s) − uθ,η (s)  c η1 (s) − η2 (s) = c μ1 (s) − μ2 (s) . (23) 1 2 V V Q Using (23) in (22), we obtain

Λθ (μ1 )(t) − Λθ (μ2 )(t)  c Q

t



μ1 (s) − μ2 (s) ds Q

∀t ∈ [0, T ],

0

where here and in what follows, c is a positive constant whose value may change from line to line. By reiterating the previous expression we find that there exists a number p ∈ N such that p p

p

Λ (μ1 )(t) − Λp (μ2 )(t)  c T μ1 − μ2 C ([0,T ];Q) θ θ Q p!

∀t ∈ [0, T ],

1280

M. Campo et al. / Applied Numerical Mathematics 58 (2008) 1274–1290

and cp T p /p! < 1. Thus, applying the Banach fixed point theorem, it follows that there exists a unique μ∗ ∈ C([0, T ]; Q) such that Λθ (μ∗ )(t) = μ∗ (t) for all t ∈ [0, T ]. As a consequence, uθ = uθ,η∗ ∈ C([0, T ]; V ) is the unique solution of (18) and the pair (uθ , ζθ ) is the unique solution to Problem VPθ . We introduce the operator Θ : C([0, T ]; Y ) → C([0, T ]; Y ) defined by     Θ(θ )(t) = φ ε uθ (t) , ζθ (t) ∀t ∈ [0, T ], where the pair (uθ , ζθ ) is the solution to Problem VPθ . Now, given θ1 , θ2 ∈ C([0, T ]; Y ), we find that

       

Θ(θ1 )(t) − Θ(θ2 )(t) 2 = φ ε uθ (t) , ζθ (t) − φ ε uθ (t) , ζθ (t) 2 1 1 2 2 Y Y

2

2   ∀t ∈ [0, T ].  L2φ uθ1 (t) − uθ2 (t) V + ζθ1 (t) − ζθ2 (t) Y

(24)

On the other hand, by taking in (19) θ = θ1 and ξ = ζθ2 (t) − ζθ1 (t) and then θ = θ2 and ξ = ζθ1 (t) − ζθ2 (t), it follows that

    1 d

ζθ (t) − ζθ (t) 2 + a ζθ (t) − ζθ (t), ζθ (t) − ζθ (t) = θ1 (t) − θ2 (t), ζθ (t) − ζθ (t) . 1 2 1 2 1 2 1 2 Y Y 2 dt By integrating in [0, t], we obtain

1

ζθ (t) − ζθ (t) 2 + 1 2 Y 2

t

  a ζθ1 (s) − ζθ2 (s), ζθ1 (s) − ζθ2 (s) ds =

0

t

  θ1 (s) − θ2 (s), ζθ1 (s) − ζθ2 (s) Y ds.

0

Thus,

1

ζθ (t) − ζθ (t) 2  1 2 Y 2

t



ζθ (s) − ζθ (s) 2 ds + 1 2 Y

0

t



θ1 (s) − θ2 (s) 2 ds. Y

0

By using the Gronwall’s Lemma, we find that

ζθ (t) − ζθ (t) 2  c 1 2 Y

t



θ1 (s) − θ2 (s) 2 ds. Y

(25)

0

Also, by a similar argument to that of (22), it follows that

uθ (t) − uθ (t) 2  c 1 2 V

t



ζθ (s) − ζθ (s) 2 ds. 1 2 Y

(26)

0

Summarizing, using (25) and (26) in (24), it leads to the following inequality,

Θ(θ1 )(t) − Θ(θ2 )(t) 2  c Y

t



θ1 (s) − θ2 (s) 2 ds. Y

(27)

0

Reiterating the previous expression and proceeding as we did with the operator Λθ , we find that there exists a θ ∗ ∈ C([0, T ]; Y ) which is the unique fixed point of Θ. We define u = uθ ∗ , ζ = ζθ ∗ and σ by (15). Then, (u, σ , ζ ) is a solution to Problem VP. The uniqueness is a consequence of the fixed point theorem. Finally, we need to prove that ζ ∈ [ζ∗ , 1]. This is based on the following comparison theorem (see [23] for details). Theorem 4. Let the damage source function φ and the initial damage field ζ0 satisfy assumptions (10) and (12), respectively. Then, the solution, ζ , to           ζ˙ (t), ξ L2 (Ω) + a ζ (t), ξ = φ ε u(t) , ζ (t) , ξ L2 (Ω) ∀ξ ∈ E, ζ (0) = ζ0 , satisfies ζ (t, x) ∈ [ζ∗ , 1] a.e. t ∈ [0, T ] and x ∈ Ω. This concludes the proof of Theorem 3.

M. Campo et al. / Applied Numerical Mathematics 58 (2008) 1274–1290

1281

3. Fully discrete approximation: error estimates In this section we consider a numerical scheme for solving Problem VP and derive some error estimates on the approximate solutions. We use two finite dimensional spaces V h ⊂ V and E h ⊂ E to approximate the spaces V and E, respectively. Here, h > 0 denotes the spatial discretization parameter. To discretize time derivatives, let us consider a uniform partition of the time interval [0, T ], denoted by 0 = t0 < t1 < · · · < tN = T , and let k = T /N be the time step size. For a continuous function f (t), we use the notation fn = f (tn ). Moreover, for a sequence {wn }N n=0 , let δwn = (wn − wn−1 )/k denote its divided differences, and let c be a positive constant independent of the discretization parameters h and k. Remark 5. In the numerical simulations presented in the next section, V h and E h are composed of continuous and piecewise affine functions, that is,   d  d  V h = v h ∈ C(Ω) ; v h |T ∈ P1 (T ) ∀T ∈ T h , v h = 0 on ΓD ,   E h = ξ h ∈ C(Ω); ξ h |T ∈ P1 (T ) ∀T ∈ T h ,

(28) (29)

where Ω is assumed to be a polyhedral domain and we denote by T h a regular family of triangulations of Ω compatible with the partition of the boundary Γ = ∂Ω into ΓD , ΓF and ΓC . Let ζ0h be an appropriate approximation of the initial condition ζ0 . A fully discrete approximation of Problem VP, based on a hybrid combination of the forward and the backward Euler schemes, is the following. N h and a discrete damage field ζ hk = Problem VPhk . Find a discrete displacement field uhk = {uhk n }n=0 ⊂ V N hk h hk h {ζn }n=0 ⊂ E such that ζ0 = ζ0 and, for n = 1, . . . , N ,

 hk h       hk  h   δζn , ξ Y + a ζnhk , ξ h = φ ε uhk ∀ξ h ∈ E h , n−1 , ζn−1 , ξ Y   hk   h   hk h  Aε un , ε w Q + j un , w   n      hk  hk   h  h = f n, w V − kG tn − tj −1 , ε uj −1 , ζj −1 , ε w ∀wh ∈ V h , j =1

(30)

(31)

Q

h where uhk 0 ∈ V is given.

Using classical arguments of nonlinear variational equations (see [19]), we obtain the existence of a unique solution to Problem VPhk . Theorem 6. Assume that (8)–(13) hold. Then there exists a unique solution (uhk , ζ hk ) ⊂ V h × E h to Problem VPhk . Remark 7. We notice that Problem VPhk is solved at each time step as follows. First, the discrete “initial condition” uhk 0 is the solution of the problem,   hk   h      h h h h Aε u0 , ε w Q + j uhk 0 , w = f 0 , w V , ∀w ∈ V . This is a nonlinear variational equation which is solved by using a penalty-duality algorithm (see [3]). Next, for each n ∈ {1, . . . , N}, the discrete damage field is obtained solving the discrete variational equation (30) by employing Cholesky’s method since it leads to a linear system. Moreover, the discrete displacement field is the solution of the nonlinear variational equation (31). Again, the penalty-duality algorithm introduced in [3] is employed. We assume the following additional regularity of the damage field,     ζ ∈ C 1 [0, T ]; Y ∩ C [0, T ]; E .

(32)

1282

M. Campo et al. / Applied Numerical Mathematics 58 (2008) 1274–1290

hk The aim of this section is to estimate the numerical errors un − uhk n and ζn − ζn . First, let us obtain an error estimate hk = φ(ε(uhk ), ζ hk ). Writing Eq. (17) at time t = t for for the damage field. Denote by φn = φ(ε(un ), ζn ) and φn−1 n n−1 n−1 all ξ = ξ h ∈ E h and subtracting it to Eq. (30), it follows that       hk , ξ h Y = 0, ∀ξ h ∈ E h , ζ˙n − δζnhk , ξ h Y + a ζn − ζnhk , ξ h − φn − φn−1

and therefore,       hk ζ˙n − δζnhk , ζn − ζnhk Y + a ζn − ζnhk , ζn − ζnhk − φn − φn−1 , ζn − ζnhk Y       hk = ζ˙n − δζnhk , ζn − ξ h Y + a ζn − ζnhk , ζn − ξ h − φn − φn−1 , ζn − ξ h Y , for all ξ h ∈ E h . Then, after some algebra we obtain (see [9]),

    2 δζn − δζnhk , ζn − ζnhk Y + c ∇ ζn − ζnhk H



 hk hk

 c un − uhk n−1 V + ζn − ζn−1 Y ζn − ζn Y



 

 + ζ˙n − δζn Y ζn − ζnhk Y + ζn − ξ h Y + ∇ ζn − ζnhk H ζn − ξ h E



   hk h hk h

+ c un − uhk n−1 V + ζn − ζn−1 Y ζn − ξ Y + δζn − δζn , ζn − ξ Y , where δζn = (ζn − ζn−1 )/k. Since

   1 

ζn − ζ hk 2 − ζn−1 − ζ hk 2 , δζn − δζnhk , ζn − ζnhk Y  n Y n−1 Y 2k using the Cauchy’s inequality 1 2 b , a, b,  ∈ R,  > 0, 4 several times, it leads to the following estimate,









 

ζn − ζ hk 2 + k ∇ ζn − ζ hk 2  ck un − uhk 2 + ζn − ζ hk 2 + ζ˙n − δζn 2 + ζn − ξ h 2 n Y n n−1 V n−1 Y H Y E

    hk hk h hk 2

+ c ζn − ζn − ζn−1 − ζn−1 , ζn − ξ Y + ζn−1 − ζn−1 Y , ab  a 2 +

(33)

for all ξ h ∈ E h . By induction we obtain n n  



 

uj −1 − uhk 2 + ζj −1 − ζ hk 2

∇ ζj − ζ hk 2  ck

ζn − ζ hk 2 + k n Y j j −1 V j −1 Y H j =1

j =1

2  + ζ˙j − δζj 2Y + uj − uj −1 2V + ζj − ζj −1 2Y + ζj − ξjh E +c

n 

2     ζj − ζjhk − ζj −1 − ζjhk−1 , ζj − ξjh Y + ζ0 − ζ0h Y ,

(34)

j =1

for all ξ h = {ξjh }nj=0 ⊂ E h . Next, let us estimate the numerical errors for the displacement field. Plugging equation(15) at time t = tn into Eq. (16) and subtracting it to Eq. (31), it follows that     h      h Aε(un ) − Aε uhk + j un , w h − j uhk n ,ε w n ,w Q  tn n        hk  hk   h  + G tn − s, ε u(s) , ζ (s) ds − kG tn − tj −1 , ε uj −1 , ζj −1 , ε w = 0, 0

for all w h ∈ V h . Therefore, we have

j =1

Q

M. Campo et al. / Applied Numerical Mathematics 58 (2008) 1274–1290

1283

    hk      hk hk Aε(un ) − Aε uhk + j un , un − uhk n , ε un − un n − j un , un − un Q tn  n        hk  hk    hk + G tn − s, ε u(s) , ζ (s) ds − kG tn − tj −1 , ε uj −1 , ζj −1 , ε un − un 0

j =1

Q

0

j =1

Q

         h h = Aε(un ) − Aε uhk + j un , un − wh − j uhk n , ε un − w n , un − w Q  tn n         hk  hk   h G tn − s, ε u(s) , ζ (s) ds − kG tn − tj −1 , ε uj −1 , ζj −1 , ε un − w . + From the property (13), it is easy to check that   hk   hk j un , un − uhk n − j un , un − un  0,



    h h hk

j un , un − wh − j uhk n , un − w  c un − w V un − un V , and using the properties (8), (9) and applying repeatedly the Cauchy’s inequality (33), we find that  n   

2

2

2 

2 hk hk hk 2 h

un − u  c k uj −1 − u + ζj −1 − ζ + I + un − w , ∀wh ∈ V h , n

V

j −1 V

j =1

j −1 Y

G ,n

V

(35)

where IG ,n is the integration error given by

tn n

      

IG ,n = G tn − s, ε u(s) , ζ (s) ds − kG tn − tj −1 , ε(uj −1 ), ζj −1 .

j =1

0

Q

Taking into account that n      ζj − ζjhk − ζj −1 − ζjhk−1 , ζj − ξjh Y j =1

n−1          = ζn − ζnhk , ζn − ξnh Y − ζ0 − ζ0h , ζ1 − ξ1h Y + ζj − ζjhk , ζj − ξjh − ζj +1 − ξjh+1 Y j =1



2

2

2

2   ζn − ζnhk Y + c ζn − ξnh Y + c ζ0 − ζ0h Y + c ζ1 − ξ1h Y n−1 n−1 

  1  hk 2

ζj − ξ h − ζj +1 − ξ h 2 , ζj − ζj Y + +k j j +1 Y k j =1

j =1

where  > 0 is assumed sufficiently small, combining Eqs. (34) and (35), it leads to the following estimate for all w h ∈ V h and {ξjh }nj=1 ⊂ E h , n 







∇ ζj − ζ hk 2

un − uhk 2 + ζn − ζ hk 2 + k n V n Y j H j =1

 ck

n  j =1





uj −1 − uhk 2 + ζj −1 − ζ hk 2 + ζ˙j − δζj 2 + uj − uj −1 2 + ζj − ζj −1 2 Y V Y j −1 V j −1 Y



2 

2

2 + ζj − ξjh E + cIG2 ,n + cun − w h 2V + c ζn − ξnh Y + c ζ1 − ξ1h Y

  1 

ζj − ξ h − ζj +1 − ξ h 2 + c ζ0 − ζ h 2 . +c j j +1 Y 0 Y k n−1 j =1

Applying a discrete version of Gronwall’s lemma (see [21]), we obtain the following main error estimates result.

1284

M. Campo et al. / Applied Numerical Mathematics 58 (2008) 1274–1290

Theorem 8. Assume that (8)–(13) and the additional regularity (32) hold. Let (u, ζ ) ∈ V × E and (uhk , ζ hk ) ⊂ V h × E h be the respective solutions to Problems VP and VPhk . Then, we have the following error estimates for all h N h h w h = {w hj }N j =1 ⊂ V and {ξj }j =1 ⊂ E , max

0nN

N 



  

un − uhk 2 + ζn − ζ hk 2 + k

∇ ζj − ζ hk 2 n V n Y j H

 ck

j =1

N  j =1

2   ζ˙j − δζj 2Y + uj − uj −1 2V + ζj − ζj −1 2Y + ζj − ξjh E

2

2 + c max IG2 ,n + c max un − w hn V + c max ζn − ξnh Y 0nN 0nN 0nN N −1



  1 

ζj − ξ h − ζj +1 − ξ h 2 + c ζ0 − ζ h 2 + c u0 − uhk 2 , +c j j +1 Y 0 Y 0 V k

(36)

j =1

h where uhk 0 ∈ V is given.

Error estimates (36) are the basis for the analysis of the convergence rate of the algorithm. As an example, let Ω be a polyhedral domain and denote by T h a regular triangulation of Ω compatible with the partition of the boundary Γ = ∂Ω into ΓD , ΓF and ΓC . Let V h and E h be defined by (28) and (29), respectively, and assume that the discrete initial condition ζ0h is obtained by ζ0h = π h ζ0 , where π h : C(Ω) → E h is the standard finite element interpolation operator (see, e.g., [13]). Moreover, the “discrete hk h d h h d h initial condition” uhk 0 is defined as u0 = Π u(0) where Π = (πi )i=1 : [C(Ω)] → V . The following corollary is then derived, establishing the linear convergence of the algorithm with respect to the discretization parameters h and k. Corollary 9. Let the assumptions of Theorem 8 hold. Under the following regularity conditions d    u ∈ C [0, T ]; H 2 (Ω) ∩ W 1,∞ (0, T ; V ),   ζ ∈ H 2 (0, T ; Y ) ∩ C [0, T ]; H 2 (Ω) ∩ H 1 (0, T ; E), the numerical algorithm introduced in Problem VPhk is linearly convergent; that is, there exists c > 0, independent of h and k, such that,

  hk

max un − uhk (37) n V + ζn − ζn Y  c(h + k). 0nN

Proof. First, we need to estimate the errors provided by the approximation of the finite element spaces V h and E h . Since u ∈ C([0, T ]; [H 2 (Ω)]d ) and ζ ∈ C([0, T ]; H 2 (Ω)), we have (see [13]),

max inf un − w hn V  chuC ([0,T ];[H 2 (Ω)]d ) , 0nN wnh ∈V h

max

inf ζn − ξnh E  chζ C ([0,T ];H 2 (Ω)) .

0nN ξnh ∈E h

From the definition of the operators π h and Π h we obtain,



u0 − uhk  ch u(0) 2 , 0 V [H (Ω)]d

ζ0 − ζ h  ch2 ζ0  2 . H (Ω) 0 Y  2 2 2 ˙ Since ζ ∈ H 2 (0, T ; Y ), it follows that N j =1 kζj − δζj Y  ck ζ H 2 (0,T ;Y ) , and proceeding as in [12] we find that

M. Campo et al. / Applied Numerical Mathematics 58 (2008) 1274–1290

1285

  IG ,n  ck uW 1,∞ (0,T ;V ) + ζ H 2 (0,T ;Y ) , N      k uj − uj −1 2V + ζj − ζj −1 2Y  ck uW 1,∞ (0,T ;V ) + ζ H 2 (0,T ;Y ) . j =1

Finally, in [20] the following estimate was obtained, N −1   1 

ζj − ξ h − ζj +1 − ξ h 2  ch2 ζ 2 1 , j j +1 Y H (0,T ;E) k j =1

which concludes the proof of (37). 4. Numerical examples In this section, we report some numerical results on two-dimensional test problems to show the performance of the numerical method discussed in the previous section. In the examples presented below the long memory operator has the form, t

    G t − s, ε u(s) , ζ (s) ds =

0

t

   ζ (s)Φ Aε u(s) ds,

0

where A is the two-dimensional elasticity tensor under the plane stress hypothesis, that is Er E ταβ , 1  α, β  2, τ ∈ Sd . (τ11 + τ22 )δαβ + 1+r 1 − r2 Here, E and r are Young’s modulus and Poisson’s ratio of the material which occupies Ω, respectively, and δαβ is the Kronecker symbol. Moreover, Φ : Sd → Sd is a truncation operator defined by ⎧ if τij > L, ⎨L   if τij ∈ [−L, L], ∀τ = (τij )di,j =1 ∈ Sd , Φ(τ ) ij = τij ⎩ −L if τ < −L. (Aτ )αβ =

ij

Here, L > 0 is a given constant. We notice that the existence of L is justified taking into account that the small displacement theory is used. In all the examples, value L = 1000 is employed. Note that this choice of function G shows that the more damage the material undergoes (ζ decreases), the less memory effects arise and the more the material behaves like a purely elastic material. Recall also the definition of the von Mises norm for a plane stress field τ = (ταβ ) as,   2 2 2 1/2 + τ22 − τ11 τ22 + 3τ12 . τ V M = τ11 The damage source function is given by [9],       1 − η(ζ ) 1 − λ2 Ψq ∗ (u) − λ3 , φ ε(u), ζ = − λ1 η(ζ ) 2 +

(38)

where Ψq ∗ (u) = min{ε(u) · ε(u), q ∗ } for some constant q ∗ > 0 and η : R → R is a function defined by  ζ if ζ > ζ∗ , η(ζ ) = ζ∗ if ζ  ζ∗ . Here, λ1 , λ2 and λ3 are process parameters whose values in the simulations below are λ1 = 10−2 ,

λ2 = 10,

λ3 = 0.

Also, the value κ = 10−2 is chosen. The following normal compliance function is employed in the simulations, p(r) =

1 r+ , μ

(39)

1286

M. Campo et al. / Applied Numerical Mathematics 58 (2008) 1274–1290

where r+ = max{0, r} and μ is a positive constant which represents a deformability coefficient (value μ = 10−10 is taken). Finally, we use affine elements for the finite element spaces V h and E h (see the definitions (28)–(29)). 4.1. First example In order to verify the asymptotic behavior of the numerical method presented in the previous section, a sequence of numerical solutions was computed following the physical setting described in Fig. 2. Uniform partitions were considered of both the time interval and the spatial domain (each side of the square [0, 1] × [0, 1] was divided into n equal parts), and the corresponding solutions compared with the “exact solution” obtained with k = 0.0005 and n = 256. The boundary ΓD = {0} × [0, 1] was supposed to be fixed, ΓC = [0, 1] × {0} was in frictionless contact with a deformable foundation and ΓF was divided into two parts: {1} × [0, 1], where a density of surface tractions f 2 acted, and [0, 1] × {1}, which was traction-free. No volume forces were supposed to act over the body. The following data were used in this case: T = 1 s,

f 0 = 0 N/m3 ,

f 2 = (−10, −10) N/m2 ,

r = 0.3,

E = 100 N/m2 ,

ζ∗ = 0.01.

The numerical errors given by

  max un − uhk + ζn − ζ hk n

0nN

n

V

Y

are shown in Table 1. Moreover, the linear convergence of the algorithm with respect to h + k, stated in Corollary 9, is clearly observed in Fig. 3. 4.2. Second example As a second example we considered a rectangular domain [0, 2] × [0, 3]. The body was fixed on ΓD = [0, 3] × {0}, the contact occurred on ΓC = {0} × [0, 1] and surface tractions were prescribed on {2} × [0, 3], while {0} × [1, 2] ∪ [0, 3] × {2} was traction-free (see Fig. 4). The following data were used in this example:

Fig. 2. Example 1: Physical setting and mesh for n = 8. Table 1 Example 1: Numerical errors obtained with different n and k n 4 8 16 32 64 128

k 0.05

0.02

0.01

0.005

0.002

0.001

6.779e − 2 3.204e − 2 1.574e − 2 9.222e − 3 6.767e − 3 5.879e − 3

6.665e − 2 3.032e − 2 1.353e − 2 6.560e − 3 3.793e − 3 2.682e − 3

6.628e − 2 2.978e − 2 1.285e − 2 5.757e − 3 2.863e − 3 1.671e − 3

6.610e − 2 2.951e − 2 1.254e − 2 5.390e − 3 2.429e − 3 1.186e − 3

6.610e − 2 2.935e − 2 1.235e − 2 5.187e − 3 2.197e − 3 9.101e − 4

6.596e − 2 2.930e − 2 1.299e − 2 5.123e − 3 2.128e − 3 8.289e − 4

M. Campo et al. / Applied Numerical Mathematics 58 (2008) 1274–1290

Fig. 3. Example 1: Evolution of the error with respect to k + h.

1287

Fig. 4. Example 2: Physical setting.

Fig. 5. Example 2: von Mises stress norm at initial and final times over the deformed configuration.

T = 2 s,

f 0 = 0 N/m3 ,

E = 1000 N/m2 ,

r = 0.3,

f 2 = (−30, 0) N/m2 , ζ∗ = 0.01.

In Fig. 5 the von Mises stress norm is plotted over the deformed configuration at initial time (left-hand side) and at final time (right-hand side). Since the force is assumed to be constant, the effect of the memory can be seen: at final time the body tends to recover the stress-free state. Also, in Fig. 6 the damage field at final time is shown over the deformed configuration. As we expected, the most stressed areas coincide with the most damaged ones. 4.3. Third example In this last example we computed the setting depicted in Fig. 7 using a different memory function than in the previous examples, t 0

    G t − s, ε u(s) , ζ (s) ds = −

t 0

   ζ (s)e−(t−s) Φ Aε u(s) ds.

1288

M. Campo et al. / Applied Numerical Mathematics 58 (2008) 1274–1290

Fig. 6. Example 2: Damage at final time over the deformed configuration.

Fig. 7. Example 3: Physical setting.

Fig. 8. Example 3: von Mises stress norm at initial time and at t = 0.5 s.

Fig. 9. Example 3: von Mises stress norm and damage field at final time.

In order to see the effect of this memory function, the applied surface tractions vanished during the process, that is,  (0, −10) N/m2 if 0 < t  0.5, f 2 (x1 , x2 , t) = if 0.5 < t  1. 0 N/m2 The following data were employed in this example: T = 1 s,

f 0 = 0 N/m3 ,

E = 1000 N/m2 ,

r = 0.3,

ζ∗ = 0.01.

The von Mises stress norm is plotted at initial time and at time t = 0.5 over the corresponding deformed configurations in Fig. 8. We notice that, with this memory function, the body tends to keep the previous deformed configuration. Due to this fact, at time t = 0.5 the forces are the same than at initial time but the deformation has grown. Finally, we can observe in Fig. 9 the von Mises stress norm (left-hand side) and the damage field (right-hand side) at final time. Although there are not forces, at final time the configuration is not the stress-free one because of the influence of the

M. Campo et al. / Applied Numerical Mathematics 58 (2008) 1274–1290

1289

memory, trying to keep the previous deformations. Again, the most damaged areas coincide with the most stressed ones. Acknowledgement This work was partially supported by MCYT-Spain (Project BFM2003-05357). It is also part of the project “New materials, Adaptive Systems and their Nonlinearities; Modelling, Control and Numerical Simulation” carried out in the framework of the European community program “Improving the Human Research Potential and the Socio-Economic Knowledge Base” (Contract No. HPRN-CT-2002-00284). References [1] K.T. Andrews, J.R. Fernández, M. Shillor, Numerical analysis of dynamic thermoviscoelastic contact with damage of a rod, IMA J. Appl. Math. 70 (6) (2005) 768–795. [2] T.A. Angelov, On a rolling problem with damage and wear, Mech. Res. Comm. 26 (1999) 281–286. [3] A. Bermúdez, C. Moreno, Duality methods for solving variational inequalities, Comp. Math. Appl. 7 (1981) 43–58. [4] E. Bonetti, M. Frémond, Damage theory: microscopic effects of vanishing macroscopic motions, Comp. Appl. Math. 22 (3) (2003) 313–333. [5] E. Bonetti, G. Schimperna, Local existence for Frémond’s model of damage in elastic materials, Comp. Mech. Thermodyn. 16 (4) (2004) 319–335. [6] M. Campo, J.R. Fernández, W. Han, M. Sofonea, A dynamic viscoelastic contact problem with normal compliance and damage, Finite Elem. Anal. Des. 42 (2005) 1–24. [7] M. Campo, J.R. Fernández, T.-V. Hoarau-Mantel, Analysis of two frictional viscoplastic contact problems with damage, J. Comput. Appl. Math. 196 (1) (2006) 180–197. [8] M. Campo, J.R. Fernández, K.L. Kuttler, An elastic-viscoplastic contact problem with damage, Comput. Methods Appl. Mech. Engrg. 196 (33–34) (2007) 3219–3229. [9] M. Campo, J.R. Fernández, K.L. Kuttler, M. Shillor, Quasistatic evolution of damage in an elastic body: numerical analysis and computational experiments, Appl. Numer. Math. 57 (9) (2007) 975–988. [10] O. Chau, J.R. Fernández, A convergence result in elastic-viscoplastic contact problems with damage, Math. Comp. Modelling 37 (2003) 301–321. [11] O. Chau, J.R. Fernández, W. Han, M. Sofonea, A frictionless contact problem for elastic-viscoplastic materials with normal compliance and damage, Comput. Methods Appl. Mech. Engrg. 191 (2002) 5007–5026. [12] J. Chen, W. Han, M. Sofonea, Numerical analysis of a nonlinear evolutionary system with applications in viscoplasticity, SIAM J. Numer. Anal. 38 (4) (2000) 1171–1199. [13] P.G. Ciarlet, The finite element method for elliptic problems, in: P.G. Ciarlet, J.L. Lions (Eds.), Handbook of Numerical Analysis, vol. II, Part 1, North-Holland, 1991, pp. 17–352. [14] G. Duvaut, J.L. Lions, Inequalities in Mechanics and Physics, Springer-Verlag, Berlin, 1976. [15] I. Figueiredo, L. Trabucho, A class of contact and friction dynamic problems in thermoelasticity and in thermoviscoelasticity, Internat. J. Engrg. Sci. 33 (1) (1995) 45–66. [16] M. Frémond, Non-smooth Thermomechanics, Springer, Berlin, 2002. [17] M. Frémond, B. Nedjar, Damage in concrete: the unilateral phenomenon, Nuclear Eng. Design 156 (1995) 323–335. [18] M. Frémond, B. Nedjar, Damage, gradient of damage and principle of virtual work, Internat. J. Solids Structures 33 (1996) 1083–1103. [19] R. Glowinski, Numerical Methods for Nonlinear Variational Problems, Springer-Verlag, New York, 1984. [20] W. Han, M. Shillor, M. Sofonea, Variational and numerical analysis of a quasistatic viscoelastic problem with normal compliance, friction and damage, J. Comput. Appl. Math. 137 (2001) 377–398. [21] W. Han, M. Sofonea, Quasistatic Contact Problems in Viscoelasticity and Viscoplasticity, American Mathematical Society–Intl. Press, 2002. [22] A. Klarbring, A. Mikeli´c, M. Shillor, Frictional contact problems with normal compliance, Internat. J. Engrg. Sci. 26 (1988) 811–832. [23] K.L. Kuttler, Quasistatic evolution of damage in an elastic-viscoplastic material, Electron. J. Differential Equations 147 (2005) 1–25. [24] K.L. Kuttler, M. Shillor, Quasistatic evolution of damage in an elastic body, Nonlinear Anal. Real World. Appl. 7 (4) (2006) 674–699. [25] J. Lemaitre, J.L. Chaboche, Mécanique des matériaux solides, Dunod, Paris, 1988. [26] R. Liebe, P. Steinmann, A. Benallal, Theoretical and numerical aspects of a thermodynamically consistent framework for geometrically linear gradient damage, Comput. Methods Appl. Mech. Engrg. 190 (2001) 6555–6576. [27] J.L. Lions, E. Magenes, Problèmes aux limites non homogènes et applications, vol. 1, Dunod, Paris, 1968. [28] C. Miehe, Discontinuous and continuous damage evolution in Ogden-type large-strain elastic materials, Eur. J. Mech. A Solids 14 (3) (1995) 697–720. [29] B. Nedjar, Elastoplastic-damage modelling including the gradient of damage. Formulation and computational aspects, Internat. J. Solids Structures 38 (2001) 5421–5451. [30] A.C. Pipkin, Lectures in Viscoelasticity Theory, Applied Mathematical Sciences, vol. 7, George Allen & Unwin Ltd., Springer-Verlag, London, New York, 1972. [31] A. Rodríguez-Arós, J.M. Viaño, M. Sofonea, ‘A class of evolutionary variational inequalities with Volterra-type term, Math. Models Methods Appl. Sci. 14 (4) (2004) 557–577.

1290

M. Campo et al. / Applied Numerical Mathematics 58 (2008) 1274–1290

[32] M. Shillor, M. Sofonea, J.J. Telega, Models and Variational Analysis of Quasistatic Contact, Lecture Notes in Physics, vol. 655, Springer, Berlin, Heidelberg, 2004. [33] M. Sofonea, A. Rodríguez-Arós, J.M. Viaño, A class of integro-differential variational inequalities with applications to viscoelastic contact, Math. Comput. Modelling 41 (11–12) (2005) 1355–1369. [34] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method, vol. 2, Solid Mechanics, McGraw-Hill, London, 1989.