Numerical analysis of a frictionless contact problem for elastic–viscoplastic materials

Comput. Methods Appl. Mech. Engrg. 190 (2000) 179±191

www.elsevier.com/locate/cma

Numerical analysis of a frictionless contact problem for elastic±viscoplastic materials Weimin Han a,*, Mircea Sofonea b a

Department of Mathematics, University of Iowa, Iowa City, IA 52242-1410, USA b Laboratoire de Th eorie des Syst emes, University of Perpignan, France Received 16 September 1998

Abstract We consider a mathematical model which describes the unilateral quasistatic contact of two elastic±viscoplastic bodies. The contact is without friction and it is modeled by the classical Signorini boundary conditions. The model consists of an evolution equation coupled with a time-dependent variational inequality. It has been shown that the variational problem of the model has a unique solution. Here we consider numerical approximations of the problem. We use the ®nite element method to discretize the spatial domain. Spatially semi-discrete and fully discrete schemes are studied. For both schemes, we show the existence of a unique solution, and derive error estimates. Under appropriate regularity assumptions of the solution, we have the optimal order convergence. Ó 2000 Elsevier Science S.A. All rights reserved. Keywords: Quasistatic frictionless contact problem; Elastic±viscoplastic material; Time-dependent variational ineqaulity; Semi-discrete approximation; Fully discrete approximation; Finite element method; Convergence; Error estimate

1. Introduction In this paper we give numerical analysis of a mathematical model which describes the unilateral quasistatic contact of two elastic±viscoplastic bodies. The contact is without friction and it is modeled by the classical Signorini boundary conditions. Frictionless contact between two elastic bodies is a topic considered by many authors; see, e.g., [10±12,14] on well-posedness and numerical approximations of these contact problems. Frictionless contact problems involving elastic perfectly plastic bodies are considered in, e.g., [2,3,9]. In this paper we consider frictionless contact problems of elastic±viscoplastic bodies modeled by a constitutive law of the form r_ ˆ E_e ‡ G…r; e†;

…1:1†

in which r denotes the stress tensor, e is the small strain tensor and E, G are the given constitutive functions. Here and below the dot above represents the derivative with respect to the time variable t. A concrete example of an elastic±viscoplastic constitutive law of this form in the case when G does not depend on e is the Perzyna law given by 1 e_ ˆ Eÿ1 r_ ‡ …r ÿ PK r†; l *

Corresponding author. Tel.: +1-319-335-0770; fax: +1-319-335-0627. E-mail address: [email protected] (W. Han).

0045-7825/00/$ - see front matter Ó 2000 Elsevier Science S.A. All rights reserved. PII: S 0 0 4 5 - 7 8 2 5 ( 9 9 ) 0 0 4 2 0 - X

180

W. Han, M. Sofonea / Comput. Methods Appl. Mech. Engrg. 190 (2000) 179±191

in which l > 0 is a viscosity constant, K is a nonempty, closed, convex set in the space of symmetric tensors and PK represents the projection mapping. A relatively simple one-dimensional example of constitutive law of the form (1.1) in the case when a full coupling in stress and strain is involved in G is obtained by taking 8 < ÿk1 F1 …r ÿ f …e†† if r > f …e†; if g…e† 6 r 6 f …e†; G…r; e† ˆ 0 : if r < g…e†; k2 F2 …g…e† ÿ r† where k1 ; k2 > 0 are viscosity constants, f and g are Lipschitz continuous functions with g…e† 6 f …e†, and F1 ; F2 : R‡ ! R are increasing functions with F1 …0† ˆ F2 …0† ˆ 0 (see [6, p. 35]). Rate-type viscoplastic models of the form (1.1) are considered in order to describe the behavior of real materials like rubbers, metals, pastes, rocks and so on. Various results and mechanical interpretation concerning models of this form may be found for instance in [6] (see also the references quoted there). Existence and uniqueness results for initial and boundary value problems involving (1.1) were obtained in [13] in the case of displacement-tractions conditions and in [1,17] in the case of contact conditions with/ without friction. The mathematical model describing the unilateral contact without friction between two elastic±viscoplastic bodies, each having a constitutive law of the form (1.1) was analyzed in [16]. Assuming the quasistatic case, this model consists of an evolution equation coupled with a time-dependent variational inequality. The well-posedness of the problem has been proved in [16] using standard arguments of elliptic variational inequalities and a ®xed point property. The aim of this paper is to consider numerical approximations of the mechanical problem analyzed in [16]. In Section 2 we state the contact problem and its variational formulation. We also recall the assumptions on the problem data together with the main existence and uniqueness result (Theorem 2.1). In Sections 3 and 4, we analyze spatially semi-discrete and fully discrete schemes, respectively. We use the ®nite element method to discretize the spatial domain. For both schemes, we show the existence of a unique solution. We also show the convergence of the schemes and derive error estimates. Under appropriate regularity assumptions of the solution, optimal order error estimates are obtained. We remark that it is also possible to analyze another semi-discrete scheme, namely the temporally semi-discrete scheme, following the arguments presented in [7] for the analysis of temporally semi-discrete schemes in solving elastoplasticity problems. 2. The contact problem We consider two elastic±viscoplastic bodies whose material particles occupy bounded domains X1 and X of Rd (d 6 3 in applications). We put a superscript m to indicate that the quantity is related to the domain Xm ; m ˆ 1; 2. In the following, the superscript m ranges between 1 and 2. For each domain Xm , we assume its boundary Cm is Lipschitz continuous, and is partitioned into three disjoint measurable parts Cm1 , Cm2 and Cm3 , with meas…Cm1 † > 0. The unit outward normal to Cm , de®ned almost everywhere, is denoted by m m ˆ …mmi †. We are interested in the resulting contact process of evolution of the bodies on the time interval ‰0; T Š; T > 0. Assume that the bodies are clamped on Cm1  …0; T † and so the displacement ®elds vanish there. We also assume that volume forces of density um1 act on Xm  …0; T † and that surface tractions of density um2 act on Cm2  …0; T †. The two bodies are in contact along the common part C13 ˆ C23 , which will be denoted C3 below. The contact is frictionless and it is modelled by the Signorini's conditions on Cm3  …0; T † in the form with a zero gap function. Finally we assume the process is quasistatic and we shall use (1.1) as constitutive law. With these assumptions, the mechanical problem we study here may be formulated as follows: For m ˆ 1; 2, ®nd the displacement ®eld um ˆ …umi † : Xm  ‰0; T Š ! Rd and the stress ®eld rm ˆ …rmij † : m X  ‰0; T Š ! Sd which satisfy, for t 2 …0; T †, 2

