**VARIATIONAL AND NUMERICAL ANALYSIS** **OF THE SIGNORINI’S CONTACT PROBLEM** **IN VISCOPLASTICITY WITH DAMAGE**

**J. R. FERNÁNDEZ AND M. SOFONEA**

*Received 11 February 2002*

We consider the quasistatic Signorini’s contact problem with damage for elastic-viscoplastic bodies. The mechanical damage of the material, caused by excessive stress or strain, is described by a damage function whose evolution is modeled by an inclusion of parabolic type. We pro- vide a variational formulation for the mechanical problem and sketch a proof of the existence of a unique weak solution of the model. We then in- troduce and study a fully discrete scheme for the numerical solutions of the problem. An optimal order error estimate is derived for the approx- imate solutions under suitable solution regularity. Numerical examples are presented to show the performance of the method.

**1. Introduction**

We consider a mathematical model for a quasistatic process of friction- less contact between an elastic-viscoplastic body and an obstacle within the framework of small deformation theory. The contact is modeled with the classical Signorini’s conditions in a form with a gap function. The eﬀect of damage due to mechanical stress or strain is included in the model. Such situations are common in many engineering applications where the forces acting on the system vary periodically, leading to the appearance and growth of microcracks which may deteriorate the mech- anism of the system. Because of the importance of the safety issue of me- chanical equipments, considerable eﬀort has been devoted to modeling and numerically simulating damage.

Early models for mechanical damage derived from thermomechan- ical considerations appeared in[15, 16], where numerical simulations

Copyright^{}^{c}2003 Hindawi Publishing Corporation
Journal of Applied Mathematics 2003:2(2003)87–114

2000 Mathematics Subject Classification: 74M15, 74S05, 74C10, 74R20 URL:http://dx.doi.org/10.1155/S1110757X03202023

were included. One-dimensional damage problems have been studied in[13,14]. Recently, the existence of weak solutions of viscoelastic prob- lems with friction and damage have been provided in[20,27]. A qua- sistatic frictionless contact problem for elastic-viscoplastic materials with normal compliance and damage was studied in[4].

In the present paper, we consider a rate-type elastic-viscoplastic ma- terial with constitutive relation

* σ*˙ =Eε(

**u) +**˙

*G*

**σ,ε(u), β**

*,* (1.1)

where * σ* represents the stress tensor field,

**u**denotes the displacement field, and

*is the linearized strain tensor field. Here,Eis a fourth- order tensor,*

**ε(u)***G*is a nonlinear constitutive function, and

*β*is the damage field. The latter is related to the inelastic part of the stress and its values are restricted to the interval[0,1]. When

*β*=1, the material is undam- aged, while the value

*β*=0 indicates the stage of complete damage, and for 0

*< β <*1, there is partial damage. In(1.1)and everywhere in what fol- lows, the dot above a variable represents the time derivative. Note that for particular forms of the function

*G, the constitutive law*(1.1)may de- scribe a viscoelastic behavior.

Rate-type viscoplastic constitutive laws of the form(1.1)in which the
function*G*does not depend on*β*were considered by many authors, see
for instance,[8,23]and the references therein. Frictionless contact prob-
lems for such kind of materials were studied in[10,12,22,28]. A fric-
tionless contact problem for materials of the form(1.1)in which*β*is an
internal state variable whose evolution is described by an ordinary dif-
ferential equation has been recently considered in[11].

One of the traits of novelty of this paper consists in the fact that, fol- lowing[15,16], the evolution of the microscopic cracks responsible for the damage is modeled by the following diﬀerential inclusion:

*β*˙−*κβ*+*∂ψ*_{K}(β)*φ*

**σ,ε(****u**), β

*.* (1.2)

In(1.2) and below,*κ >*0 is a constant, *∂ψ*_{K} denotes the subdiﬀerential
of the indicator function*ψ*KofK, which represents the set of admissible
damage functions satisfying 0≤*β*≤1, and*φ*is a given constitutive func-
tion which describes damage sources in the system. In[13,14,15,16],
the damage was assumed irreversible and therefore the condition ˙*β*≤0
was imposed. In this paper, we assume that the material may recover
from damage and cracks may close, and thus, we do not impose this
restriction. In[15,16], the damage-source function was chosen to be un-
bounded when *β*→0, a condition which is not allowed under our as-
sumptions on*φ. Therefore, we may consider the global solution, which*
we establish below for our problem as local solution of a problem with

the damage source used in[15,16], valid as long as an inequality of the
form*β*≥*β*∗*>*0 holds.

Our purpose of this paper is threefolds. First, we provide a variational analysis of the mechanical problem and briefly show the existence of a unique weak solution for the model. We then introduce a fully discrete scheme and derive error estimates. Finally, we report some numerical re- sults on the performance of the numerical method considered. Literature on the study of variational inequalities is rather extensive, see, for exam- ple, the monographs[1,18,24]. In particular, some results on numerical analysis of variational inequalities we use here can be found in[19,21].

The rest of the paper is organized as follows. In Section2, we present the mechanical problem and provide its variational analysis including an existence and uniqueness result, Theorem 2.3. The proof of Theo- rem2.3is based on classical results of elliptic and parabolic variational inequalities and Banach fixed point theorem. In Section3, we analyze a fully discrete scheme for the problem. We use the finite-element method to discretize the spatial domain and a backward Euler finite diﬀerence to discretize the time derivative. We obtain an optimal order error es- timate under appropriate regularity assumptions on the exact solution and data. Finally, in Section4, we give some numerical examples to show the performance of the scheme.

**2. Mechanical problem and variational formulation**

The physical setting is as follows. A viscoplastic body occupies the do-
mainΩ⊂R* ^{d}* (d=1,2,3 in applications) with outer Lipschitz surfaceΓ
that is divided into three disjoint measurable partsΓ1,Γ2, andΓ3 such
that meas(Γ1)

*>*0. Let [0, T] be the time interval of interest. The body is clamped on Γ1, a volume force of density

**f**0 acts in Ω and surface tractions of density

**f**

_{2}act onΓ2. The functions

**f**

_{0}and

**f**

_{2}can depend on the time variable. The body may come in frictionless contact onΓ3with an obstacle, the so-called foundation. We assume that the foundation is rigid and therefore we model the contact with the classical Signorini’s conditions in a form with a gap function. Finally, we use(1.1)and(1.2) to describe the viscoplastic behaviour of the material. Then, the classical form of the mechanical problem is as follows.

*Problem 2.1.* Find a displacement field**u**:Ω×[0, T]→R* ^{d}*, a stress field

*:Ω×[0, T]→*

**σ***S*

*d*, and a damage field

*β*:Ω×[0, T]→Rsuch that

* σ*˙ =Eε(

**u) +**˙

*G*

**σ,ε(u), β**

inΩ×(0, T), (2.1)

*β*˙−*κβ*+*∂ψ*K(β)*φ*

**σ,ε(u), β**

inΩ×(0, T), (2.2)

Divσ+**f**0=**0** inΩ×(0, T), (2.3)

**u**=**0** onΓ1×(0, T), (2.4)

* σν*=

**f**2 onΓ2×(0, T), (2.5)

*u**ν*≤*g,* *σ**ν*≤0, *σ**ν*

*u**ν*−*g*

=0, **σ***τ* =**0** onΓ3×(0, T), (2.6)

*∂β*

*∂ν* =0 onΓ×(0, T), (2.7)

**u**(0) =**u**_{0}*,* * σ(0) =σ*0

*,*

*β(0) =β*0 inΩ. (2.8)

Here,*S**d*denotes the space of second-order symmetric tensors onR* ^{d}*,
(2.3) represents the equilibrium equation in which Div

*denotes the divergence of the stress, while(2.4)and(2.5)are the displacement and traction boundary conditions, respectively. In contact conditions (2.6), we used*

**σ***u*

*ν*,

*σ*

*ν*, and

**σ***τ*to denote the normal displacement, the normal stress, and the tangential stress, respectively, and

*g*represents the initial gap between the potential contact surfaceΓ3 and the foundation, mea- sured along the outward normal vector

*to Γ. The fourth relation in (2.6)indicates that the friction force on the contact surface vanishes, that is, the contact is frictionless. Equation(2.7)describes the homogeneous Neumann boundary condition for the damage field which we use here for simplicity, according to[20,27]. Finally, in(2.8),*

**ν****u**

_{0},

*0, and*

**σ***β*0are the initial data for the displacement, stress, and damage field, respectively.

We introduce the notation to be used in the rest of the paper. Further
details can be found in[23,24,26]. In the sequel, “·” and| · |represent the
inner product and the Euclidean norm on both*S**d* andR* ^{d}*, respectively,
and we use the following spaces:

*H*=

*L*^{2}(Ω)*d*

*,* *Q*=

* σ*=

**σ***ij*

|**σ***ij*=**σ***ji*∈*L*^{2}(Ω)
*,*
*H*1=

*H*^{1}(Ω)*d*

*,* *Q*1=

* σ*∈

*Q*|

**σ***ij,j*∈

*H*

*.* (2.9)

Here and below,*i, j* =1, . . . , d, summation over repeated indices is im-
plied, and the index that follows a comma indicates a partial derivative
while*H,Q,H*1, and*Q*1 are real Hilbert spaces endowed with the inner
products given by

(u,**v)***H*=

Ω*u**i**v**i**dx,* (σ,**τ)***Q*=

Ω*σ**ij**τ**ij**dx,*
(u,**v)***H*1= (u,**v)***H*+

**ε(u),ε(v)**

*Q**,*
(σ,τ)*Q*1= (σ,τ)*Q*+ (Divσ,Divτ)*H**,*

(2.10)

respectively. Here* ε*:

*H*1→

*Q*and Div :

*Q*1→

*H*are the

*deformation*and

*divergence*operators

**ε(v) =****ε***ij*(v)

*,* *ε**ij*(v) =1
2

*v**i,j*+*v**j,i*

*,* Div* σ*=

**σ***ij,j*

*.* (2.11)
The associated norms on the spaces*H,Q,H*1, and*Q*1 are denoted by

| · |*H*,| · |*Q*,| · |*H*1, and| · |*Q*1, respectively.

Since the boundaryΓis Lipschitz continuous, the unit outward nor-
mal vector* ν* is defined a.e. on Γ. For

**v**∈

*H*1, we again write

**v**for the trace

*γv*of

**v**onΓ, and we denote by

*v*

*ν*and

**v**

*the*

_{τ}*normal*and

*tangential*components of

**v**on the boundary given by

*v*

*ν*=

**v**·

*and*

**ν****v**

*τ*=

**v**−

*v*

*ν*

**ν.**For a regular(say*C*^{1})tensor field* σ*:Ω→

*S*

*d*, we define its normal and tangential components by

*σ*

*ν*= (σν)·

*and*

**ν**

**σ***τ*=

*−*

**σν***σ*

*ν*

*and we recall that the following Green’s formula holds:*

**ν****σ,ε(****v**)

*Q*+ (Divσ,v)*H*=

Γ* σν*·

**vda,**∀

**v**∈

*H*1

*.*(2.12) Let

*V*denote the closed subspace of

*H*1defined by

*V*=

**v**∈*H*1|**v**=**0**onΓ1

*.* (2.13)

Since meas(Γ1)*>*0, Korn’s inequality holds that there exists*C**K**>*0 de-
pending only onΩandΓ1such that

**ε(v)*** _{Q}*≥

*C*

*K*|v|

*H*1

*,*∀v∈

*V.*(2.14) A proof of Korn’s inequality may be found in[25, page 79]. On

*V*, we consider the inner product given by

(**u,v**)*V* =

**ε(****u**),**ε(****v**)

