In this paper, we study a mathematical model in the form of a hemivariational inequality for a quasistatic bilateral frictional contact problem between a viscoelastic body and a rigid foundation. The friction law is given in the form of subdifferential condition. Damage of the material is incorporated. Modeling, variational analysis and numerical solution of contact problems have been studied extensively; in this regard, a few comprehensive references are [22, 21, 16, 27] in the context of variational inequalities (VIs), and [23, 28] in the context of hemivariational inequalities (HVIs).
The notion of hemivariational inequalities (HVIs) was introduced in early 1980s to model mechanical problems involving non-smooth, non-monotone or multi-valued relations (). Early results on modeling, mathematical analysis and engineering applications of HVIs are summarized in [26, 24]; recent summarized accounts include [4, 23, 28]. Since there are no solution formulas for HVIs in applications, numerical simulation is the only feasible approach to solving HVIs. Detailed discussion of the finite element method for solving HVIs can be found in . More recently, there has been substantial progress in numerical analysis of HVIs, especially on optimal order error estimates for numerical solutions of HVIs, starting with the paper , followed by a sequence of papers, e.g., [2, 18, 13, 19]; the reader is referred to  for a recent survey.
Many contact processes are accompanied with material damage. In applications, it is very important to consider the damage effect. General mathematical models for damage were derived in [10, 11]; see also . In , a quasistatic contact problem for a viscoelastic material is studied variationally and numerically, where the damage effect of the viscoelastic material is taken into account. Systematic variational analysis and numerical analysis of contact problems with damage effect is summarized in . The mathematical problems investigated in these references are in the form of VIs. For studies of contact problems with damage in the form of HVIs, the reader is referred to .
This is the first paper devoted to numerical analysis of an HVI arising in a contact problem with damage. The rest of the paper is organized as follows. In Section 2, we introduce the contact problem, present its weak formulation as an HVI and comment on the solution existence and uniqueness. In Section 3, we consider a fully discrete numerical scheme for the contact problem and derive an optimal order error estimate under appropriate solution regularity assumptions. In Section 4, we report computer simulation results on a numerical example and illustrate numerical convergence orders that match the theoretical error bound.
2 The contact problem
We first introduce the pointwise formulation of the quasistatic contact problem for the contact between a viscoelastic body and a rigid foundation. The initial configuration of the body is , a Lipschitz bounded domain in ( in applications). The body is subject to the action of volume forces of a total density . The boundary of the domain is split into three disjoint measurable parts, such that is non-trivial. We will assume the body is fixed along , is subject to the action of surface tractions with a total density . Along the contact boundary , the body and the foundation are in bilateral contact and the frictional process is described by a generalized subdifferential inclusion.
where is the displacement field, is the damage function, is the stress field, and are the viscosity operator and the elasticity operator. These operators are allowed to depend on the spatial location. For convenience, we use the shorthand notation and for and , respectively. The symbol denotes the time derivative of . The time interval of interest is for some .
The notion of the damage function was introduced in [10, 11] to quantify the damage to the material. It is defined to be the ratio between the elastic modulus of the damaged material and that of the original material. The value of the damage function lies in . The value indicates that there is no damage in the material, whereas the value corresponds to a completely damaged material. When , there is a partial damage and the system has a reduced load carrying capacity. A popular model for the evolution of the damage function is given by a parabolic differential inclusion:
where is a constant microcrack diffusion coefficient, is the indicator function of the interval , is the convex subdifferential of , and is the mechanical source of damage, depending on the strain and the damage itself. On the boundary , a homogeneous Neumann condition is described for .
For a vectordefined on , we let be its normal component, and let
be its tangential component. For a stress tensordefined on , we let and be its normal and tangential components, respectively.
Denote by and the initial values of the displacement and the damage function. The pointwise formulation of the contact problem is as follows.
Find a displacement field , a stress field , and a damage field such that
We already know that (2.1) is the viscoelastic constitutive law with damage, and (2.2) is the evolution relation for the damage function. We consider a quasistatic contact process and (2.3) is the corresponding equilibrium equation. The initial conditions for the displacement field and the damage function are given by (2.8). The relations (2.4)–(2.7) are the boundary condition for the damage function, the displacement boundary condition on , the traction boundary condition on , and the bilateral friction contact condition on . Here, the friction dissipation pseudopotential will be assumed to be Lipschitz continuous, and represents the generalized subdifferential in the sense of Clarke (cf. [6, 7]). We will also need the notion of the generalized directional derivative in the sense of Clarke. Let be a Banach space and let be a locally Lipschitz continuous functional. Recall that the generalized directional derivative of at in the direction is
whereas the generalized subdifferential of at is
We note the following properties:
Problem 2.1 will be studied in its weak form. For this purpose, we first need to introduce some function spaces. Let
which is a Hilbert space with the inner product
This will be the space for stress and strain fields. The function space for the displacement field is the Hilbert space
with the inner product and the associated norm . The space for the damage field is . For convenience, we let . The spaces and are endowed with their canonical inner products and norms.
In the study of the contact problem, we assume that the operator satisfies the following conditions:
Similarly, we assume the operator has the following properties:
As an example of the viscoelastic constitutive law with damage, we consider
where the viscosity tensor satisfies (2.13), is a positive coefficient, is a damage dependent elasticity set, which is assumed to be convex and is the projection operator onto the set . We require the properties and implies . The second property implies that as the damage of the material increases, i.e., the value of the damage function decreases, the elasticity convex set expands, and the material resembles a purely viscous one. A concrete example is given by the von Mises convex set
where is the deviatoric part of , and is the yield limit of the damage-free material. Since the projection operator is a contraction, it can be verified that satisfies (2.14).
On the damage source function , the assumptions are
On the friction dissipation pseudopotential ,
Moreover, we assume
on the microcrack diffusion coefficient,
on the densities of forces and tractions,
on the initial data. Here represents the set of admissible damage functions defined by
By the Riesz representation theorem, we can define by
Then conditions (2.20) imply
Let be the bilinear form
Find a displacement field , a stress field , and a damage field such that for all ,
Denote by the smallest constant in the trace inequality
Similar to [12, Theorem 5.1], we can prove the following result.
3 Numerical analysis of the weak formulation
For the approximation of the time derivative of the damage function, we use finite difference. We divide the time interval uniformly and comment that much of the discussion of the numerical method below can be extended straightforward to the case of general partition of the time interval. Thus, let be a positive integer, and define the step-size. Then is a uniform partition of with the nodes , . For a function continuous on , we write . We use the backward difference approximation
For the spatial discretization, we use the finite element method. For simplicity, we assume is a polygonal/polyhedral domain, and express the three parts of the boundary, , , as unions of closed flat components with disjoint interiors:
Let be a regular family of finite element partitions of into triangular/tetrahedral elements, compatible with the partition of the boundary into , , , in the sense that if the intersection of one side/face of an element with one set has a positive measure with respect to , then the side/face lies entirely in . Here denotes the finite element mesh-size. Corresponding to the partition , we introduce the linear finite element space
for the displacement field, the piecewise constant finite element space
for the stress field, and the linear finite element space
for the damage field. Define the constrained subset of :
Let and be appropriate approximations of and such that
These conditions are valid if, e.g., , , and we define
to be the interpolant or- or -projection of onto , define to be the -projection of onto . The smoothness conditions and will follow from the solution regularities (3.14) and (3.16) below.
The discrete velocity and displacement approximations are denoted by and , whereas the discrete stress and damage function approximations are denoted by and . Let be the orthogonal projection from to , defined by
Then a fully discrete scheme for Problem 2.2 is the following.
Find a discrete displacement field , a discrete stress field , and a discrete damage field such that for ,
Here and are related by the equalities
The solution existence and uniqueness of Problem 3.1 can be proved by an induction argument. The focus of the rest of this section is to bound the numerical solution errors. For this purpose, we assume the following additional solution regularities: