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

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 provide a variational formulation for the mechanical problem and sketch a proof of the existence of a unique weak solution of the model. We then introduce and study a fully discrete scheme for the numerical solutions of the problem. An optimal order error estimate is derived for the approximate solutions under suitable solution regularity. Numerical examples are presented to show the performance of the method.


Introduction
We consider a mathematical model for a quasistatic process of frictionless 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 effect 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 mechanism of the system. Because of the importance of the safety issue of mechanical equipments, considerable effort has been devoted to modeling and numerically simulating damage.
Early models for mechanical damage derived from thermomechanical considerations appeared in [15,16], where numerical simulations were included. One-dimensional damage problems have been studied in [13,14]. Recently, the existence of weak solutions of viscoelastic problems with friction and damage have been provided in [20,27]. A quasistatic 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 material with constitutive relatioṅ σ = Eε(u) + G σ, ε(u), β , (1.1) where σ represents the stress tensor field, u denotes the displacement field, and ε(u) is the linearized strain tensor field. Here, E is a fourthorder tensor, 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 undamaged, 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 follows, the dot above a variable represents the time derivative. Note that for particular forms of the function G, the constitutive law (1.1) may describe 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 problems for such kind of materials were studied in [10,12,22,28]. A frictionless contact problem for materials of the form (1.1) in which β is an internal state variable whose evolution is described by an ordinary differential equation has been recently considered in [11].
One of the traits of novelty of this paper consists in the fact that, following [15,16], the evolution of the microscopic cracks responsible for the damage is modeled by the following differential inclusion: (1. 2) In (1.2) and below, κ > 0 is a constant, ∂ψ K denotes the subdifferential of the indicator function ψ K of K, which represents the set of admissible damage functions satisfying 0 ≤ β ≤ 1, and φ is a given constitutive function 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 unbounded when β → 0, a condition which is not allowed under our assumptions 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 results on the performance of the numerical method considered. Literature on the study of variational inequalities is rather extensive, see, for example, 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 Section 2, we present the mechanical problem and provide its variational analysis including an existence and uniqueness result, Theorem 2.3. The proof of Theorem 2.3 is based on classical results of elliptic and parabolic variational inequalities and Banach fixed point theorem. In Section 3, 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 difference to discretize the time derivative. We obtain an optimal order error estimate under appropriate regularity assumptions on the exact solution and data. Finally, in Section 4, we give some numerical examples to show the performance of the scheme.

Mechanical problem and variational formulation
The physical setting is as follows. A viscoplastic body occupies the domain Ω ⊂ 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 Γ 3 with 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.

2)
Div 3) Here, S d denotes the space of second-order symmetric tensors on R 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, measured 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 β 0 are 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 and R d , respectively, and we use the following spaces: (2.9) Here and below, i, j = 1, . . . , d, summation over repeated indices is implied, 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 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 normal 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: (2.12) Let V denote the closed subspace of H 1 defined by Since meas(Γ 1 ) > 0, Korn's inequality holds that there exists C K > 0 depending only on Ω and Γ 1 such that (2.14) A proof of Korn's inequality may be found in [25, page 79]. On V , we consider the inner product given by

15)
and let | · | V be the associated norm, that is, It follows that | · | H 1 and | · | 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 If X 1 and X 2 are real Hilbert spaces, then X 1 × X 2 denotes the product space which is endowed with the canonical inner product (·, ·) X 1 ×X 2 . Finally, 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 In the study of the mechanical Problem 2.1, we assume that the elas- (2.18) The damage source function φ : The body forces and surface tractions have the regularity The gap function g is such that

21)
and the initial data satisfy

27)
and let K denote the set of admissible damage functions By a standard procedure, we can derive the following variational formulation of the mechanical Problem 2.1.
The main steps of the proof are stated as follows.
(iii) Consider the Banach space X = L 2 (0, T; Q × L 2 (Ω)) with the norm | · | X given by and define the operator Λ : X → X by for any η ∈ X, t ∈ [0, T]. Then, the operator Λ has a unique fixed point