*Q**,* ∀**u,v**∈*V,* (2.15)
and let | · |*V* be the associated norm, that is,|v|*V* =|ε(v)|*Q* for **v**∈*V*. It
follows that| · |*H*1and| · |*V* are equivalent norms on*V*and therefore(V,| ·

|*V*)is a real Hilbert space.

Denote by*U*the convex subset of admissible displacements defined
by

*U*=

**v**∈*V* |*v**ν*≤*g*onΓ3

*.* (2.16)

If*X*1and*X*2are real Hilbert spaces, then*X*1×*X*2denotes the product
space which is endowed with the canonical inner product(·,·)*X*1×X2. Fi-
nally, if*X*is a real Hilbert space, we denote by| · |*X* the norm on*X. For*

*T >*0, we use the standard notation for *L** ^{p}*(0, T;

*X)*and Sobolev spaces

*W*

*(0, T;X),*

^{k,p}*k*∈N, 1≤

*p*≤ ∞.

In the study of the mechanical Problem2.1, we assume that the elas-
ticity tensorE= (E*ijkh*):Ω×*S**d*→*S**d*satisfies

E*ijkh*∈*L*^{∞}(Ω),

Eσ·* τ*=

*· Eτ*

**σ***,*∀σ,τ∈

*S*

*d*

*,*a.e. inΩ, Eσ·

*≥*

**σ***α|σ|*

^{2}

*,*∀σ∈

*S*

*d*

*,*for some

*α >*0.

(2.17)

The viscoplastic function*G*:Ω×*S**d*×*S**d*×R→*S**d*satisfies
(a)There exists*L >*0 such that

*G*

**x,σ**1*, ε*1

*, β*1

−*G*

**x,*** σ*2

*,*2

**ε***, β*2≤

*L*1−

**σ***2+*

**σ***1−*

**ε***2+*

**ε***β*1−

*β*2

*,*

∀σ1*, σ*2

*,ε*1

*,*2∈

**ε***S*

*d*

*,*

*β*1

*, β*2∈R

*,*a.e.

**x**∈Ω,

(b)**x**−→*G(***x,σ,*** ε, β)*is a Lebesgue measurable function onΩ,

∀σ,* ε*∈

*S*

*d*

*,*

*β*∈R

*,*(c)

**x**−→

*G(x,*

**0,0,**0)∈

*Q.*

(2.18)

The damage source function*φ*:Ω×*S**d*×*S**d*×R→Rsatisfies
(a)There exists ˜*L >*0 such that

*φ*

*x,σ*1*, ε*1

*, β*1

−*φ*

*x,σ*2*, ε*2

*, β*2≤

*L*˜

*1−*

**σ***2+*

**σ***1−*

**ε***2+*

**ε***β*1−

*β*2

*,*

∀σ1*,σ*2*, ε*1

*,ε*2∈

*S*

*d*

*,*

*β*1

*, β*2∈R

*,*a.e.

**x**∈Ω,

(b)**x**−→*φ(***x,*** σ,ε, β)*is a Lebesgue measurable function onΩ,

∀σ,ε∈*S**d**,* *β*∈R*,*

(c)**x**−→*φ(***x,0,0,**0)∈*L*^{2}(Ω).

(2.19)

The body forces and surface tractions have the regularity
**f**_{0}∈*W*^{1,2}(0, T;*H),* **f**_{2}∈*W*^{1,2} 0, T;

*L*^{2}
Γ2

*d*

*.* (2.20)

The gap function*g*is such that
*g*∈*L*^{2}

Γ3

*,* *g*(x)≥0, a.e.**x**∈Γ3*,* (2.21)

and the initial data satisfy

**u**0∈*V,* * σ*0∈

*Q*1

*,*(2.22)

* σ*0

*,ε*

**v**−**u**0

*Q*≥

**f(0),v**−**u**0

*V**,* ∀v∈*U,* (2.23)
*β*0∈*H*^{1}(Ω), 0*< β*_{∗}≤*β*0≤1, a.e. inΩ. (2.24)
Here,**f**:[0, T]→*V* is the function defined by

**f**(t),v

*V* =

**f**_{0}(t),v

*H*+

**f**_{2}(t),v

[L^{2}(Γ2)]^{d}*,* ∀**v**∈*V,*∀t∈[0, T]. (2.25)

Notice that conditions(2.20)imply

**f**∈*W*^{1,2}(0, T;*V*). (2.26)

Let*a*:*H*^{1}(Ω)×*H*^{1}(Ω)→Rbe the bilinear form
*a(ξ, ψ) =κ*

Ω∇ξ· ∇ψ dx, ∀ξ, ψ∈*H*^{1}(Ω), (2.27)
and letKdenote the set of admissible damage functions

K=

*ξ*∈*H*^{1}(Ω)|0≤*ξ*≤1 inΩ

*.* (2.28)

By a standard procedure, we can derive the following variational for- mulation of the mechanical Problem2.1.

*Problem 2.2.* Find a displacement field **u**:[0, T]→*V*, a stress field * σ*:
[0, T]→

*Q*1, and a damage field

*β*:[0, T]→

*H*

^{1}(Ω)such that

* σ(t) =*˙ Eε

˙
**u**(t)

+*G*
**σ(t),****ε**

**u**(t)
*, β(t)*

*,* a.e.*t*∈(0, T), (2.29)
**u**(t)∈*U,*

**σ(t),****ε**

**v**−**u**(t)

*Q*≥

**f**(t),**v**−**u**(t)

*V**,*

∀**v**∈*U,* ∀t∈[0, T], (2.30)

*β(t)*∈ K, a.e.*t*∈(0, T), (2.31)

*β(t), ξ*˙ −*β(t)*

*L*^{2}(Ω)+*a*

*β(t), ξ*−*β(t)*

≥
*φ*

**σ(t),ε****u(t)**

*, β(t)*

*, ξ*−*β(t)*

*L*^{2}(Ω)*,*

