On a dynamic contact problem for elastic-visco-plastic materials

On a dynamic contact problem for elastic-visco-plastic materials

Applied Numerical Mathematics 57 (2007) 498–509 www.elsevier.com/locate/apnum On a dynamic contact problem for elastic-visco-plastic materials Weimin...

194KB Sizes 0 Downloads 28 Views

Applied Numerical Mathematics 57 (2007) 498–509 www.elsevier.com/locate/apnum

On a dynamic contact problem for elastic-visco-plastic materials Weimin Han a,∗ , M. Sofonea b a Department of Mathematics, University of Iowa, Iowa City, IA 52242, USA b Laboratoire de Mathématiques et Physique pour les Systèmes, Université de Perpignan,

52 Avenue Paul Alduy, 66 860 Perpignan, France Available online 2 August 2006

Abstract In this paper, we study a mathematical problem for dynamic contact between an elastic-visco-plastic body and an obstacle. The contact is frictionless, modelled with a normal compliance condition involving adhesion effect of contact surfaces. Evolution of the bonding field is described by a first order differential equation. We provide a weak formulation of the contact problem, in the form of an integro-differential system, and present an existence and uniqueness result for its solution. We then introduce and study a fully discrete scheme for solving the problem. While it is possible to show the convergence of the scheme without assuming additional solution regularities, we focus on the derivation of optimal order error estimates under certain solution regularity assumptions. © 2006 IMACS. Published by Elsevier B.V. All rights reserved. MSC: 65M60; 74M15; 74H20 Keywords: Elastic-visco-plastic material; Dynamic process; Frictionless contact; Normal compliance; Adhesion; Weak formulation; Fully discrete scheme; Finite element method; Error estimate

1. Introduction The aim of this paper is to study a frictionless contact problem for elastic-visco-plastic materials with a constitutive relation of the form     ˙ σ (t) = Aε u(t) + Eε u(t) +

t

     ˙ G σ (s) − Aε u(s) , ε u(s) ds,

(1.1)

0

where u is the displacement field, σ and ε(u) are the stress and the linearized strain tensor, respectively. Here A and E are nonlinear operators describing the purely viscous and the elastic properties of the material, respectively, and G is a nonlinear constitutive function describing the visco-plastic behaviour of the material. We use dots for derivatives with respect the time variable t. It follows from (1.1) that at each time moment t, the stress tensor σ (t) is split into two

* Corresponding author.

E-mail address: [email protected] (W. Han). 0168-9274/$30.00 © 2006 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2006.07.003

W. Han, M. Sofonea / Applied Numerical Mathematics 57 (2007) 498–509

499

˙ parts: σ (t) = σ V (t) + σ R (t), where σ V (t) = Aε(u(t)) represents the purely viscous part of the stress whereas σ R (t) satisfies a rate-type elastic-visco-plastic relation,   σ (t) = Eε u(t) +

t

R

   G σ R (s), ε u(s) ds.

(1.2)

0

When G = 0, (1.1) reduces to the Kelvin–Voigt viscoelastic constitutive relation, ˙ + Eε(u). σ = Aε(u)

(1.3)

Examples and mechanical interpretation of elastic-visco-plastic materials of the form (1.2) can be found in [5,13]. Quasistatic contact problems for materials of the form (1.2) or (1.3) are the topic of numerous papers, e.g. [1,2,7,20], and the comprehensive reference [12]. Dynamic contact problems with Kelvin–Voigt materials of the form (1.3) are studied in [14–16]. More detail can be found in the recent monograph [6]. The normal compliance contact condition was first considered in [18] in the study of dynamic problems with linearly elastic and viscoelastic materials and then it was used in various references, see e.g. [17,20]. This condition allows the interpenetration of the body’s surface into the obstacle and it was justified by considering the interpenetration and deformation of surface asperities. Processes of adhesion are important in many industrial settings where parts, usually nonmetallic, are glued together. For this reason, adhesive contact between bodies, when a glue is added to prevent the surfaces from relative motion, has recently received increased attention in the literature. General models with adhesion can be found in [8,9]. Results on the mathematical analysis of various adhesive contact problems can be found in [21,22]. The idea is the introduction of a surface internal variable, the bonding field β ∈ [0, 1], which describes the fractional density of active bonds on the contact surface. At a point on the bonding contact surface Γ3 , when β = 1, the adhesion is complete and all the bonds are active; when β = 0 all the bonds are inactive, severed, and there is no adhesion; when 0 < β < 1 the adhesion is partial and only a fraction β of the bonds is active. In this paper we study a dynamic frictionless contact problem for rate-type materials of the form (1.1). The contact is frictionless and is modelled with normal compliance and adhesion. The paper is organized as follows. In Section 2 we describe the contact problem, list assumptions on the data and provide a weak formulation; then we state and prove an existence and uniqueness result for solution of the weak formulation; the proof is based on arguments of nonlinear evolutionary equations and fixed point. In Section 3 we study a fully-discrete scheme for numerical solution; we derive optimal order error estimates under certain solution regularity assumptions. 2. The contact problem We first describe the physical setting of the contact problem. An elastic-visco-plastic body of density ρ occupies a bounded domain Ω ⊂ Rd (d = 2, 3) with a Lipschitz boundary Γ that is partitioned into three disjoint measurable parts Γ1 , Γ2 and Γ3 , such that meas(Γ1 ) > 0. We are interested in the deformation of the body in a time interval [0, T ]. At any time, the body is clamped on Γ1 , a volume force of density f0 acts in Ω, and a surface traction of density f2 acts on Γ2 . In the reference configuration the body is in adhesive contact on Γ3 with an obstacle/foundation; the contact is frictionless and it is modelled with normal compliance. The classical formulation of the dynamic contact problem is the following. Problem 2.1. Find a displacement field u : Ω × [0, T ] → Rd , a stress field σ : Ω × [0, T ] → Sd , and a bonding field β : Γ3 × [0, T ] → R such that     ˙ σ (t) = Aε u(t) + Eε u(t) +

t

     ˙ G σ (s) − Aε u(s) , ε u(s) ds

in Ω × (0, T ),

(2.1)

ρ u¨ = Div σ + f0

in Ω × (0, T ),

(2.2)

u=0

on Γ1 × (0, T ),

(2.3)

σ ν = f2

on Γ2 × (0, T ),

(2.4)

0

500

W. Han, M. Sofonea / Applied Numerical Mathematics 57 (2007) 498–509

−σν = pν (uν ) − γν β 2 Rν (uν )

on Γ3 × (0, T ),

(2.5)

−σ τ = pτ (β)Rτ (uτ ) 2      β˙ = − β γν Rν (uν )2 + γτ Rτ (uτ ) − a

onΓ3 × (0, T ),

(2.6)

on Γ3 × (0, T ),

(2.7)

in Ω,

(2.8)

on Γ3 .

(2.9)

u(0) = u0 ,

+

˙ u(0) = v0

β(0) = β0

Here and below Sd denotes the space of second order symmetric tensors on Rd , whereas “·” and  ·  represent the inner product and the Euclidean norm on Sd and Rd , respectively; ν is the unit outer normal vector on Γ , and the subscripts ν and τ refer to normal and tangential components of vectors and tensors; r+ = max{r, 0} denotes the positive part of r and, finally, u0 , v0 , β0 are given initial displacement, velocity and bonding fields. We briefly comment on the contact conditions (2.5)–(2.7). Condition (2.5) represents the normal compliance condition with adhesion. Condition (2.6) is the tangential boundary condition on the contact surface Γ3 , showing that the shear on the contact surface depends on the adhesion field and on the tangential displacement. pν and pτ are given functions,   L if s < −L, v if v  L, Rτ (v) = L v if v > L Rν (s) = −s if − L  s  0, v 0 ifs > 0, with L > 0 being a characteristic length of the bond, beyond which there is no any additional traction (see, e.g., [19]). Eq. (2.7) describes the evolution of the bonding field with given material parameters γν , γτ and a . It is easy to see that if 0  β0  1 a.e. on Γ3 , then 0  β  1 a.e. on Γ3 during the process. Before introducing a weak formulation for Problem 2.1, we present some notation and preliminary material. For further details, we refer the reader to [12,22]. Everywhere below the indices i and j range between 1 to d, summation over repeated indices is implied and the index following a comma indicates the partial derivative with respect to the corresponding component of the independent spatial variable. For the stress and strain variables, we use the real spaces   Q1 = σ ∈ Q: Div σ = (σij,j ) ∈ L2 (Ω)d ; Q = σ = (σij ): σij = σj i ∈ L2 (Ω) , these are Hilbert spaces with the inner products  (σ , τ )Q1 = (σ , τ )Q + (Div σ , Div τ )L2 (Ω)d . (σ , τ )Q = σij τij dx, Ω

The associated norms are denoted by  · Q and  · Q1 . For the displacement variable, we use the real Hilbert space H1 = {u = (ui ) ∈ L2 (Ω)d : ε(u) ∈ Q} endowed with the inner product (u, v)H1 = (u, v)L2 (Ω)d + (ε(u), ε(v))Q and associated norm  · H1 . Here ε is the deformation operator: ε(u) = (εij (u)), εij (u) = (ui,j + uj,i )/2. For the damage field, we use the space B = L2 (Γ3 ) with the inner product (·, ·)B = (·, ·)L2 (Γ3 ) and the norm  · B =  · L2 (Γ3 ) . For v ∈ H1 we also use v to denote the trace of v on Γ and denote by vν = v · ν and vτ = v − vν ν the normal and the tangential components of v on Γ . Similarly, σν and σ τ denote the normal and the tangential traces of a function σ ∈ Q1 . When σ is a regular function, then σν = (σ ν) · ν, σ τ = σ ν − σν ν, and the following Green’s formula holds:    (2.10) σ , ε(v) Q + (Div σ , v)L2 (Ω)d = σ ν · v da ∀v ∈ H1 . Γ

Define V = {v ∈ H1 : v = 0 on Γ1 }. Since meas(Γ1 ) > 0, Korn’s inequality holds:   ε(v)  cK vH ∀v ∈ V , Q

1

(2.11)

where cK > 0 is a constant depending only on Ω and Γ1 . A proof of Korn’s inequality can be found in, e.g. [6, pp. 16– 17]. The space V is a real Hilbert space endowed with the inner product (u, v)V = (ε(u), ε(v))Q and the associated norm  · V . By (2.11), the norms  · H1 and  · V are equivalent on V . Given a real Banach space (X,  · X ) and T > 0, we will use real Banach spaces C([0, T ]; X), C 1 ([0, T ]; X), p L (0, T ; X), W k,p (0, T ; X) (k  0 integer, p ∈ [1, ∞]), and H 1 (0, T ; X) ≡ W 1,2 (0, T ; X). Also, we will use the set  Q = θ ∈ L∞ (0, T ; B): 0  θ (t)  1 ∀t ∈ [0, T ], a.e. on Γ3 .

W. Han, M. Sofonea / Applied Numerical Mathematics 57 (2007) 498–509

501

We now list assumptions on the data. Assume the operators A, E and G satisfy the following conditions (LA , mA , LE and LG being positive constants). ⎫ (a) A : Ω × Sd → Sd . ⎪ ⎪ ⎪ ⎪ (b) A(x, ε1 ) − A(x, ε2 )  LA ε1 − ε 2  ∀ε1 , ε2 ∈ Sd , a.e. x ∈ Ω. ⎬ (2.12) (c) (A(x, ε 1 ) − A(x, ε 2 )) · (ε 1 − ε2 )  mA ε1 − ε2 2 ∀ε 1 , ε2 ∈ Sd , a.e. x ∈ Ω. ⎪ d ⎪ ⎪ (d) For any ε ∈ S , x → A(x, ε)is measurable on Ω. ⎪ ⎭ (e) The mapping x → A(x, 0) belongs to Q. ⎫ (a) E : Ω × Sd → Sd . ⎪ ⎪ ⎬ (b) E(x, ε1 ) − E(x, ε 2 )  LE ε1 − ε2  ∀ε 1 , ε2 ∈ Sd , a.e. x ∈ Ω. (2.13) (c) For any ε ∈ Sd , x → E(x, ε) is measurable on Ω. ⎪ ⎪ ⎭ (d) The mapping x → E(x, 0) belongs to Q. ⎫ (a) G : Ω × Sd × Sd → Sd . ⎪ ⎪ ⎪ (b) G(x, σ 1 , ε 1 ) − G(x, σ 2 , ε 2 )  LG (σ 1 − σ 2  + ε1 − ε2 ) ⎪ ⎬ d (2.14) ∀σ 1 , σ 2 , ε1 , ε 2 ∈ S , a.e. x ∈ Ω. ⎪ ⎪ (c) For any σ , ε ∈ Sd , x → G(x, σ , ε) is measurable on Ω. ⎪ ⎪ ⎭ (d) The mapping x → G(x, 0, 0) belongs to Q. The normal compliance function pν and the tangential function pτ satisfy the assumptions (Lν , Lτ and Mτ being positive constants) ⎫ (a) pν : Γ3 × R → R+ . ⎪ ⎪ ⎪ ⎪ (b) |pν (x, r1 ) − pν (x, r2 )|  Lν |r1 − r2 | ∀r1 , r2 ∈ R, a.e. x ∈ Γ3 . ⎬ (c) (pν (x, r1 ) − pν (x, r2 ))(r1 − r2 )  0 ∀r1 , r2 ∈ R, a.e. x ∈ Γ3 . (2.15) ⎪ ⎪ (d) For any r ∈ R, x → pν (x, r) is measurable on Γ3 . ⎪ ⎪ ⎭ (e) pν (x, r) = 0 for all r  0, a.e. x ∈ Γ3 . ⎫ (a) pτ : Γ3 × R −→ R+ . ⎪ ⎪ ⎪ (b) |pτ (x, β1 ) − pτ (x, β2 )|  Lτ |β1 − β2 | ∀β1 , β2 ∈ R, a.e. x ∈ Γ3 . ⎪ ⎬ (c) |pτ (x, β)|  Mτ ∀β ∈ R, a.e. x ∈ Γ3 . (2.16) ⎪ ⎪ (d) For any β ∈ R, x → pτ (x, β) is measurable on Γ3 . ⎪ ⎪ ⎭ (e) The mapping x → pτ (x, 0) belongs to L2 (Γ3 ). We suppose that the mass density satisfies ρ ∈ L∞ (Ω), there exists ρ ∗ > 0 such that ρ(x)  ρ ∗ a.e. x ∈ Ω,

(2.17)

and the body forces and surface tractions have the regularity     f0 ∈ L2 0, T ; L2 (Ω)d , f2 ∈ L2 0, T ; L2 (Γ2 )d .

(2.18)

The adhesion coefficients γν , γτ and a satisfy the conditions γν , γτ ∈ L∞ (Γ3 ),

a ∈ L2 (Γ3 ),

γ ν , γτ , a  0

a.e. on Γ3

(2.19)

and, finally, the initial data satisfy u0 ∈ V ,

v0 ∈ L2 (Ω)d ,

β0 ∈ B,

0  β0  1 a.e. on Γ3 .

(2.20) 1/2

We will use a modified inner product on H = L2 (Ω)d : ((u, v))H = (ρu, v)H . Let |v|H = (ρ v, v)H for v ∈ H . By assumption (2.17), | · |H and  · H are equivalent norms on H . The embedding of (V ,  · V ) into (H, | · |H ) is continuous and dense. We denote by V the dual space of V and identify H with its own dual. Then u, v V ×V = ((u, v))H ∀u ∈ H, v ∈ V . Assumptions (2.18) allow us, for a.e. t ∈ (0, T ), to define f(t) ∈ V by     f(t), v V ×V = f0 (t) · v dx + f2 (t) · v da ∀v ∈ V . (2.21) Ω

Γ2

502

W. Han, M. Sofonea / Applied Numerical Mathematics 57 (2007) 498–509

Also, define j : L∞ (Γ3 ) × V × V → R by the formula     j (β, u, v) = pν (uν ) − γν β 2 Rν (uν ) vν + pτ (β)Rτ (uτ ) · vτ da.

(2.22)

Γ3

The integral is well-defined due to the assumptions (2.15), (2.16) and (2.19). To derive a weak formulation for Problem 2.1, assume u, σ and β are smooth functions satisfying (2.1)–(2.9) and let t ∈ [0, T ]. We take the dot product of Eq. (2.2) with an arbitrary w ∈ V , integrate over Ω, and use Green’s formula (2.10) to obtain       ¨ ρ u(t), w H + σ (t), ε(w) Q = f0 (t) · w dx + σ (t)ν · w da. (2.23) Ω

Γ

Applying the boundary conditions (2.4)–(2.6) and noting that w = 0 on Γ1 , we have       pν (uν ) − γν β 2 Rν (uν ) wν da − pτ (β) Rτ (uτ ) · wτ da. σ (t)ν · w da = f2 (t) · w dx − Γ

Γ2

Γ3

Γ3

Use this relation in (2.23) and recall (2.21) and (2.22) to find         ¨ u(t), w V ×V + σ (t), ε(w) Q + j β(t), u(t), w = f(t), w V ×V . Thus, we have the following weak formulation of Problem 2.1. Problem 2.2. Find a displacement field u : [0, T ] → V , a stress field σ : [0, T ] → Q, and a bonding field β : [0, T ] → L∞ (Γ3 ) such that, for a.e. t ∈ (0, T ),     ˙ σ (t) = Aε u(t) + Eε u(t) +

t

     ˙ G σ (s) − Aε u(s) , ε u(s) ds,

(2.24)

0



       ¨ u(t), w V ×V + σ (t), ε(w) Q + j β(t), u(t), w = f(t), w V ×V       ˙ = − β γν Rν (uν )2 + γτ R(uτ )2 − a , β(t) +

∀w ∈ V ,

(2.25) (2.26)

and u(0) = u0 ,

˙ u(0) = v0 ,

β(0) = β0 .

(2.27)

We have the following existence and uniqueness result for Problem 2.2. Theorem 2.3. Assume (2.12)–(2.20). Then Problem 2.2 has a unique solution and   u¨ ∈ L2 (0, T ; V ), u ∈ H 1 (0, T ; V ) ∩ C 1 [0, T ]; H , σ ∈ L (0, T ; Q), 2

β ∈W

1,∞



Div σ ∈ L (0, T ; V ), 2

(2.29)

(0, T ; B) ∩ Q.

(2.30)

Proof. The arguments are similar to those presented in [22], and so we only give a sketch. (i) Given η ∈ L2 (0, T ; V ), consider the problem of finding uη : [0, T ] → V such that         u¨ η (t), w V ×V + Aε u˙ η (t) , ε(w) Q + η(t), w V ×V = f(t), w V ×V ∀w ∈ V , a.e. t ∈ (0, T ), uη (0) = u0 ,

(2.28)

u˙ η (0) = v0 .

(2.31) (2.32)

By a result of [3, p. 140], there is a unique solution uη with the regularity (2.28). Next, we prove that there exists a unique solution σ η ∈ H 1 (0, T ; Q) of the equation   σ η (t) = Eε uη (t) +

t 0

   G σ η (s), ε uη (s) ds

∀t ∈ [0, T ].

(2.33)

W. Han, M. Sofonea / Applied Numerical Mathematics 57 (2007) 498–509

503

Consider the operator Λη : L2 (0, T ; Q) → L2 (0, T ; Q) given by 

 Λη σ (t) = Eε uη (t) +

t

   G σ (s), ε uη (s) ds

∀σ ∈ L2 (0, T ; Q), t ∈ [0, T ].

0

For σ 1 , σ 2

∈ L2 (0, T ; Q)

we use (2.14) to obtain

  Λη σ 1 (t) − Λη σ 2 (t)  LG Q

t

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

∀t ∈ [0, T ].

0

It follows from this inequality that a high enough power of the operator Λη is a contraction on the Banach space L2 (0, T ; Q) and therefore an application of a variant of Banach’s fixed point theorem [12, Corollary 1.63] shows that (2.33) has a unique solution σ η ∈ L2 (0, T ; Q). Moreover, the regularity of uη and the properties of the operators E and G imply that σ η ∈ H 1 (0, T ; Q). Then, consider the problem of finding βη : [0, T ] → B such that       2 2  β˙η (t) = − βη (t) γν Rν uην (t) + γτ Rτ uητ (t)  − a + a.e. t ∈ (0, T ), (2.34) βη (0) = β0 .

(2.35)

It has a unique solution βη ∈ W 1,∞ (0, T ; B) ∩ Q by the Cauchy–Lipschitz theorem. (ii) Now define an operator Λ : L2 (0, T ; V ) → L2 (0, T ; V ) by the formula  t           G σ η (s), ε uη (s) ds, ε(w) Λη(t), w V ×V = Eε uη (t) , ε(w) Q + Q

0

  + j βη (t), uη (t), w

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

(2.36)

L2 (0, T ; V ),

i = 1, 2, denote by ui , σ i , βi , i = 1, 2, the corresponding functions introduced in (i). For given ηi ∈ Then, for some constant c > 0 and for any t ∈ [0, T ], t

  u˙ 1 (s) − u˙ 2 (s)2 ds  c V

0

t

  η1 (t) − η2 (t)2 ds, V

0

  t       σ 1 (t) − σ 2 (t)  c u1 (t) − u2 (t) + u1 (s) − u2 (s) ds , Q V V 0

  β1 (t) − β2 (t)  c B

t

  u1 (s) − u2 (s) ds. V

0

We use these bounds and the assumptions on the data to show that a power of the operator Λ is a contraction on the space L2 (0, T ; V ). So by the variant of Banach’s fixed point theorem, Λ has a unique fixed point η∗ ∈ L2 (0, T ; V ). ˙ + σ η∗ , and β = βη∗ . Clearly, equalities (2.24) and (2.26) hold from (2.33) and (2.34). Let u = uη∗ , σ = Aε(u) Moreover, since η∗ = Λ η∗ , it follows from (2.31), (2.33) and (2.36) that (2.25) holds, too. The regularity of the solution expressed in (2.28) follows from step (i). Since u ∈ H 1 (0, T ; V ) and σ η∗ ∈ H 1 (0, T ; Q), the definition of σ and (2.12) yield σ ∈ L2 (0, T ; Q). Choosing now w = ϕ in (2.25), where ϕ ∈ C0∞ (Ω)d , and using (2.21), (2.22), we obtain ¨ = Div σ (t) + f0 (t) ρ u(t)

a.e. t ∈ (0, T ). u¨ ∈ L2 (0, T ; V )

(2.37) ∈ L2 (0, T ; V ).

and (2.37) imply that Div σ Also, the Now, assumptions (2.17), (2.18), the fact that regularity β ∈ W 1,∞ (0, T ; B) ∩ Q follows from step (i). So (u, σ , β) is a solution of Problem 2.2 and it satisfies (2.28)–(2.30). This proves the existence of a solution. Uniqueness of the solution follows from the uniqueness of the fixed point of Λ and from the unique solvability of problems in step (i). 2

504

W. Han, M. Sofonea / Applied Numerical Mathematics 57 (2007) 498–509

3. Numerical approximation In this section we study a fully discrete scheme for the numerical approximation of Problem 2.2. Introduce t ˙ the velocity variable v(t) = u(t). Then u(t) = u0 + 0 v(s) ds, t ∈ [0, T ]. It follows from Theorem 2.3 that v ∈ L2 (0, T ; V ) ∩ C([0, T ]; H ), v˙ ∈ L2 (0, T ; V ), and         v˙ (t), w V ×V + σ (t), ε(w) Q + j β(t), u(t), w = f(t), w V ×V ∀w ∈ V , a.e. t ∈ (0, T ), (3.1) v(0) = v0 .

(3.2)

For simplicity, assume that Ω is a polygonal/polyhedral domain and its boundary Γ is split into three relatively ij (i) Γj such that on closed subsets Γ1 , Γ2 , and Γ3 , with mutually disjoint interiors. For 1  j  3, write Γj = i=1 (i)

(i)

each Γj , the unit outward normal vector is constant (i.e., Γj is a straight or planar part of the boundary). Let {T h } be a regular family of finite element partitions of Ω, assumed to be compatible with the decomposition Γ = 3 ij (i) (i) j =1 i=1 Γj , i.e., if S is an element side such that for some j and i, S ∩ Γj contains an interior point of S, then (i)

S ⊂ Γj . Corresponding to a partition T h , we use  V h = vh ∈ C(Ω)d : vh |K linear ∀K ∈ T h , vh = 0 on Γ1

(3.3)

to approximate the space V . For w ∈ H 2 (Ω)d ∩ V , let Π h w ∈ V h be the finite element interpolant of w in V h . Then, we have the following interpolation error estimate:     w − Π h w 2 d + hw − Π h w 1 d  ch2 |w| 2 d . (3.4) H (Ω) L (Ω) H (Ω) We use

 Qh = τ h ∈ Q: τ h |K ∈ Sd ∀K ∈ T h

(3.5)

to approximate the space Q. Let PQh : Q → Qh be the orthogonal projection operator:     PQh q, qh Q = q, qh Q ∀q ∈ Q, qh ∈ Qh . Then we have the following error estimate: τ − PQh τ Q  ch|τ |H 1 (Ω)d×d

∀τ ∈ H 1 (Ω)d×d .

(3.6)

Denote by TΓh3 the partition of Γ3 induced by the triangulation T h . Then the space B is approximated by  B h = θ h ∈ B: θ h |γ ∈ R ∀γ ∈ TΓh3 .

(3.7)

We define PB h : B → B h to be the orthogonal projection operator in B:     PB h θ, θ h B = θ, θ h B ∀θ ∈ B, θ h ∈ B h . Then we have the following error estimate: θ − PB h θ L2 (Γ3 )  ch|θ |H 1 (Γ3 )

∀θ ∈ H 1 (Γ3 ).

(3.8)

Note that the orthogonal projection operators are nonexpansive: PQh qQ  qQ

∀q ∈ Q;

PB h θ B  θ B

∀θ ∈ B.

Now we introduce a fully discrete numerical scheme for Problem 2.2. The discussion will be given only for the case of uniform partitions of the time interval, although it is readily extendable to the case with an arbitrary time interval partition. Let N be a positive integer, and let 0 = t0 < t1 < · · · < tN = T be a uniform partition of the time interval [0, T ]. Denote k = T /N the stepsize, tn = nk, 0  n  N . Let uh0 , vh0 ∈ V h be approximations of u0 and v0 , and let β0h ∈ B h be an approximation of β. We make the following additional assumptions on the data: f0 ∈ C([0, T ]; L2 (Ω)d ), f2 ∈ C([0, T ]; L2 (Γ2 )d ). Then f, defined in (2.21), satisfies f ∈ C([0, T ]; V ). For a continuous function z(t), we let zn = z(tn ). For a sequence {zn }N n=0 , we let δzn = (zn − zn−1 )/k, 1  n  N .

W. Han, M. Sofonea / Applied Numerical Mathematics 57 (2007) 498–509

505

N h hk = {σ hk }N ⊂ Qh , and a Problem 3.1. Find a discrete velocity field vhk = {vhk n }n=0 ⊂ V , a discrete stress field σ n n=1 h discrete bonding field β hk = {βnhk }N n=0 ⊂ B with h vhk 0 = v0 ,

β0hk = β0h

(3.9)

such that for n = 1, 2, . . . , N , n−1  hk   hk     hk   hk  + P + , = P Aε v Eε u kPQh G σ hk σ hk h h Q Q n n j − Aε vj , ε uj n−1

(3.10)