r_ m ˆ Em e…u_ m † ‡ Gm …rm ; e…um †† in Xm ; m

Divr ‡

um1

ˆ0

m

in X ;

…2:1† …2:2†

W. Han, M. Sofonea / Comput. Methods Appl. Mech. Engrg. 190 (2000) 179±191

um ˆ 0 m m

on Cm1 ;

r m ˆ um2

u1m ‡ u2m 6 0;

181

…2:3†

on Cm2 ;

…2:4†

r1m ˆ r2m 6 0;

r1m …u1m ‡ u2m † ˆ 0;

rms ˆ 0 on C3 ;

…2:5†

and the initial value conditions um …0† ˆ um0 ;

in Xm :

rm …0† ˆ rm0

…2:6†

Here and below, i; j ˆ 1; 2; . . . ; d, Sd represents the space of second-order symmetric tensors on Rd (or equivalently, the space of the symmetric matrices of order d). Throughout the paper, the summation convention is adopted over a repeated index, unless stated otherwise. In (2.1)±(2.5), Div rm represents the vector-valued ®eld of the divergence of the tensor-valued function m r given by Div rm ˆ …rmij;j †; e…um † denotes the small strain tensor given by 1 eij …um † ˆ …umi;j ‡ umj;i † 2

e…um † ˆ …eij …um ††;

and the index that follows a comma indicates a partial derivative with respect to the corresponding component of the spatial variable. Moreover, we used the notation umm ˆ umi mmi ;

rmm ˆ rmij mmi mmj ;

rms ˆ …rmsi † with rmsi ˆ rmij mmj ÿ rmm mmi

for the normal displacement, normal stress and tangential stress, respectively. As usual, the inner products and the corresponding norms on Rd and Sd are u
v ˆ u i vi ;

jvj ˆ …v
v†1=2

r
s ˆ rij sij ;

jsj ˆ …s
s†

1=2

8u; v 2 Rd ; 8r; s 2 Sd :

We introduce the spaces V m ˆ fv ˆ …vi † j vi 2 H 1 …Xm †; vi ˆ 0 on Cm1 ; 1 6 i 6 dg; Hm ˆ fs ˆ …sij † j sij 2 L2 …Xm †; 1 6 i; j 6 dg; Hm1 ˆ fs 2 Hm j Div s 2 …L2 …Xm ††d g: These are Hilbert spaces with their canonical inner products. We will also use vector-valued function spaces, such as W 1;1 …0; T ; X † for a Sobolev space X; a detailed discussion of these spaces can be found in [18]. Since meas…Cm1 † > 0, Korn's inequality holds ke…v†kHm P ckvk…H 1 …Xm ††d

8v 2 V m :

…2:7†

Here c > 0 is a constant which depends only on Xm and Cm1 (for a proof, see for instance [15, p. 79]). Therefore, the map v7!ke…v†kHm de®nes a norm on V m , which is equivalent to the norm kvk…H 1 …Xm ††d . In the study of the mechanical problem (2.1)±(2.6) we make the following assumptions with m ˆ 1; 2. Em : Xm  Sd ! Sd is a bounded, symmetric, positive de®nite fourth-order tensor: …a† Emijkl 2 L1 …Xm †;

1 6 i; j; k; l 6 d:

…b† E r
s ˆ r
E s 8r; s 2 Sd ; a:e: in Xm : …c† There exists an am > 0 such that m

m

Em s
s P am jsj

2

8s 2 Sd ; a:e: in Xm :

…2:8†

182

W. Han, M. Sofonea / Comput. Methods Appl. Mech. Engrg. 190 (2000) 179±191

Gm : Xm  Sd  Sd ! Sd has the properties: …a† There exists an Lm > 0 such that 8r1 ; r2 ; e1 ; e2 2 Sd ; a:e: in Xm ; jGm …x; r1 ; e1 † ÿ Gm …x; r2 ; e2 †j 6 Lm …jr1 ÿ r2 j ‡ je1 ÿ e2 j†: …b† For any r; e 2 Sd ; x 7! Gm …x; r; e† is measurable: …c† The mapping x 7! Gm …x; 0; 0† 2 Hm :

…2:9†

For the given force densities, we assume um1 2 W 1;1 …0; T ; …L2 …Xm ††d †;

um2 2 W 1;1 …0; T ; …L2 …Cm2 ††d †:

…2:10†