∀ξ∈ K, a.e.*t*∈(0, T),

(2.32)
**u(0) =u**0*,* * σ(0) =σ*0

*,*

*β(0) =β*0 inΩ. (2.33)

Concerning the well-posedness of Problem2.2, we have the following result.

Theorem2.3. *Assume (2.17), (2.18), (2.19), (2.20), (2.21), (2.22), (2.23), and*
*(2.24). Then there exists a unique solution* (**u,****σ, β)***of Problem* *2.2with the*
*regularity*

**u**∈*W*^{1,2}(0, T;V), * σ*∈

*W*

^{1,2}0, T;Q1

*,*
*β*∈*W*^{1,2}

0, T;*L*^{2}(Ω)

∩*L*^{2}

0, T;*H*^{1}(Ω)

*.* (2.34)

*Proof.* The proof or Theorem2.3is based on fixed-point type arguments
similar to those used in [4] but with a diﬀerent choice of the opera-
tors. Since the modifications are straightforward, we omit the details.

The main steps of the proof are stated as follows.

(i)For any* η*= (η

^{1}

*, η*

^{2})∈

*L*

^{2}(0, T;

*Q*×

*L*

^{2}(Ω)), let

**z**^{1}* _{η}*(t) =

*t*

0

**η**^{1}(s)ds+* σ*0− Eε

**u**

_{0}

*.* (2.35)

Then,**z**^{1}* _{η}*∈

*W*

^{1,2}(0, T;Q)and there exists a unique solution(

**u**

_{η}*,*

**σ***η*)of the problem

**σ***η*(t) =Eε
**u***η*(t)

+**z**^{1}* _{η}*(t), ∀t∈[0, T], (2.36)

**u**(t)∈

*U,*

**σ***η*(t),**ε**

**v**−**u**(t)

*Q*≥

**f**(t),**v**−**u**(t)

*V**,*

∀**v**∈*U,* ∀t∈[0, T], (2.37)

**u*** _{η}*(0) =

**u**

_{0}

*,*

**σ***η*(0) =

*0*

**σ***.*(2.38) Moreover, the solution satisfies

**u**

*η*∈

*W*

^{1,2}(0, T;V)and

**σ***η*∈

*W*

^{1,2}(0, T;

*Q*1).

(ii)For any* η*= (η

^{1}

*, η*

^{2})∈

*L*

^{2}(0, T;

*Q*×

*L*

^{2}(Ω)), there exists a unique so- lution

*β*

*η*of the problem

*β**η*(t)∈ K, a.e.*t*∈(0, T), (2.39)
*β*˙*η*(t), ξ−*β(t)*

*L*^{2}(Ω)+*a*

*β**η*(t), ξ−*β**η*(t)

≥

*η*^{2}(t), ξ−*β**η*(t)

*L*^{2}(Ω)*,*

∀ξ∈ K, a.e.*t*∈(0, T), (2.40)

*β**η*(0) =*β*0*.* (2.41)

Moreover, the solution satisfies*β**η*∈*W*^{1,2}(0, T;*L*^{2}(Ω))∩*L*^{2}(0, T;*H*^{1}(Ω)).

(iii)Consider the Banach space*X*=*L*^{2}(0, T;*Q*×*L*^{2}(Ω))with the norm

| · |*X*given by

|η|^{2}* _{X}*=

_{T}0

**η**^{1}(s)^{2}

*Q*+*η*^{2}(s)^{2}

*L*^{2}(Ω)

*ds,* ∀η=
**η**^{1}*, η*^{2}

∈*X,* (2.42)

and define the operatorΛ:*X*→*X*by
Λη(t) =

*G*

**σ***η*(t),ε
**u*** _{η}*(t)

*, β**η*(t)
*, φ*

**σ***η*(t),ε
**u*** _{η}*(t)

*, β**η*(t)

*,* (2.43)
for any* η*∈

*X,t*∈[0, T]. Then, the operatorΛhas a unique fixed point

**η**^{∗}∈

*X.*

(iv)Let**η**^{∗}= (η^{∗1}*, η*^{∗2})∈*X*be the fixed point ofΛand denote**u**=**u***η*^{∗},
* σ*=

**σ***η*

^{∗}, and

*β*=

*β*

*η*

^{∗}, where(

**u**

_{η}^{∗}

*,*

**σ***η*

^{∗})is the solution of problem(2.36), (2.37), and(2.38) for

*=*

**η**

**η**^{∗}and

*β*

*η*

^{∗}is the solution of problem (2.39), (2.40), and(2.41)for

*=*

**η**

**η**^{∗}. Then,(u,

*is the unique solution of Prob-*

**σ, β)**lem2.2which satisfies(2.34).

We conclude by Theorem2.3 that, under assumptions(2.17),(2.18),
(2.19),(2.20),(2.21), (2.22),(2.23), and(2.24), the mechanical problem
(2.1), (2.2), (2.3), (2.4),(2.5),(2.6),(2.7), and (2.8)has a unique weak
solution(**u,σ, β), with regularity**(2.34).

**3. Numerical approximation**

We analyze in this section a fully discrete approximation scheme for
Problem2.2. To this end, we suppose in the sequel that conditions(2.17),
(2.18),(2.19),(2.20),(2.21),(2.22),(2.23), and(2.24)hold. We consider
arbitrary finite-dimensional spaces*V** ^{h}*⊂

*V*and

*Q*

*⊂*

^{h}*Q, and let*K

*⊂ K be a nonempty, finite-dimensional closed convex set. Here*

^{h}*h >*0 is a dis- cretization parameter and we assume that

**ε(V***)⊂*

^{h}*Q*

*. This assumption is not a restriction for actual implementation of the method since, usu- ally,*

^{h}*V*

*and*

^{h}*Q*

*are constructed to be finite-element spaces and*

^{h}*V*

*con- sists of continuous piecewise polynomials of a degree one higher than that of*

^{h}*Q*

*. Finally, denote by*

^{h}*U*