j =0

 hk h    h   hk    h h δvn , w H + σ hk + j βn−1 , uhk n ,ε w n−1 , w = fn , w V ×V Q   hk 2   hk   2   − a .  γν Rν uhk δβnhk = −PB h βn−1 n−1,ν + γτ Rτ un−1,τ + n h hk hk Here, un = u0 + j =1 kvj , 0  n  N .

∀wh ∈ V h ,

(3.11) (3.12)

We can combine (3.10) and (3.11) to obtain   n−1       hk   h     hk h  hk , ε uhk kG σ hk δvn , w H + Aε vn , ε w Q = − Eε(uhk j − Aε vj j n−1 ) + 

+ fn , w

h

 V ×V

j =0

  hk h − j βn−1 , uhk n−1 , w

Q

∀w ∈ V . h

h

hk hk hk h For given fn , uhk n−1 , and βn−1 , this equation has a unique solution vn and (3.12) has a unique solution βn ∈ B . h N hk So Problem 3.1 admits a unique solution. Since β0 L∞ (Γ3 )  1, it follows from (3.12) that {βn L∞ (Γ3 ) }n=0 is uniformly bounded. hk We now turn to an error analysis of the numerical solution by estimating the numerical errors un − uhk n , vn − vn , hk hk σ n − σ n , and βn − βn , for 1  n  N . For this purpose, we make additional assumptions on the regularity of the solution:   (3.13) v ∈ H 1 (0, T ; V ) ∩ H 2 (0, T ; H ) ∩ C [0, T ]; H 2 (Ω)d ,   1,1 1 d×d σ ∈ W (0, T ; Q) ∩ C [0, T ]; H (Ω) , (3.14)   β ∈ W 2,1 (0, T ; B) ∩ C 1 [0, T ]; H˜ 1 (Γ3 ) . (3.15)

i3 (i) (i) Γ3 , where on each Γ3 the unit outward normal vecThe space H˜ 1 (Γ3 ) is defined as follows. Recall that Γ3 = i=1 (i) tor is constant. Then β ∈ H˜ 1 (Γ3 ) if and only if β|Γ (i) ∈ H 1 (Γ3 ), 1  i  i3 . Thus, β belongs to the space if its restric3  i3 (i) (i) tion to Γ3 lies in H 1 (Γ3 ) for each i. It is natural to choose the product norm βH˜ 1 (Γ3 ) = [ i=1 β2 1 (i) ]1/2 . Note that (3.1) can be replaced by         (˙v(t), w) H + σ (t), ε(w) Q + j β(t), u(t), w = f(t), w V ×V

H (Γ3 )

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

(3.16)

From the assumed regularity of the solution we have u0 , v0 ∈ H 2 (Ω)d , β0 ∈ H˜ 1 (Γ3 ). For definiteness, choose be the finite element interpolants of u0 and v0 in V h , and β0h ∈ B h is the orthogonal projection of β0

uh0 , vh0 ∈ V h to on B h . Then

      u0 − uh  + v0 − vh  + β0 − β h   ch. 0 V 0 V 0 B

(3.17)

We recall some relations found in [22] that are useful for error estimation. n    2    2 2  un − uhk 2  c k vj − vhk n V j V +c h +k ,

(3.18)

n−1     2   2 2 un − uhk 2  c  k vj − vhk j V +c h +k , n−1 V

(3.19)

j =1

j =1

506

W. Han, M. Sofonea / Applied Numerical Mathematics 57 (2007) 498–509 n−1    2  2   2   2 βn − β hk 2  c  k βj − βjhk B + vj − vhk n B j V +c h +k ,

(3.20)

n−2    2  2   2   2 βn − β hk 2  c  k βj − βjhk B + vj − vhk j V +c h +k . n−1 B

(3.21)

j =1

j =1

The first two relations are found on pages 60 and 59, and the other two on page 66 of [22]. We will also use the following result [11, Lemma 4.1].  Lemma 3.2. Assume gn  0, en  0 satisfy en  cgn + c n−1 j =1 kej , 1  n  N . Then   n−1  en  c gn + kgj , 1  n  N, max en  c max gn . 1nN

j =1

1nN

In particular, we apply Lemma 3.2 on (3.20) to obtain n−1     2  2  2 βn − β hk 2  c  k vj − vhk n B j V +c h +k .

Denote θ G ,n =

 tn 0

j =1

G(σ (s) − Aε(v(s)), ε(u(s))) ds −

(3.22) n−1

hk j =0 kG(σ j

hk − Aε(vhk j ), ε(uj )). Then

t n−1 j +1         σ (s) − σ j  + v(s) − vj  + u(s) − uj  ds θ G ,n Q  c Q V V j =0 tj

+c

n−1         hk  hk     k σ j − σ hk j Q + vj − vj V + uj − uj V . j =0

˙ L1 (0,T ;V ) . The first sum can be bounded by ck with a constant c proportional to σ˙ L1 (0,T ;Q) + ˙vL1 (0,T ;V ) + u In the following we will allow c to be this kind of solution dependent constant, allowed by the solution regularities (3.13)–(3.15). Hence, θ G ,n 2Q  ck 2 + c

n−1   2      hk 2 hk 2    k σ j − σ hk j Q + vj − vj V + uj − uj V . j =0

We can apply (3.17) and (3.18) to obtain n−1  2       hk 2   θ G ,n 2Q  c h2 + k 2 + c k σ j − σ hk j Q + vj − vj V .

(3.23)

j =1

We subtract (3.11) from (3.16) at t = tn with w = wh ∈ V h to obtain        h   hk h h h v˙ n − δvhk = − σ n − σ hk + j βn−1 , uhk n ,w n ,ε w n−1 , w − j βn , un , w H Q From (3.10) and (2.24) at t = tn we obtain     h     h   hk = Aε vn − vhk σ n − σ hk n ,ε w n + Eε un − un−1 + θ G ,n , ε w Q Q