We then de®ne product spaces V ˆ V 1  V 2 ; H ˆ H1  H2 and H1 ˆ H11  H21 . They are all Hilbert spaces endowed with the canonical inner products denoted by …
;
†V , …
;
†H and …
;
†H1 , respectively. The associated norms will be denoted by k
kV ; k
kH and k
kH1 , respectively. Moreover, we denote by …V 0 ; k
kV 0 † the strong dual of V and h
;
i will represent the duality between V 0 and V. For all t 2 ‰0; T Š let f…t† denote the element of V 0 given by hf…t†; vi ˆ …u11 …t†; v1 †…L2 …X1 ††d ‡ …u21 …t†; v2 †…L2 …X2 ††d ‡ …u12 …t†; v1 †…L2 …C1 ††d ‡ …u22 …t†; v2 †…L2 …C2 ††d 2

1

2

2

8v ˆ …v ; v † 2 V : Using (2.10), we see that f 2 W 1;1 …0; T ; V 0 †: We de®ne the set U of admissible displacement ®elds by U ˆ fv ˆ …v1 ; v2 † 2 V j v1m ‡ v2m 6 0 on C3 g;

…2:11†

and we suppose that u0 ˆ …u10 ; u20 † 2 U ;

…r0 ; e…u0 ††H ˆ hf…0†; u0 i:

…2:12†

Finally, for simplicity we shall use the notation e…v† ˆ …e…v1 †; e…v2 †† for all v ˆ …v1 ; v2 † 2 V and Ee ˆ …E1 e1 ; E2 e2 †; G…r; e† ˆ …G…r1 ; e1 †; G…r2 ; e2 †† for all e ˆ …e1 ; e2 † 2 H and all r ˆ …r1 ; r2 † 2 H. With these notations, in [16] the following weak formulation of the contact problem was derived. Problem P. Find the displacement ®eld u : ‰0; T Š ! U and the stress ®eld r : ‰0; T Š ! H1 such that u…0† ˆ u0 ;

r…0† ˆ r0 ;

…2:13†

and for a.a. t 2 …0; T †, _ _ ˆ Ee…u…t†† r…t† ‡ G…r…t†; e…u…t†††; …r…t†; e…v ÿ u…t†††H P hf…t†; v ÿ u…t†i

8v 2 U :

…2:14† …2:15†

The well-posedness of the Problem P has been studied in [16] and the main existence and uniqueness result is the following. Theorem 2.1. Under the assumptions (2.8)±(2.10) and (2.12), Problem P has a unique solution …u; r† 2 W 1;1 …0; T ; U  H1 †. Remark 2.2. A careful examination of the arguments in [16] shows that if the space W 1;1 in the condition (2.10) is replaced by W 1;p ; 1 6 p < 1, then Theorem 2.1 is still valid with the conclusion that Problem P has a unique solution …u; r† 2 W 1;p …0; T ; U  H1 †.

W. Han, M. Sofonea / Comput. Methods Appl. Mech. Engrg. 190 (2000) 179±191

183

3. Spatially semi-discrete approximation We will discretize the spatial domain by the ®nite element method. Let Th be a regular ®nite element partition of the domain X in such a way that if a side of an element lies on the boundary, then the side is m m entirely on one of the subsets C1 ; C2 and C3 . We then choose a ®nite element space V h  V for the approximation of the displacement variable u, and another ®nite element space Qh for the approximation of the stress variable r. The ®nite element spaces consist of piecewise (images of) polynomials. At our general level of discussions below, we do not need to be speci®c about the detailed structures of the ®nite element spaces. We will only assume that for these ®nite element spaces, there hold the relations e…V h †  Qh ;

G…Qh ; Qh †  Qh :

…3:1†

Such a requirement is satis®ed if, e.g., we use linear elements for V h and piecewise constants for Qh . Then we de®ne the discrete admissible set 2;h U h ˆ fvh ˆ …v1;h ; v2;h † 2 V h j v1;h m ‡ vm 6 0 on C3 g:

Obviously, U h  U . Now the spatially semi-discrete scheme is the following. Problem Ph . Find the displacement ®eld uh : ‰0; T Š ! U h and the stress ®eld rh : ‰0; T Š ! Qh such that uh …0† ˆ uh0 ;

rh …0† ˆ rh0 ;

…3:2†

and for a.a. t 2 …0; T †, r_ h …t† ˆ Ee…u_ h …t†† ‡ G…rh …t†; e…uh …t†††; h

h

h

h

h

…r …t†; e…v ÿ u …t†††H P hf…t†; v ÿ u …t†i

…3:3† h

h

8v 2 U :

…3:4†

Here, uh0 2 U h ; rh0 2 Qh are appropriate approximations of u0 and r0 . Remark 3.1. The assumption G…Qh ; Qh †  Qh in (3.1) looks unnatural for higher order elements. In [4], we consider numerical approximations of a family of similar contact problems involving elastic±viscoplastic materials; we demonstrate that it is possible to develop and analyze numerical schemes without the assumption G…Qh ; Qh †  Qh . Using the arguments presented in [16], it can be shown that Problem Ph has a unique solution …u ; rh † 2 W 1;1 …0; T ; U h  Qh †. In this section, we are interested in deriving estimates for the errors u ÿ uh and r ÿ rh . h

Theorem 3.2. Let …u; r† 2 W 1;1 …0; T ; U  H1 † be the solution of Problem P, and …uh ; rh † 2 W 1;1 …0; T ; U h  Qh † be the solution of the semi-discrete Problem Ph . Assume (3.1) holds. Then we have the error estimate kr ÿ rh kL1 …0;T ;H† ‡ ku ÿ uh kL1 …0;T ;V † 6 c…kr0 ÿ rh0 kH ‡ ku0 ÿ uh0 kV † n o 1=2 1=2 ku ÿ vh kL1 …0;T ;V † ‡ …kfkL1 …0;T ;V 0 † ‡ krkL1 …0;T ;H† † ku ÿ vh kL1 …0;T ;V † : ‡c inf vh 2L1 …0;T ;U h †