*⊂*

^{h}*V*

*an approximation for the convex set*

^{h}*U*for which we assume the following conformity condition:

*U** ^{h}*⊂

*U.*(3.1)

LetP*Q** ^{h}*:

*Q→Q*

*be the orthogonal projection operator defined through the relation*

^{h}P*Q*^{h}**q,q**^{h}

*Q*=
**q,q**^{h}

*Q**,* ∀**q**∈*Q,* **q*** ^{h}*∈

*Q*

^{h}*.*(3.2)

The orthogonal projection operator is nonexpansive, that is,

P*Q*^{h}**q*** _{Q}*≤ |q|

*Q*

*,*∀q∈

*Q.*(3.3) This property will be used on several occasions.

We use a uniform partition of time interval[0, T]with the step-size*k*=
*T/N*and the nodes*t**n*=*nk*for*n*=0,1, . . . , N. The extension of the discus-
sion here to the case of nonuniform partition does not present any diﬃ-
culty. For a continuous function*z(t), we use the notationz**n*=*z(t**n*). For
a sequence{z*n*}^{N}* _{n=0}*, we denote

*δz*

*n*= (z

*n*−

*z*

*n−1*)/kfor the corresponding divided diﬀerence. No summation is implied over the repeated index

*n.*

In the rest of this section,*c*will denote positive constants which are in-
dependent of the discretization parameters*h*and*k.*

Let**u**^{h}_{0} ∈*U** ^{h}*,

**σ**

^{h}_{0}∈

*Q*

*, and*

^{h}*β*

_{0}

*∈ K*

^{h}*be chosen to approximate the initial values*

^{h}**u**

_{0},

*0, and*

**σ***β*0. A fully discrete approximation scheme for Prob- lem2.2is the following problem.

*Problem 3.1.* Find **u*** ^{hk}* ={

**u**

^{hk}*}*

_{n}

^{N}*⊂*

_{n=0}*U*

*,*

^{h}

**σ***={σ*

^{hk}

^{hk}*n*}

^{N}*⊂*

_{n=0}*Q*

*, and*

^{h}*β*

*= {β*

^{hk}

^{hk}*n*}

^{N}*⊂ K*

_{n=0}*such that*

^{h}**u**^{hk}_{0} =**u**^{h}_{0}*,* **σ**^{hk}_{0} =**σ**^{h}_{0}*,* *β*^{hk}_{0} =*β*^{h}_{0}*,* (3.4)
and for*n*=1,2, . . . , N,

*δσ*^{hk}*n* =P*Q** ^{h}*Eε

*δu*

^{hk}

_{n}+P*Q*^{h}*G*

**σ**^{hk}_{n−1}*, ε*

**u**

^{hk}

_{n−1}*, β*^{hk}_{n−1}*,*

**σ**^{hk}_{n}*,ε*

**v*** ^{h}*−

**u**

^{hk}

_{n}*Q*≥

**f***n**,***v*** ^{h}*−

**u**

^{hk}

_{n}*V**,* ∀v* ^{h}*∈

*U*

^{h}*,*

*δβ*

^{hk}

_{n}*, ξ*

*−*

^{h}*β*

^{hk}

_{n}*L*^{2}(Ω)+*a*

*β*^{hk}_{n}*, ξ** ^{h}*−

*β*

_{n}

^{hk}≥
*φ*

**σ**^{hk}_{n−1}*, ε*

**u**

^{hk}

_{n−1}*, β*_{n−1}^{hk}

*, ξ** ^{h}*−

*β*

^{hk}

_{n}*L*^{2}(Ω)*,* ∀ξ* ^{h}*∈ K

^{h}*.*(3.5) By induction, we obtain that this problem is equivalent to

**σ**^{hk}* _{n}* =

**σ**

^{hk}_{0}− P

*Q*

*Eε*

^{h}**u**

^{h}_{0}

+P*Q** ^{h}*Eε

**u**

^{hk}*+k*

_{n}

^{n}*j=1*

P*Q*^{h}*G*
**σ**^{hk}_{j−1}*,ε*

**u**^{hk}_{j−1}*, β*^{hk}_{j−1}

*,* (3.6)

**σ**^{hk}*n* *,ε*

**v*** ^{h}*−

**u**

^{hk}

_{n}*Q*≥

**f**_{n}*,***v*** ^{h}*−

**u**

^{hk}

_{n}*V**,* ∀**v*** ^{h}*∈

*U*

^{h}*,*(3.7)

*δβ*

^{hk}

_{n}*, ξ*

*−*

^{h}*β*

^{hk}

_{n}*L*^{2}(Ω)+*a*

*β*^{hk}_{n}*, ξ** ^{h}*−

*β*

_{n}

^{hk}≥
*φ*

**σ**^{hk}_{n−1}*, ε*

**u**

^{hk}

_{n−1}*, β*_{n−1}^{hk}

*, ξ** ^{h}*−

*β*

^{hk}

_{n}*L*^{2}(Ω)*,* ∀ξ* ^{h}*∈ K

^{h}*.*(3.8)

For given **σ**^{hk}* _{j−1}*,

**ε(****u**

^{hk}*),*

_{j−1}*β*

_{j−1}*, 1≤*

^{hk}*j*≤

*n, we can first determineβ*

_{n}*from (3.8)which has a unique solution by a classical result on elliptic varia- tional inequalities. Combining(3.6)and(3.7), we have*

^{hk}Eε
**u**^{hk}_{n}

*, ε*

**v*** ^{h}*−

**u**

^{hk}

_{n}*Q*≥

**f***n**,v** ^{h}*−

**u**

^{hk}

_{n}*V*−

**σ**^{h}_{0}− Eε
**u**^{h}_{0}

*,ε*

**v*** ^{h}*−

**u**

^{hk}

_{n}*Q*

−*k*
*n*
*j=1*

*G*
**σ**^{hk}_{j−1}*,ε*