Numerical approximation
We analyze in this section a fully discrete approximation scheme for Problem 2.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 h ⊂ K be a nonempty, finite-dimensional closed convex set. Here h > 0 is a discretization parameter and we assume that ε(V h ) ⊂ Q h . This assumption is not a restriction for actual implementation of the method since, usually, V h and Q h are constructed to be finite-element spaces and V h consists of continuous piecewise polynomials of a degree one higher than that of Q h . Finally, denote by U h ⊂ V h an approximation for the convex set U for which we assume the following conformity condition: (3.1) Let P Q h : Q→Q h be the orthogonal projection operator defined through the relation

A contact problem in viscoplasticity with damage
The orthogonal projection operator is nonexpansive, that is, 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 discussion here to the case of nonuniform partition does not present any difficulty. For a continuous function z(t), we use the notation z n = z(t n ). For a sequence {z n } N n=0 , we denote δz n = (z n − z n−1 )/k for the corresponding divided difference. No summation is implied over the repeated index n.
In the rest of this section, c will denote positive constants which are independent of the discretization parameters h and k.
Let u h 0 ∈ U h , σ h 0 ∈ Q h , and β h 0 ∈ K h be chosen to approximate the initial values u 0 , σ 0 , and β 0 . A fully discrete approximation scheme for Problem 2.2 is the following problem.

4)
and for n = 1, 2, . . . , N, By induction, we obtain that this problem is equivalent to 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 h . To solve variational inequality (3.9) for u hk n and variational inequality (3.8) for β hk n , a penalty-duality algorithm can be used (see [2,3]). Once u hk n is known, we can determine σ 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 n ∈ K h . We summarize this result as follows. 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 Problem 2.2 (n = 1, . . . , N): β 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 and I is the identity operator on Q. It can be verified that (cf. [5]) From (3.13), we obtain

17)
where properties (2.17) and (2.18) have been used. Taking now v = u hk n in (3.11), we obtain Rewriting (3.7) in the following form:

19)
and subtracting the above two variational inequalities, we get (3.20)

A contact problem in viscoplasticity with damage
(3.23) Then, applying the inequality to (3.23), we obtain where δ 0 is a constant parameter assumed to be small enough. Now using (3.8) and Using this bound, replacing n by j, and making the summation over j = 1, 2, . . . , n, after some algebraic manipulations, we obtain (3.28) We now combine estimates (3.17), (3.25), (3.28), and (3.15). Using inequality (3.24), after some calculus, we obtain 102 A contact problem in viscoplasticity with damage (3.29) We will now use the following Gronwall's inequality (see [21]). (3.33) Inequality (3.32) is the basis for deriving error estimates. For instance, consider an approximation by using the finite-element method. We assume the following solution regularity: For simplicity, suppose that Ω is a polyhedral domain and let T h be a regular finite-element partition of the domain Ω, compatible with the boundary splitting Γ = Γ 1 ∪ Γ 2 ∪ Γ 3 . Let X h 1 and X h 0 be the corresponding finite-element spaces of continuous piecewise linear functions and piecewise constants, respectively. Then, we define We use the same symbol Π h for either the standard finite-element interpolation operator (see [6]) when the function to be interpolated is continuous, or the Clément's interpolation operator (see [7]) when the function to be interpolated is not continuous. We choose the initial values for the fully discrete scheme to be

37)
and moreover, we consider the following approximation for the convex subset U h : Here, g h is the piecewise Lagrange interpolation function of the gap function g, assumed to be concave in order to verify (3.1).
Then, using arguments similar to those in [20], we obtain the following result.
Theorem 3.5. Assume that the conditions stated in Theorem 2.3 hold and the solution regularity (3.34), (3.35), and (3.36). Then, choosing the initial data for the fully discrete scheme by (3.37), we have the following optimal order error estimate: Finally, we remark that error estimate (3.39) is only a sample result, obtained under the above regularity conditions. If the regularity conditions are different, the error estimate needs to be changed accordingly.

Numerical results
In order to verify the performance of the numerical method described in Section 3, some numerical experiences have been done in test examples including simulations in one and two dimensions. In this section, we resume some of these numerical simulations.