…3:5†

Proof. First we integrate (2.14) and (3.3), and use the initial conditions (2.13) and (3.2) to obtain Z r…t† ˆ Ee…u…t†† ‡ rh …t† ˆ Ee…uh …t†† ‡

t 0

G…r…s†; e…u…s††† ds ‡ r0 ÿ Ee…u0 †;

Z 0

t

G…rh …s†; e…uh …s††† ds ‡ rh0 ÿ Ee…uh0 †:

…3:6† …3:7†

184

W. Han, M. Sofonea / Comput. Methods Appl. Mech. Engrg. 190 (2000) 179±191

We then subtract (3.7) from (3.6) to get the ®rst error relation Z t   G…r…s†; e…u…s††† ÿ G…rh …s†; e…uh …s††† ds ‡ r0 ÿ rh0 ÿ Ee…u0 ÿ uh0 †: r…t† ÿ rh …t† ˆ Ee…u…t† ÿ uh …t†† ‡ 0

…3:8† Now we take v ˆ uh …t† in (2.15), …r…t†; e…uh …t† ÿ u…t†††H P hf…t†; uh …t† ÿ u…t†i; and rewrite (3.4) in the form …rh …t†; e…u…t† ÿ uh …t†††H P hf…t†; vh ÿ uh …t†i ‡ …rh …t†; e…u…t† ÿ vh ††H : Adding these two inequalities, we obtain the second error relation, …r…t† ÿ rh …t†; e…u…t† ÿ uh …t†††H 6 …r…t† ÿ rh …t†; e…u…t† ÿ vh ††H ‡ R…u…t† ÿ vh †;

…3:9†

where …3:10†

R…v† ˆ hf…t†; vi ÿ …r…t†; e…v††H represents the residual. Denote e0 ˆ kr0 ÿ rh0 kH ‡ ku0 ÿ uh0 kV :

…3:11†

We then derive the following inequality from (3.8) with the use of the assumptions (2.8) and (2.9), Z t ÿ
h h kr…s† ÿ rh …s†kH ‡ ke…u…s† ÿ uh …s††kH ds ‡ ce0 : kr…t† ÿ r …t†kH 6 cke…u…t† ÿ u …t††kH ‡ c 0

…3:12† Plug (3.8) into (3.9) and again use the assumptions (2.8) and (2.9). After some elementary manipulations, we have 2

ke…u…t† ÿ uh …t††kH 6 ckr…t† ÿ rh …t†kH ke…u…t† ÿ vh †kH ‡ cjR…u…t† ÿ vh †j ‡ cke…u…t† ÿ uh …t††kH e0 Z t ÿ
‡ cke…u…t† ÿ uh …t††kH kr…s† ÿ rh …s†kH ‡ ke…u…s† ÿ uh …s††kH ds: 0

Then we use the inequality ab 6 da2 ‡

1 2 b 4d

8d > 0

to obtain ke…u…t† ÿ uh …t††kH 6 d0 kr…t† ÿ rh …t†kH ‡ cke…u…t† ÿ vh †kH ‡ ce0 ‡ cjR…u…t† ÿ vh †j1=2 Z t
ÿ ‡c kr…s† ÿ rh …s†kH ‡ ke…u…s† ÿ uh …s††kH ds;

…3:13†

0

where d0 > 0 is arbitrary. Let d0 be a suciently small positive number. Then combining (3.12) and (3.13), we obtain, for any vh ˆ vh …t† 2 U h , 1=2

kr…t† ÿ rh …t†kH ‡ ke…u…t† ÿ uh …t††kH 6 cke…u…t† ÿ vh …t††kH ‡ ce0 ‡ cjR…u…t† ÿ vh …t††j Z t
ÿ ‡c kr…s† ÿ rh …s†kH ‡ ke…u…s† ÿ uh …s††kH ds; 0

W. Han, M. Sofonea / Comput. Methods Appl. Mech. Engrg. 190 (2000) 179±191

185

or equivalently 1=2

kr…t† ÿ rh …t†kH ‡ ku…t† ÿ uh …t†kV 6 cku…t† ÿ vh …t†kV ‡ ce0 ‡ cjR…u…t† ÿ vh …t††j Z t ÿ
‡c kr…s† ÿ rh …s†kH ‡ ku…s† ÿ uh …s†kV ds 0

…3:14†

A direct estimate on the residual R…
† de®ned in (3.10) yields jR…u…t† ÿ vh …t††j 6 c…kf…t†kV 0 ‡ kr…t†kH †ku…t† ÿ vh …t†kV :

…3:15†

Since vh …t† 2 U h is arbitrary, applying Gronwall's inequality, we obtain the estimate (3.5) from (3.14).



The estimate (3.5) is the basis for convergence analysis. Indeed, arguing as in [8], we immediately have the following consequence. Corollary 3.3. Assume the initial values rh0 2 Qh and uh0 2 U h are chosen in such a way that kr0 ÿ rh0 kH ! 0;

ku0 ÿ uh0 kV ! 0

as h ! 0:

d

Also assume, the set …H 2 …X†† \ U is dense in U and for some constant c and a > 0, inf kv ÿ vh kV 6 ckvk…H 2 …X††d ha

vh 2U h

d

8v 2 …H 2 …X†† \ U :

…3:16†