**u**^{hk}_{j−1}*, β*^{hk}_{j−1}

*,ε*

**v*** ^{h}*−

**u**

^{hk}

_{n}*Q**,*

∀**v*** ^{h}*∈

*U*

^{h}*.*

(3.9)
By using classical results on variational inequalities (see, for instance,
[17, Chapter IV]), we see that(3.9)has a unique solution**u**^{hk}* _{n}* ∈

*U*

*. To solve variational inequality(3.9)for*

^{h}**u**

^{hk}*and variational inequality(3.8) for*

_{n}*β*

^{hk}*, a penalty-duality algorithm can be used(see[2,3]). Once*

_{n}**u**

^{hk}*is known, we can determine*

_{n}

**σ**

^{hk}*n*from(3.6). So, an induction argument shows that the fully discrete scheme has a unique solution. In the same way, we obtain that variational inequality (3.8) has a unique solution

*β*

^{hk}*∈ K*

_{n}*. We summarize this result as follows.*

^{h}Theorem3.2. *Assume (2.17), (2.18), (2.19), (2.20), (2.21), (2.22), (2.23), and*
*(2.24). Then there exists a unique solution (u*^{hk}*, σ*

^{hk}*, β*

^{hk}*) of Problem3.1.*

Now, we proceed to derive error estimates for the discrete solution.

Integrating(2.29)at time*t*=*t**n*, we obtain the following relations for the
solution of Problem2.2(n=1, . . . , N):

**σ***n*=* σ*0− Eε

**u**

_{0}

+Eε
**u**_{n}

+
_{t}_{n}

0

*G*

**σ(s),ε****u**(s)

*, β(s)*

*ds,* (3.10)
**σ***n**,ε*

**v**−**u**_{n}

*Q*≥

**f**_{n}*,v*−**u**_{n}

*V**,* ∀**v**∈*U,* (3.11)
*β*˙*n**, ξ*−*β**n*

*L*^{2}(Ω)+a

*β**n**, ξ*−*β**n*

≥
*φ*

**σ***n**,ε*
**u***n*

*, β**n*
*, ξ*−*β**n*

*L*^{2}(Ω)*,* ∀ξ∈K.

(3.12) Subtracting(3.10)and(3.6), we find

**σ***n*−**σ**^{hk}* _{n}* =

*0−*

**σ**

**σ**

^{h}_{0}−

*I*− P

*Q*

^{h}Eε
**u***n*−**u**0

− P*Q** ^{h}*Eε

**u**0−**u**^{h}_{0}
+D*G,n*+P*Q** ^{h}*Eε

**u*** _{n}*−

**u**

^{hk}*+k*

_{n}

^{n}*j=1*

*I*−P*Q*^{h}

*G*
**σ***j−1**,ε*

**u**_{j−1}*, β**j−1*
+*kP**Q*^{h}

*n*
*j=1*

*G*
**σ***j−1**, ε*

**u**_{j−1}*, β**j−1*

−*G*
**σ**^{hk}_{j−1}*, ε*

**u**^{hk}_{j−1}*, β*^{hk}_{j−1}

*,*
(3.13)

where
*D**G,n*=

*t**n*

0

*G*

**σ(s),****ε****u**(s)

*, β(s)*
*ds*−*k*

*n*
*j=1*

*G*
**σ***j−1**, ε*

**u**_{j−1}*, β*_{j−1}

(3.14)

and*I*is the identity operator on*Q. It can be verified that*(cf.[5])

1≤n≤Nmax*D**G,n*

*Q*≤*ck*

|* σ|*˙

*L*

^{∞}(0,T;Q)+|

**u**˙|

*L*

^{∞}(0,T;V)+|

*β|*˙

*L*

^{∞}(0,T;L

^{2}(Ω))

*.* (3.15)

Denote

*e*^{hk}* _{n}* =

**u**

*n*−

**u**

^{hk}

_{n}^{2}

*V*+**σ***n*−**σ**^{hk}_{n}^{2}

*Q*+*β**n*−*β*_{n}^{hk}^{2}

*L*^{2}(Ω)*,* *n*=0,1, . . . , N. (3.16)
From(3.13), we obtain

**σ***n*−**σ**^{hk}*n*

*Q*

≤*c*

* σ*0−

**σ**

^{h}_{0}

*+*

_{Q}*I*− P

*Q*

^{h}Eε

**u***n*−**u**0

*Q*+**u**0−**u**^{h}_{0}* _{V}*+

*D*

*G,n*

*+*

_{Q}*k*

*n*
*j=1*

*I*− P*Q*^{h}

*G*
**σ***j−1**,ε*

**u**_{j−1}*, β*_{j−1}

*Q*+**u*** _{n}*−

**u**

^{hk}

_{n}*V*

+^{n}

*j=1*

*k***u*** _{j}*−

**u**

^{hk}

_{j}*V*+**σ***j*−**σ**^{hk}_{j}

*Q*+*β**j*−*β*^{hk}_{j}

*L*^{2}(Ω)

*,*

(3.17) where properties(2.17)and(2.18)have been used.

Taking now**v**=**u**^{hk}* _{n}* in(3.11), we obtain

**σ***n*

*,*

**ε****u**^{hk}* _{n}* −

**u**

*n*

*Q*≥

**f***n**,u*^{hk}* _{n}* −

**u**

*n*

*V**.* (3.18)

Rewriting(3.7)in the following form:

**σ**^{hk}*n* *, ε*

**v*** _{n}*−

**u**

^{hk}

_{n}*Q*≥

**f**_{n}*,v** ^{h}*−

**u**

^{hk}

_{n}*V*+
**σ**^{hk}*n* *, ε*

**u*** _{n}*−

**v**

^{h}*Q**,* ∀**v*** ^{h}*∈

*U*

^{h}*,*(3.19) and subtracting the above two variational inequalities, we get

**σ***n*−**σ**^{hk}*n* *, ε*

**u*** _{n}*−

**u**

^{hk}