A one-dimensional test example
We consider a viscoplastic rod Ω = (0, L) clamped at its left end x = 0 and submitted to the action of a body force of density f 0 (x, t) in the xdirection (see where E > 0 is the Young's modulus of the material. We assume that the function G is a version of the Perzyna's law (see [9]), that is, where λ > 0 is a viscosity constant and P R(β) is the projection operator over the convex set R(β) defined as Here, σ Y is the uniaxial stress yield and β represents the damage function. Note that function (4.2) does not depend on the strain field and satisfies condition (2.18) (see, e.g., [23, page 92]). When σ ∈ R(β), G(σ, ε(u), β) = 0 and therefore (4.1) implies that only elastic deformations occur in the rod. When σ ∈ R(β), G(σ, ε(u), β) = 0 and therefore plastic deformations occur in the rod. Thus, the rod Ω = (0, L) is divided at each moment into two zones: the elastic zone (characterized by the condition σ ∈ R(β)) and the plastic zone (characterized by the condition σ ∈ R(β)). The boundaries of these zones are unknown a priori and depend on the damage field β. From (4.3), it follows that the elastic zones decrease with β since β 1 ≤ β 2 imply R(β 1 ) ⊂ R(β 2 ). This property models the effect of the damage of the material. In particular, if β = 0, then R(β) = {0} and plastic deformations occur for any stress σ = 0. We use in (2.2) the damage function which satisfies condition (2.19).

A contact problem in viscoplasticity with damage
A complete description of this problem is the following. (4.5) For the numerical computation, the following data have been considered: (4.6) The fully discrete scheme introduced in Section 3 was implemented by using continuous piecewise linear elements for the space V h and the set K h , and piecewise constant functions for the space Q h . We used the discretization parameters h = k = 0.01.
Our numerical results are plotted in

A two-dimensional test problem
As a two-dimensional example of Problem 2.1, we consider the physical situation described in Figure 4.5 during one second (i.e., T = 1second). We assume the plane stress hypothesis and therefore Ω = (0, 2) × (0, 2) is the cross section of a deformable three-dimensional body submitted to the action of vertical body forces (assumed to be linearly increasing in time). The body is clamped on Γ 1 = (0, 2) × {2} and it is in unilateral contact with a rigid foundation on Γ 3 = [0.8, 1.2] × {0}. Finally, we suppose that on Γ 2 = Γ − (Γ 1 ∪ Γ 3 ), the body is traction free.      We use a constitutive law of the form (2.1). The plane stress hypothesis allows us to define the elasticity tensor E by (Eτ) αβ = Eκ 1 − κ 2 τ 11 + τ 22 δ αβ +  where E is the Young's modulus and κ the Poisson's ratio of the material. We take G σ, ε(u), β = (β − 1)σ, (4.8) 112 A contact problem in viscoplasticity with damage which leads to the following model for the damage of the material: if β = 1 (undamaged material), the body is linearly elastic; if 0 ≤ β < 1 (partial damaged material), the body has a Maxwellian viscoelastic behavior and therefore irreversible strains occur. Here, the coefficient of relaxation r = 1 − β increases when β decreases. The maximum value r = 1 is obtained for β = 0, that is, in the case of completely damaged material. The damage function φ is defined as φ σ, ε(u), β = − 1 100 |σ| V M , (4.9) where | · | V M is the two-dimensional Von Mises norm for stress, |τ| 2 V M = τ 2 11 + τ 2 22 − τ 11 τ 22 + 3τ 2 12 , ∀τ ∈ S 2 . (4.10) The following data are used in the numerical test:  Again, we used continuous piecewise linear elements for the space V h and the set K h , and piecewise constant functions for the space Q h . Also, the value k = 0.01 for the time discretization parameter was considered.
The deformed mesh at final time and the initial boundary are shown in Figure 4.6. The Von Mises stress norm in the deformed configuration and the damage field at final time T = 1 second are plotted in Figure 4.7 (left side and right side, respectively).