Then the semi-discrete method converges: kr ÿ rh kL1 …0;T ;H† ‡ ku ÿ uh kL1 …0;T ;V † ! 0

as h ! 0:

We comment that the assumption (3.16) is satis®ed as long as the ®nite element space V h contains piecewise linear functions, and we may take a ˆ 1. Thus the assumption (3.16) is not a restriction in applications. A result on the density of C 1 …X†d \ U in U is proved in [11] in the planar case d ˆ 2. Under suitable regularity assumptions, the estimate (3.5) can be improved. We assume r1m 2 L1 …0; T ; L2 …C3 ††:

…3:17†

Then it can be veri®ed that (cf. (3.12) in [16]) Z r1m …t†…v1m ‡ v2m † ds: …r…t†; e…v ÿ u…t†††H ˆ hf…t†; v ÿ u…t†i ‡ C3

Upon using (2.5), we see that Z 1 h2 2 r1m …t†…vh1 R…u…t† ÿ vh …t†† ˆ m …t† ÿ um …t† ‡ vm …t† ÿ um …t†† ds; C3

and therefore, instead of (3.15) jR…u…t† ÿ vh …t††j 6 kr1m …t†kL2 …C3 † kvhm …t† ÿ um …t†kL2 …C3 † ; where 1 h2 2 kvhm …t† ÿ um …t†kL2 …C3 † ˆ kvh1 m …t† ÿ um …t†kL2 …C3 † ‡ kvm …t† ÿ um …t†kL2 …C3 † :

we have the next result.

…3:18†

186

W. Han, M. Sofonea / Comput. Methods Appl. Mech. Engrg. 190 (2000) 179±191

Corollary 3.4. Assume additionally the regularity condition (3.17). Then kr ÿ rh kL1 …0;T ;H† ‡ ku ÿ uh kL1 …0;T ;V † 6 c…kr0 ÿ rh0 kH ‡ ku0 ÿ uh0 kV † o n 1=2 1=2 ku ÿ vh kL1 …0;T ;V † ‡ kr1m kL1 …0;T ;L2 …C3 †† kum ÿ vhm kL1 …0;T ;L2 …C3 †† : ‡ inf

…3:19†

vh 2L1 …0;T ;U h †

The inequality (3.19) is the basis for deriving order error estimates for the semi-discrete solutions. We present one sample result. Assume u 2 L1 …0; T ; …H 2 …X††d † and um 2 L1 …0; T ; H 2 …C3 ††. Let us use linear elements for the ®nite element space V h , and piecewise constants for Qh . We choose the initial values uh0 2 U h and rh0 2 Qh in such a way that kr0 ÿ rh0 kH 6 ch; ku0 ÿ uh0 kV 6 ch: Let Ph u…t† 2 V h be the piecewise linear interpolant of u…t†. Then (cf. [5]) ku…t† ÿ Ph u…t†kV 6 chju…t†j…H 2 …X††d ; kum …t† ÿ Ph um …t†kL2 …C3 † 6 ch2 jum …t†jH 2 …C3 † : Hence, inf

vh 2L1 …0;T ;U h †

n

1=2

1=2

ku ÿ vh kL1 …0;T ;V † ‡ kr1m kL1 …0;T ;L2 …C3 †† kum ÿ vhm kL1 …0;T ;L2 …C3 †† 1=2

o

1=2

6 ku ÿ Ph ukL1 …0;T ;V † ‡ kr1m kL1 …0;T ;L2 …C3 †† kum ÿ Ph um kL1 …0;T ;L2 …C3 ††   1=2 1=2 6 ch jujL1 …0;T ;…H 2 …X††d † ‡ kr1m kL1 …0;T ;L2 …C3 †† jum jL1 …0;T ;H 1 …C3 †† : Summarizing, under the regularity assumptions d

u 2 L1 …0; T ; …H 2 …X†† †;

um 2 L1 …0; T ; H 2 …C3 †† and

r1m 2 L1 …0; T ; L2 …C3 ††;

we have the optimal order error estimate kr ÿ rh kL1 …0;T ;H† ‡ ku ÿ uh kL1 …0;T ;V † 6 ch:

4. Fully discrete approximation Now we consider a fully discrete approximation of Problem P. We use the ®nite element spaces V h and Q introduced in the last section, with the assumption (3.1). Remark 3.1 concerning the restriction (3.1) still applies and it is possible to develop and analyze fully discrete schemes without this restriction. In addition, we need a partition of the time interval ‰0; T Š : 0 ˆ t0 < t1 <


< tN ˆ T . We denote the step-size kn ˆ tn ÿ tnÿ1 for n ˆ 1; . . . ; N . We allow nonuniform partition of the time interval, and denote k ˆ maxn kn the maximal step-size. For a sequence fwn gNnˆ0 , we denote Dwn ˆ wn ÿ wnÿ1 for the di€erence, and dwn ˆ Dwn =kn the corresponding divided di€erence. In this section, no summation is implied over the repeated index n. Then a fully discrete approximation method based on a backward Euler scheme is the following. h

N

N

h hk hk h Problem Phk . Find the displacement ®eld uhk ˆ fuhk n gnˆ0  U and the stress ®eld r ˆ frn gnˆ0  Q , such that h hk h uhk 0 ˆ u0 ; r0 ˆ r0 ;

…4:1†

W. Han, M. Sofonea / Comput. Methods Appl. Mech. Engrg. 190 (2000) 179±191

187

and for n ˆ 1; . . . ; N , hk hk hk drhk n ˆ Ede…un † ‡ G…rn ; e…un ††; h …rhk n ; e…v

ÿ

uhk n ††H

h

P hf n ; v ÿ

uhk n i