_{n}*Q*≤

**f**_{n}*,***u*** _{n}*−

**v**

^{h}*V*+
**σ***n**,ε*

**v*** ^{h}*−

**u**

_{n}*Q*

+

**σ**^{hk}*n* −**σ***n**, ε*

**v*** ^{h}*−

**u**

_{n}*Q**,* ∀**v*** ^{h}*∈

*U*

^{h}*.*(3.20)

Since
P*Q** ^{h}*Eε

**u***n*−**u**^{hk}_{n}*,ε*

**u***n*−**u**^{hk}_{n}

*Q*

= Eε

**u*** _{n}*−

**u**

^{hk}

_{n}*,ε*

**u*** _{n}*−

**u**

^{hk}

_{n}*Q*

+

P*Q** ^{h}*−

*I*Eε

**u***n*−**u**^{hk}_{n}*,ε*

**u***n*−**u**^{hk}_{n}

*Q*

= Eε

**u*** _{n}*−

**u**

^{hk}

_{n}*,ε*

**u*** _{n}*−

**u**

^{hk}

_{n}*Q*

+

P*Q** ^{h}*−

*I*Eε

**u***n*−**u**^{hk}_{n}*,ε*

**u***n*−**v**^{h}

*Q**,* ∀v* ^{h}*∈

*V*

^{h}*,*

(3.21)

introducing(3.10)and(3.6)into(3.20), we obtain Eε

**u***n*−**u**^{hk}_{n}*, ε*

**u***n*−**u**^{hk}_{n}

*Q*

≤

**f**_{n}*,***u*** _{n}*−

**v**

^{h}*V*+
**σ***n**, ε*

**v*** ^{h}*−

**u**

_{n}*Q*

−

**σ***n*−**σ**^{hk}_{n}*,ε*

**v*** ^{h}*−

**u**

*n*

*Q*+*c***u***n*−**u**^{hk}_{n}_{V}**u***n*−**v**^{h}_{V}

−
*I*− P*Q*^{h}

Eε

**u*** _{n}*−

**u**

_{0}

*,ε*

**u*** _{n}*−

**u**

^{hk}

_{n}*Q*−

* σ*0−

**σ**

^{h}_{0}

*,*

**ε****u*** _{n}*−

**u**

^{hk}

_{n}*Q*

−

*D**G,n*− P*Q** ^{h}*Eε

**u**0−**u**^{h}_{0}
*,ε*

**u***n*−**u**^{hk}_{n}

*Q*

+*k*
*n*
*j=1*

P*Q*^{h}

*G*
**σ***j−1**, ε*

**u**_{j−1}*, β*_{j−1}

−*G*
**σ**^{hk}_{j−1}*,ε*

**u**^{hk}_{j−1}*, β*_{j−1}^{hk}

*, ε*

**u*** _{n}*−

**u**

^{hk}

_{n}*Q*

+*k*
*n*
*j=1*

*I*− P*Q*^{h}

*G*
**σ***j−1**,ε*

**u**_{j−1}*, β*_{j−1}

*, ε*

**u*** _{n}*−

**u**

^{hk}

_{n}*Q**,* ∀**v*** ^{h}*∈

*U*

^{h}*.*(3.22) Thus, taking norms in the above inequality and using properties(2.17) and(2.18), we obtain the following estimate:

**u*** _{n}*−

**u**

^{hk}

_{n}^{2}

*V*

≤*c***u***n*−**v**^{h}* _{V}*+

**σ***n*−

**σ**

^{hk}

_{n}

_{Q}**u**

*n*−

**v**

^{h}*+*

_{V}**u**

*n*−

**u**

^{hk}

_{n}

_{V}**u**

*n*−

**v**

^{h}*+*

_{V}*0−*

**σ**

**σ**

^{h}_{0}

*Q***u*** _{n}*−

**u**

^{hk}

_{n}*V*+*D**G,n*

*Q***u*** _{n}*−

**u**

^{hk}

_{n}*V*

+**u**_{0}−**u**^{h}_{0}

*V***u*** _{n}*−

**u**

^{hk}

_{n}*V*

+

*k*
*n*
*j=1*

**u*** _{j}*−

**u**

^{hk}

_{j}*V*+**σ***j*−**σ**^{hk}_{j}

*Q*+*β**j*−β^{hk}_{j}

*L*^{2}(Ω)

**u*** _{n}*−

**u**

^{hk}

_{n}*V*

+*I*− P*Q*^{h}

Eε

**u*** _{n}*−

**u**

_{0}

*Q***u*** _{n}*−

**u**

^{hk}

_{n}*V*

+

*k*
*n*

*j=1*

*I*− P*Q*^{h}

*G*
**σ***j−1**, ε*

**u**_{j−1}*, β**j−1*

*Q*

**u*** _{n}*−

**u**

^{hk}

_{n}*V*

*.*
(3.23)

Then, applying the inequality

*ab*≤*δa*^{2}+ 1

4δ*b*^{2}*,* *δ, a, b*∈R*, δ >*0 (3.24)
to(3.23), we obtain

**u***n*−**u**^{hk}_{n}^{2}

*V* ≤*c***u***n*−**v**^{h}

*V*+*δ*0**σ***n*−**σ**^{hk}_{n}^{2}

*Q*+**u***n*−**v**^{h}^{2}

*V*+* σ*0−

**σ**

^{h}_{0}

^{2}

*Q*

+*D**G,n*^{2}

*Q*+**u**_{0}−**u**^{h}_{0}^{2}

*V*+*I*− P*Q*^{h}

Eε

**u*** _{n}*−

**u**

_{0}

^{2}

*Q*

+*k*^{2}
*n*
*j=1*

*e*^{hk}* _{j}* +

^{n}*j=1*

*k*^{2}*I*− P*Q*^{h}

*G*
**σ***j−1**, ε*

**u**_{j−1}*, β**j−1*^{2}

*Q*

