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

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.

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)

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., ). 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 vH ∀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 , 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 + hw − Π 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 qQ  qQ

∀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  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 . 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 . 1nN

j =1

1nN

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) + ˙vL1 (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 + cvn − 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 + cvn − 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  + hvj − 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  chvj − vj +1 V  ch˙vL1 (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

1nN

 max

1nN

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). 1nN

509

n=1

1nN

(3.32)