…4:2† h

h

8v 2 U :

…4:3†

Here the point value f n ˆ f…tn † is well-de®ned by the Sobolev embedding theorem H 1 …0; T ; X † ,! C…‰0; T Š; X † for any normed space X. In the error analysis below, we will similarly use the notation un ˆ u…tn † 2 V and rn ˆ r…tn † 2 H. Theorem 4.1. There exists an k0 > 0, such that if k 6 k0 , then Problem Phk has a unique solution …uhk ; rhk †. hk Proof. We prove the result by induction. For some n between 1 and N, assume uhk nÿ1 and rnÿ1 are known. hk h hk h Then the problem for un 2 U and rn 2 Q is (4.3) and hk hk hk hk hk rhk n ˆ Ee…un † ‡ kn G…rn ; e…un †† ‡ rnÿ1 ÿ Ee…unÿ1 †;

…4:4†

which is obtained from (4.2). Let us use the Banach ®xed-point theorem to prove the unique solvability of the problem (4.3) and (4.4). Given gh 2 Qh , we consider the problem of ®nding ugh 2 U h and rgh 2 Qh such that hk rgh ˆ Ee…ugh † ‡ kn gh ‡ rhk nÿ1 ÿ Ee…unÿ1 †; h

h

…rgh ; e…v ÿ ugh ††H P hf n ; v ÿ ugh i

h

…4:5† h

8v 2 U :

…4:6†

We can eliminate the variable rgh to get an elliptic variational inequality hk h …Ee…ugh †; e…vh ÿ ugh ††H P hf n ; vh ÿ ugh i ÿ …kn gh ‡ rhk nÿ1 ÿ Ee…unÿ1 †; e…v ÿ ugh ††H

8vh 2 U h ;

which has a unique solution ugh 2 U h . From (4.5), we then determine a unique rgh 2 Qh . Now let us de®ne an operator Kh : Qh ! Qh by the formula Kh gh ˆ G…rgh ; e…ugh ††; and show that Kh is a contraction. Thus let gh1 ; gh2 2 Qh , and denote the corresponding solutions of problem hk hk hk (4.5) and (4.6) by …uhk g1 ; rg1 † and …ug2 ; rg2 †. Then subtracting the corresponding relations (4.5), we get hk hk hk h h rhk g1 ÿ rg2 ˆ Ee…ug1 ÿ ug2 † ‡ kn …g1 ÿ g2 †:

…4:7†

Take the norm, hk hk hk h h krhk g1 ÿ rg2 kH 6 cke…ug1 ÿ ug2 †kH ‡ kn kg1 ÿ g2 kH :

…4:8†

h h h hk We then use the relation (4.6) for gh ˆ gh1 with vh ˆ uhk g2 and that for g ˆ g2 with v ˆ ug1 , and add the two inequalities, hk hk hk …rhk g1 ÿ rg2 ; Ee…ug1 ÿ ug2 †† 6 0:

This inequality, together with (4.7), implies hk hk hk h h hk hk …Ee…uhk g1 ÿ ug2 †; e…ug1 ÿ ug2 †† 6 ÿ kn …g1 ÿ g2 ; e…ug1 ÿ ug2 ††:

Hence, hk h h ke…uhk g1 ÿ ug2 †kH 6 ckn kg1 ÿ g2 kH :

…4:9†

Then from (4.8), hk h h krhk g1 ÿ rg2 kH 6 ckn kg1 ÿ g2 kH :

…4:10†

188

W. Han, M. Sofonea / Comput. Methods Appl. Mech. Engrg. 190 (2000) 179±191

Now

  hk hk hk hk hk hk hk ; e…u †† ÿ G…r ; e…u ††k 6 c ke…u †k ‡ kr ÿ r k ÿ u kKh gh1 ÿ Kh gh2 kH ˆ kG…rhk H H g1 g1 g2 g2 g1 g2 g1 g2 H :

By (4.9) and (4.10), kKh gh1 ÿ Kh gh2 kH 6 ckn kgh1 ÿ gh2 kH : Therefore, if kn 6 k is suciently small, Kh is a contraction on Qh . By the Banach ®xed-point theorem, Kh has a unique ®xed-point gh in Qh , the corresponding solution …ugh ; rgh † 2 U h  Qh is the unique solution hk h h …uhk  n ; rn † 2 U  Q of the problem (4.3) and (4.4). Now let us derive an error estimate for the fully discrete solution …uhk ; rhk †. Apply (4.4) recursively, hk rhk n ˆ Ee…un † ‡

n X jˆ1

hk hk hk kj G…rhk j ; e…uj †† ‡ r0 ÿ Ee…u0 †:

…4:11†

Subtracting (4.11) from (3.6) at t ˆ tn , we obtain the ®rst error relation n h i Z tn X hk hk hk † ˆ r ÿ r ÿ k G…r ; e…u †† ÿ G…r ; e…u †† ÿ G…r; e…u†† dt Ee…un ÿ uhk n j j j n n j j 0

jˆ1

‡

n X jˆ1

  hk kj G…rj ; e…uj †† ÿ r0 ÿ rhk 0 ÿ Ee…u0 ÿ u0 † :

…4:12†

Using (4.3) and (2.15) at t ˆ tn with v ˆ uhk n , we have hk h hk h …rn ÿ rhk n ; e…un ÿ un ††H 6 hf n ; un ÿ v i ÿ …rn ; e…un ÿ v ††H :

Thus we have the second error relation hk hk h h …rn ÿ rhk n ; e…un ÿ un ††H 6 …rn ÿ rn ; e…un ÿ v ††H ‡ Rn …un ÿ v †;