*,*
(3.25)
where*δ*0is a constant parameter assumed to be small enough.

Now using(3.8)and(3.12)with*ξ*=*β*^{hk}* _{n}* and

*ξ*

*=*

^{h}*ξ*

_{n}*, we have*

^{h}*δ*

*β**n*−*β*_{n}^{hk}

*, β**n*−*β*^{hk}_{n}

*L*^{2}(Ω)

+*a*

*β**n*−*β*^{hk}_{n}*, β**n*−*β*^{hk}_{n}

≤

*δβ**n*−*β*˙*n**, β**n*−*β*_{n}^{hk}

*L*^{2}(Ω)

+
*δ*

*β**n*−*β*^{hk}_{n}

*, β**n*−*ξ*_{n}^{h}

*L*^{2}(Ω)−

*δβ**n**, β**n*−*ξ*^{h}_{n}

*L*^{2}(Ω)−*a*

*β**n**, β**n*−*ξ*_{n}* ^{h}*
+

*φ*
**σ***n**,ε*

**u**_{n}*, β**n*

*, β**n*−*ξ*_{n}^{h}

*L*^{2}(Ω)+*a*

*β**n*−*β*^{hk}_{n}*, β**n*−*ξ*^{h}* _{n}*
+

*φ*
**σ***n**,ε*

**u***n*
*, β**n*

−*φ*
**σ**^{hk}_{n−1}*,ε*

**u**^{hk}_{n−1}*, β*^{hk}_{n−1}

*, ξ*_{n}* ^{h}*−

*β*

_{n}

^{hk}*L*^{2}(Ω)*.*
(3.26)

Then, we bound the first term from below by

*δ*

*β**n*−*β*^{hk}_{n}

*, β**n*−*β*^{hk}_{n}

*L*^{2}(Ω)≥ 1
2k

*β**n*−*β*^{hk}_{n}^{2}

*L*^{2}(Ω)−*β** _{n−1}*−

*β*

_{n−1}

^{hk}^{2}

*L*^{2}(Ω)

*.*
(3.27)

Using this bound, replacing*n*by*j*, and making the summation over
*j*=1,2, . . . , n, after some algebraic manipulations, we obtain

*β**n*−*β*_{n}^{hk}^{2}_{L}_{2}_{(Ω)}+*k*
*n*
*j=1*

∇

*β**j*−*β*^{hk}_{j}^{2}_{[L}_{2}_{(Ω)]}_{d}

≤*c***u**0−**u**^{h}_{0}^{2}* _{V}*+

*0−*

**σ**

**σ**

^{h}_{0}

^{2}

*+*

_{Q}*β*0−

*β*

^{h}_{0}

^{2}

_{L}_{2}

_{(Ω)}+

*β*1−

*ξ*

^{h}_{1}

^{2}

_{L}_{2}

_{(Ω)}+

*k*

*n*
*j=1*

*β**j*−*β*^{hk}_{j}^{2}

*L*^{2}(Ω)+*k*
*n−1*

*j=1*

**u***j*−**u**^{hk}_{j}^{2}

*V*+**σ***j*−**σ**^{hk}_{j}^{2}

*Q*

+*k*^{2}+*k*
*n*
*j=1*

*δβ**j*−*β*˙*j*^{2}

*L*^{2}(Ω)+*k*
*n*

*j=1*

∇

*β**j*−*ξ*^{h}_{j}^{2}

[L^{2}(Ω)]^{d}

+*k*
*n*

*j=1*

*β**j*−*ξ*^{h}_{j}^{2}

*L*^{2}(Ω)+*β**n*−*ξ*^{h}_{n}^{2}

*L*^{2}(Ω)

+1
*k*

*n−1*

*j=1*

*β**j+1*−*ξ*^{h}_{j+1}

−

*β**j*−*ξ*_{j}^{h}^{2}

*L*^{2}(Ω)

+*k*
*n*

*j=1*

*φ*
**σ***j**,ε*

**u**_{j}*, β**j*

−*δβ**j*+*κ∆β**j*

*L*^{2}(Ω)·*β**j*−*ξ*_{j}^{h}

*L*^{2}(Ω)

*.*

(3.28)

We now combine estimates(3.17), (3.25), (3.28), and (3.15). Using in- equality(3.24), after some calculus, we obtain

**u*** _{n}*−

**u**

^{hk}

_{n}^{2}

*V*+**σ***n*−**σ**^{hk}*n* ^{2}

*Q*+*β**n*−*β*^{hk}_{n}^{2}

*L*^{2}(Ω)+*k*
*n*
*j=1*

∇

*β**j*−*β*_{j}^{hk}^{2}

[L^{2}(Ω)]^{d}

≤*c*

**u**_{0}−**u**^{h}_{0}^{2}

*V*+* σ*0−

**σ**

^{h}_{0}

^{2}

*Q*+*β*0−*β*_{0}^{h}^{2}

*L*^{2}(Ω)+*β*1−*ξ*_{1}^{h}^{2}

*L*^{2}(Ω)

+^{n}

*j=1*

*k*^{2}*e*^{hk}* _{j−1}*+

**u**

*n*−

**v**

^{h}^{2}

*+*

_{V}**u**

*n*−

**v**

^{h}*+*

_{V}*k*

^{2}+

*k*

+^{n}

*j=1*

*k*^{2}*I*− P*Q*^{h}

*G*
**σ***j−1**, ε*

**u**_{j−1}*, β**j−1*^{2}

*Q*+*k*
*n*

*j=1*

*δβ**j*−*β*˙*j*^{2}

*L*^{2}(Ω)

+*k*
*n*

*j=1*

∇

*β**j*−*ξ*^{h}_{j}^{2}

[L^{2}(Ω)]* ^{d}*+

*k*

*n*

*j=1*

*β**j*−*ξ*_{j}^{h}^{2}

*L*^{2}(Ω)+*β**n*−*ξ*_{n}^{h}^{2}

*L*^{2}(Ω)