∀wh ∈ V h .

hk Write σ n − σ hk n = (I − PQh )σ n + PQh σ n − σ n . Note that (I − PQh )σ n Q  ch, and       hk hk PQh σ n − σ hk n = PQh Aε vn − vn + Eε un − un−1 + θ G ,n .

Use the non-expansiveness of the projection operator,          P h σ n − σ hk   c vn − vhk  + un − uhk  + θ G ,n  . Q n Q n V n−1 V Q

∀wh ∈ V h .

(3.24)

(3.25)

W. Han, M. Sofonea / Applied Numerical Mathematics 57 (2007) 498–509

Then,

507

        σ n − σ hk 2  c h2 + vn − vhk 2 + un − uhk 2 + θ G ,n 2 . n Q n V Q n−1 V

So applying (3.19) and (3.23), we have n−1       2      hk 2  σ n − σ hk 2  c h2 + k 2 + vn − vhk 2 + c  k σ j − σ hk n Q n V j Q + vj − vj V . j =1

Applying Lemma 3.2, we obtain n−1        2  σ n − σ hk 2  c h2 + k 2 + vn − vhk 2 + c  k vj − vhk n Q n V j V.

(3.26)

j =1

hk hk hk Consider the quantity An = ((δ(vn − vhk n ), vn − vn ))H + (Aε(vn − vn ), ε(vn − vn ))Q . Easily, we have the following lower bound:        1  vn − vhk 2 − vn−1 − vhk 2 + cvn − vhk 2 . (3.27) An  n H n V n−1 H 2k For an upper bound, we write An as            h hk h ˙ n , whn − vhk δ vn − vhk n , vn − wn H + Aε vn − vn , ε vn − wn Q + δvn − v n H    h     h hk hk hk , ε w + v˙ n − δvhk , w − v + Aε v − v − v . n n n n n n n H Q

Apply (3.24) and (3.25),    h     h hk hk + Aε vn − vhk v˙ n − δvhk n , w n − vn n , ε w n − vn H Q     hk     h   hk h hk hk h hk + θ . = − Eε un − uhk , ε w − v − j β , u n n , wn − vn − j βn−1 , un−1 , wn − vn G ,n n n n−1 Q From the definition of j and the assumptions, we have            j βn , un , vhk − wh − j β hk , uhk , vhk − wh   c βn − β hk  + un − uhk  vhk − wh  . n n n n n V n−1 n−1 n n−1 B n−1 V So we can bound An as follows:        hk   h hk    An  c δvn − v˙ n H + un − uhk n−1 V + θ G ,n Q + βn − βn−1 B wn − vn V         h hk   h  + δ vn − vhk n , vn − w n H + c vn − vn V vn − w n V .

(3.28)

h h hk We then combine (3.27) and (3.28), bound vhk n − wn V by vn − wn V + vn − vn V , and use the inequality 2 2 ab  a + b /(4 ) ∀ > 0 to obtain

       1  vn − vhk 2 − vn−1 − vhk 2 + cvn − vhk 2 n H n V n−1 H 2k   2 2    2 hk 2    c vn − whn V + δvn − v˙ n 2H + un − uhk n−1 V + θ G ,n Q + βn − βn−1 B     h + δ vn − vhk n , vn − w n H . Using (3.19), (3.21) and (3.23), we have        1  vn − vhk 2 − vn−1 − vhk 2 + vn − vhk 2 n H n V n−1 H k 2       h 2 2  c vn − whn V + δvn − v˙ n 2H + δ vn − vhk n , vn − w n H + h + k n−1   2      hk 2 hk 2    k σ j − σ hk +c j Q + vj − vj V + βj − βj B . j =1

Multiply this inequality by k, change n to j , and sum over j = 1 to n,

508

W. Han, M. Sofonea / Applied Numerical Mathematics 57 (2007) 498–509 n       2 vn − vhk 2 − v0 − vh 2 +  k vj − vhk n H j V 0 H j =1

c

n  j =1

 2 k vj − whj V + c

n 

kδvj − v˙ j 2H + c

j =1

n      h k δ vj − vhk j , vj − w j H j =1

j n−1    2         + vi − vhk 2 + βi − β hk 2 . + c h2 + k 2 + c k k σ i − σ hk i i i Q V B j =1

(3.29)

i=1

t  Note that δvj − v˙ j H  tnn+1 ¨v(t)H dt, and so nj=1 kδvj − v˙ j 2H  ck 2 due to the regularity v ∈ H 2 (0, T ; H ). From now on, we take whj = Π h vj , 1  j  N . Then     vj − wh  + hvj − wh   ch2 , 1  j  N. j H j V Using these bounds as well as (3.26) and (3.22) in (3.29), we have n n             h vn − vhk 2 + vj − vhk 2  c h2 + k 2 + c k k δ vj − vhk n H j V j , vj − w j H j =1

j =1

+c

j n−1    2  . k k vi − vhk i V j =1

(3.30)

i=1

Write n          h hk h h h k δ vj − vhk j , vj − w j H = vn − vn , vn − w n H − v0 − v0 , v1 − w 1 H j =1

n−1       h h + vj − vhk j , vj − wj − vj +1 − wj +1 H j =1

        h  h   h      vn − vhk n H vn − w n H + v0 − v0 H v1 − w 1 H +

n−1        vj − vhk   vj − wh − vj +1 − wh  . j j j +1 H H j =1