…4:13†

where Rn …v† ˆ hf n ; vi ÿ …rn ; e…v††:

…4:14†

Denote

Z tn n X G…r; e…u†† dt ÿ kj G…rj ; e…uj †† ; In ˆ 0 jˆ1

n ˆ 1; . . . ; N

…4:15†

for the errors of the numerical integration, and hk en ˆ krn ÿ rhk n kH ‡ kun ÿ un kV ;

n ˆ 0; . . . ; N

…4:16†

for the error of the numerical solution. From (4.12), we see that hk kun ÿ uhk n kV 6 ckrn ÿ rn kH ‡ cIn ‡ ce0 ‡ c

n X jˆ1

hk kj …krj ÿ rhk j kH ‡ kuj ÿ uj kV †:

Combining (4.12) and (4.13), we obtain ÿ1 hk hk h h hk …rn ÿ rhk n ; E …rn ÿ rn †† 6 …rn ÿ rn ; e…un ÿ v †† ‡ Rn …un ÿ v † ‡ ckrn ÿ rn kH " # n X hk  kj …krj ÿ rhk j kH ‡ kuj ÿ uj kV † ‡ In ‡ e0 ; jˆ1

…4:17†

W. Han, M. Sofonea / Comput. Methods Appl. Mech. Engrg. 190 (2000) 179±191

189

from which n X   h h 1=2 hk ÿ v ‡ c k 6 c ku k ‡ j R …u ÿ v † j ‡ I ‡ e kj …krj ÿ rhk krn ÿ rhk n n n n 0 V n H j kH ‡ kuj ÿ uj kV †: jˆ1

…4:18† Inequalities (4.17) and (4.18) can be used together to yield, for any vhn 2 V h , h i hk h h 1=2 k ‡ ku ÿ u k 6 c ku k ‡ jR …u ÿ v †j ‡ I ‡ e ÿ v krn ÿ rhk n n n n n 0 n H n V n V n n X

‡c

jˆ1

hk kj …krj ÿ rhk j kH ‡ kuj ÿ uj kV †:

…4:19†

We now derive an error estimate from (4.19). Theorem 4.2. Let …u; r† 2 W 1;1 …0; T ; U  H1 † be the solution of Problem P, and uhk  U h and rhk  Qh be the solution of the fully discrete Problem Phk . Assume (3.1) holds. Then we have the error estimate ÿ
hk max krn ÿ rhk n kH ‡ kun ÿ un kV 16n6N

1=2

6 c…kr0 ÿ rh0 kH ‡ ku0 ÿ uh0 kV † ‡ c max …kun ÿ vhn kV ‡ jRn …un ÿ vhn †j 16n6N   _ L1 …0;T ;V † : _ L1 …0;T ;H† ‡ kuk ‡ ck krk

† …4:20†

Proof. Denote n X

En ˆ

kj ej ;

…4:21†

jˆ1

and gn ˆ kun ÿ vhn kV ‡ j Rn …un ÿ vhn † j1=2 ‡In ‡ e0 :

…4:22†

Then (4.19) can be rewritten as en 6 cgn ‡ cEn ;

n ˆ 1; . . . ; N :

…4:23†

Now En ÿ Enÿ1 ˆ kn en 6 ckn gn ‡ ckn En ; and hence, …1 ÿ ckn †En ÿ Enÿ1 6 ckn gn : We introduce a sequence of numbers zn ˆ

…4:24† fzn gNnˆ0

by z0 ˆ 1 and zn ˆ …1 ÿ ckn †znÿ1 ; 1 6 n 6 N . Obviously,

n Y …1 ÿ ckj †: jˆ1

Using the inequalities eÿ2ckj 6 1 ÿ ckj 6 eÿckj

if kj 6 k is sufficiently small

we have the following bounds on zn : eÿ2ctn 6 zn 6 eÿctn

if k is sufficiently small:

…4:25†

190

W. Han, M. Sofonea / Comput. Methods Appl. Mech. Engrg. 190 (2000) 179±191

Now we rewrite (4.24), zn En ÿ znÿ1 Enÿ1 6 ckn znÿ1 gn : A simple induction argument yields z n En 6 c

n X

kj zjÿ1 gj :

jˆ1

Using the bounds (4.25), we obtain En 6 c

n X

k j gj :

jˆ1

Then from (4.23), e n 6 c gn ‡

n X

! k j gj ;

n ˆ 1; . . . ; N :

jˆ1

Therefore,

ÿ

max

16n6N


hk krn ÿ rhk n kH ‡ kun ÿ un kV 6 c max gn :

…4:26†

16n6N

Notice that the error bound (4.26) for the fully discrete solution contains an additional term In , resulted from numerical integration of the constitutive relation (2.14). This term can be estimated as follows. n Z tj X jG…r…t†; e…u…t††† ÿ G…rj ; e…uj ††j dt In 6 tjÿ1

jˆ1

6c

n X jˆ1



Z

tj

tjÿ1

ÿ


kr…t† ÿ rj kH ‡ ke…u…t† ÿ uj †kH dt

 _ L1 …0;T ;H† : _ L1 …0;T ;H† ‡ ke…u†k 6 ctn k krk Therefore, we obtain the error estimate (4.20) from the estimate (4.26).



Theorem 4.2 is the basis for convergence analysis and error estimation. The estimate (3.15) is now modi®ed to jRn …un ÿ vhn †j 6 c…kf n kV 0 ‡ krn kH †kun ÿ vhn kV :

…4:27†

Then again arguing as in [8], we have the following convergence result. h hk h Corollary 4.3. Assume the initial values rhk 0 2 Q and u0 2 U are chosen in such a way that

kr0 ÿ rhk 0 kH ! 0;

ku0 ÿ uhk 0 kV ! 0

as h; k ! 0:

Also assume the set …H 2 …X††d \ U is dense in U and the condition (3.16) is satisfied. Then the fully discrete method converges: ÿ
hk as h; k ! 0: max krn ÿ rhk n kH ‡ kun ÿ un kV ! 0 16n6N

Again we comment that assumption (3.16) is not a restriction in applications. Under the regularity assumption (3.17), we can improve the estimate (4.27). Indeed, as in (3.18), jRn …un ÿ vh †j 6 kr1m …tn †kL2 …C3 † kvhm ÿ um …tn †kL2 …C3 † : Thus we have proved the following error estimate.

W. Han, M. Sofonea / Comput. Methods Appl. Mech. Engrg. 190 (2000) 179±191

Corollary 4.4. Assume additionally the regularity condition (3.17). Then ÿ
hk max krn ÿ rhk n kH ‡ kun ÿ un kV 16n6N   hk _ _ k ‡ ku ÿ u k † ‡ ck k rk ‡ k uk 6 c…kr0 ÿ rhk 1 1 0 L …0;T ;H† L …0;T ;V † 0 H 0 V n o 1=2 1=2 inf ku…tn † ÿ vh kV ‡ kr1m kL1 …0;T ;L2 …C3 †† kum …tn † ÿ vhm kL2 …C3 † : ‡ max 1 6 n 6 N vh 2L1 …0;T ;U h †

191

…4:28†

Inequalities (4.20) and (4.28) are the basis for deriving order error estimates for the fully discrete solution. We present one sample result. Let us use linear elements for the ®nite element space V h , and piecewise h hk h constants for Qh . We choose the initial values uhk 0 2 U and r0 2 Q in such a way that kr0 ÿ rhk 0 kH 6 ch; ku0 ÿ uhk 0 kV 6 ch: Then as in the case of the spatially semi-discrete approximations discussed in the last section, under the regularity assumptions d

u 2 L1 …0; T ; …H 2 …X†† †;

um 2 L1 …0; T ; H 2 …C3 †† and

r1m 2 L1 …0; T ; L2 …C3 ††;

we have the optimal order error estimate ÿ
hk max krn ÿ rhk n kH ‡ kun ÿ un kV 6 c…k ‡ h†: 16n6N

Acknowledgements The work of W. Han was supported by the Carver Scienti®c Initiative Grant, 1998. References [1] A. Amassad, M. Sofonea, Analysis of a quasistatic viscoplastic problem involving Tresca friction law, Discrete and Continuous Dynamical Systems 4 (1997) 55±72. [2] M. Burguera, J.M. Via~ no, Numerical method for stress approximation in unilateral contact between two perfect plastic bodies, in: D.R.J. Owen, E. Hinton, E. O~ nate (Eds.), Computational Plasticity II, Models Software and Applications, Part II, Pineridge Press, Swansea, 1989, pp. 1099±1110. [3] M. Burguera, J.M. Via~ no, Numerical solving of frictionless contact problems in perfectly plastic bodies, Comput. Meth. Appl. Mech. Engrg. 121 (1995) 303±322. [4] J. Chen, W. Han, M. Sofonea, Numerical analysis of a contact problem in rate-type viscoplasticity, submitted. [5] P.G. Ciarlet, The Finite Element Method for Elliptic Problems, North-Holland, Amsterdam, 1978. [6] N. Cristescu, I. Suliciu, Viscoplasticity, Martinus Nijho€ Publishers, Editura Tehnica, Bucarest, 1982. [7] W. Han, B.D. Reddy, Computational plasticity: the variational basis and numerical analysis, Comput. Mech. Adv. 2 (1995) 283±400. [8] W. Han, B.D. Reddy, Plasticity: Mathematical Theory and Numerical Analysis, Springer-Verlag, New York, 1999. [9] J. Haslinger, I. Hlavacek, Contact between elastic perfectly plastic bodies, Appl. Mat. 27 (1982) 27±45. [10] J. Haslinger, I. Hlavacek, J. Necas, Numerical methods for unilateral problems in solid mechanics, in: P.G. Ciarlet, J.L. Lions (Eds.), Handbook of Numerical Analysis, Elsevier, Amsterdam, 1996, pp. 313±485. [11] I. Hlavacek, J. Haslinger, J. Nec as, J. Lovisek, Solution of Variational Inequalities in Mechanics, Springer, New York, 1988. [12] P. Hild, Problemes de contact unilateral et maillages elements ®nis incompatibles, These, Universite Paul Sabatier de Toulouse, 1998. [13] I. Ionescu, M. Sofonea, Functional and Numerical Methods in Viscoplasticity, Oxford University Press, Oxford, 1993. [14] N. Kikuchi, J.T. Oden, Contact Problems in Elasticity, SIAM, Philadelphia, 1988. [15] J. Necas, I. Hlavacek, Mathematical Theory of Elastic and Elastoplastic Bodies: An Introduction, Elsevier, Amsterdam, 1981. [16] M. Rochdi, M. Sofonea, On frictionless contact between two elastic±viscoplastic bodies, Q. J. Mech. Appl. Math. 50 (1997) 481±496. [17] M. Sofonea, On a contact problem for elastic±viscoplastic bodies, Nonlinear Anal. TMA 29 (1997) 1037±1050. [18] E. Zeidler, Nonlinear Functional Analysis and its Applications, Volume II/A: Linear Monotone Operators, Springer, New York, 1990.