For j = 1, . . . , N − 1,        vj − wh − vj +1 − wh  = (vj − vj +1 ) − Π h (vj − vj +1 ) j j +1 H H  chvj − vj +1 V  ch˙vL1 (tj ,tj +1 ;V ) . Hence, c

n n−1     2      1 h  . vn − vhk 2 + ch2 + c , v k δ vj − vhk − w  k vj − vhk j n H j j H j H 2 j =1

j =1

Using this in (3.30), we have   j n n−1              2 hk 2 2 hk 2 hk 2     vn − vhk 2 + k vj − vj V  c h + k + c k vj − vj H + k vi − vi V . n H j =1

Applying Lemma 3.2, we obtain

j =1

i=1

W. Han, M. Sofonea / Applied Numerical Mathematics 57 (2007) 498–509

 max

1nN

 max

1nN

n    2   vn − vhk 2 + k vj − vhk n H j V



j =1



n    2   vn − vhk  + k vj − vhk n H j V

   c h2 + k 2 ,

1/2   c(h + k).

(3.31)

j =1

Using (3.18), (3.26) and (3.22), we obtain  N 1/2     2   hk hk max un − un V + k σ n − σ n Q + max βn − βnhk B  c(h + k). 1nN

509

n=1

1nN

(3.32)

Summarizing the estimates above establishes the following theorem. Theorem 3.3. Assume the conditions stated in Theorem 2.3, and the solution regularities (3.13), (3.14), (3.15). Then the optimal order error estimates (3.31) and (3.32) hold. In the error estimation, we need to make additional solution regularity assumptions. It is possible to show the convergence of the numerical solution to the true solution, without assuming the additional solution regularity assumptions. Such a convergence analysis is deemed to be delicate and lengthy. Due to the space limitation, we do not provide such a convergence analysis in this paper, but refer the interested reader to [10, Section 11.4], [4], and [12, Section 17.6] for basic ideas of the convergence analysis. References [1] A. Amassad, C. Fabre, Existence for viscoplastic contact with Coulomb friction problems, Int. J. Math. Math. Sci. 32 (2002) 411–437. [2] A. Amassad, C. Fabre, M. Sofonea, A quasistatic viscoplastic contact problem with normal compliance and friction, IMA J. Appl. Math. 69 (2004) 463–482. [3] V. Barbu, Nonlinear Semigroups and Differential Equations in Banach Spaces, Editura Academiei, Bucharest, Noordhoff, Leyden, 1976. [4] J. Chen, W. Han, M. Sofonea, Numerical analysis of a class of evolution systems with applications in viscoplasticity, SIAM J. Numer. Anal. 38 (2000) 1171–1199. [5] N. Cristescu, I. Suliciu, Viscoplasticity, Martinus Nijhoff Publishers, Editura Tehnica, Bucharest, 1982. [6] C. Eck, J. Jarušek, M. Krbeˇc, Unilateral Contact Problems: Variational Methods and Existence Theorems, Pure and Applied Mathematics, vol. 270, Chapman/CRC Press, New York, 2005. [7] J.R. Fernández-García, M. Sofonea, J.M. Viaño, A frictionless contact problem for elastic-viscoplastic materials with normal compliance, Numer. Math. 90 (2002) 689–719. [8] M. Frémond, Adhérence des solides, J. Mécanique Théorique et Appliquée 6 (1987) 383–407. [9] M. Frémond, Non-Smooth Thermomechanics, Springer, Berlin, 2002. [10] W. Han, B.D. Reddy, Plasticity: Mathematical Theory and Numerical Analysis, Springer, Berlin, 1999. [11] W. Han, M. Sofonea, Evolutionary variational inequalities arising in viscoelastic contact problems, SIAM J. Numer. Anal. 38 (2000) 556–579. [12] W. Han, M. Sofonea, Quasistatic Contact Problems in Viscoelasticity and Viscoplasticity, Studies in Advanced Mathematics, vol. 30, American Mathematical Society–International Press, Providence, RI, 2002. [13] I.R. Ionescu, M. Sofonea, Functional and Numerical Methods in Viscoplasticity, Oxford University Press, Oxford, 1994. [14] J. Jarušek, Contact problem with given time-dependent friction force in linear viscoelasticity, Comment. Math. Univ. Carolinae 31 (1990) 257–262. [15] J. Jarušek, Dynamic contact problems with given friction for viscoelastic bodies, Czechoslovak Math. J. 46 (1996) 475–487. [16] J. Jarušek, C. Eck, Dynamic contact problems with small Coulomb friction for viscoelastic bodies. Existence of solutions, Math. Models Methods Appl. Sci. 9 (1999) 11–34. [17] N. Kikuchi, J.T. Oden, Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods, SIAM, Philadelphia, PA, 1988. [18] J.A.C. Martins, J.T. Oden, Existence and uniqueness results for dynamic contact problems with nonlinear normal and friction interface laws, Nonlinear Anal. TMA 11 (1987) 407–428. [19] M. Raous, L. Cangémi, M. Cocu, A consistent model coupling adhesion, friction and unilateral contact, Comput. Methods Appl. Engrg. 177 (1999) 383–399. [20] M. Rochdi, M. Shillor, M. Sofonea, A quasistatic viscoelastic contact problem with normal compliance and friction, J. Elasticity 51 (1998) 105–126. [21] M. Shillor, M. Sofonea, J.J. Telega, Models and Analysis of Quasistatic Contact, Lecture Notes in Physics, vol. 655, Springer, Berlin, 2004. [22] M. Sofonea, W. Han, M. Shillor, Analysis and Approximation of Contact Problems with Adhesion or Damage, Pure and Applied Mathematics, vol. 276, Chapman-Hall/CRC Press, New York, 2006.