# Numerical Studies of a Hemivariational Inequality for a Viscoelastic Contact Problem with Damage

This paper is devoted to the study of a hemivariational inequality modeling the quasistatic bilateral frictional contact between a viscoelastic body and a rigid foundation. The damage effect is built into the model through a parabolic differential inclusion for the damage function. A solution existence and uniqueness result is presented. A fully discrete scheme is introduced with the time derivative of the damage function approximated by the backward finite different and the spatial derivatives approximated by finite elements. An optimal order error estimate is derived for the fully discrete scheme when linear elements are used for the velocity and displacement variables, and piecewise constants are used for the damage function. Simulation results on numerical examples are reported illustrating the performance of the fully discrete scheme and the theoretically predicted convergence orders.

## Authors

• 6 publications
• 3 publications
• 3 publications
03/27/2020

### Solving Quasistatic Contact Problems Using Nonsmooth Optimization Approach

This paper is devoted to a study of time-dependent hemivariational inequ...
02/21/2022

### Semi-discrete and fully discrete weak Galerkin finite element methods for a quasistatic Maxwell viscoelastic model

This paper considers weak Galerkin finite element approximations for a q...
01/22/2021

### Mixed Finite Element Discretization for Maxwell Viscoelastic Model of Wave Propagation

This paper considers semi-discrete and fully discrete mixed finite eleme...
03/07/2022

### A DG/CR discretization for the variational phase-field approach to fracture

Variational phase-field models of fracture are widely used to simulate n...
04/03/2020

### Finite difference approach to fourth-order linear boundary-value problems

Discrete approximations to the equation L_contu = u^(4) + D(x) u^(3) ...
02/11/2022

### Identifying Processes Governing Damage Evolution in Quasi-Static Elasticity. Part 2 – Numerical Simulations

We investigate numerically a quasi-static elasticity system of Kachanov-...
09/04/2019

### The polarization process of ferroelectric materials analyzed in the framework of variational inequalities

We are concerned with the mathematical modeling of the polarization proc...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

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 ([25]). 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 [20]. 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 [14], followed by a sequence of papers, e.g., [2, 18, 13, 19]; the reader is referred to [17] 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 [9]. In [15], 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 [27]. 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 [12].

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.

Following [8, 16], we consider a viscoelastic constitutive law with the damage effect in the form

 \boldmath{σ}=A\boldmath{ε}(˙% \boldmath{u})+B(\boldmath{ε}(\boldmath{u% }),ζ),

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:

 ˙ζ−κ△ζ+∂I[0,1](ζ)∋ϕ(% \boldmath{ε}(\boldmath{u}),ζ),

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 vector

defined on , we let be its normal component, and let

be its tangential component. For a stress tensor

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

###### Problem 2.1

Find a displacement field , a stress field , and a damage field such that

 σ =A\boldmath{ε}(˙\boldmath{u% })+B(\boldmath{ε}(\boldmath{u}),ζ) in Ω×(0,T), (2.1) ˙ζ−κ△ζ+∂I[0,1](ζ) ∋ϕ(\boldmath{ε}(\boldmath{u}),ζ) in Ω×(0,T), (2.2) Div\boldmath{σ}+\boldmath{f}0 =\boldmath{0} in Ω×(0,T), (2.3) ∂ζ∂ν =0 on Γ×(0,T), (2.4) u =\boldmath{0} on Γ1×(0,T), (2.5) σν =\boldmath{f}2 on Γ2×(0,T), (2.6) uν =0, −\boldmath{σ}τ∈∂j(˙% \boldmath{u}τ) on Γ3×(0,T), (2.7) \boldmath{u}(0) =\boldmath{u}0, ζ(0)=ζ0 in Ω. (2.8)

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

 ψ0(u;v):=limsupw→u,λ↓0ψ(w+λv)−ψ(w)λ,

whereas the generalized subdifferential of at is

 ∂ψ(u):={ξ∈V∗∣ψ0(u;v)≥⟨ξ,v⟩V∗×V ∀v∈V}.

We note the following properties:

 ψ0(u;tv)=tψ0(u;v)∀u,v∈V,t≥0, (2.9) ψ0(u;v1+v2)≤ψ0(u;v1)+ψ0(u;v2)∀u,v1,v2∈V, (2.10) ψ0(u;v)=max{⟨ζ,v⟩V∗×V∣ζ∈∂ψ(u)}∀u,v∈V, (2.11) un→u and vn→v in V⟹limsupn→∞ψ0(un;vn)≤ψ0(u;v). (2.12)

Problem 2.1 will be studied in its weak form. For this purpose, we first need to introduce some function spaces. Let

 Q=L2(Ω)d×dsym,

which is a Hilbert space with the inner product

 (\boldmath{σ},\boldmath{τ})Q=∫Ωσijτijdx,\boldmath{σ},\boldmath{τ}∈Q.

This will be the space for stress and strain fields. The function space for the displacement field is the Hilbert space

 V={\boldmath{v}∈H1(Ω)d∣\boldmath{v}=\boldmath{0} on Γ1, vν=0 on Γ3}

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:

 (a) There exists LA>0 such that∥A(\boldmath{x},\boldmath{ε}1)−A(\boldmath{x},\boldmath{ε}2)∥≤LA∥\boldmath{ε}1−\boldmath{ε}2∥∀\boldmath{ε}1,\boldmath{ε}2∈Sd, a.e. \boldmath{x}∈Ω.(b) There exists mA>0 such that(A(\boldmath{x},\boldmath{ε}1)−A(\boldmath{x},\boldmath{ε}2))⋅(%\boldmath$ε$1−\boldmath{ε}2)≥mA∥\boldmath{ε}1−\boldmath{ε}2∥2∀\boldmath{ε}1,\boldmath{ε}2∈Sd, a.e. \boldmath{x}∈Ω.(c) For any \boldmath{ε}∈Sd, % \boldmath{x}↦A(\boldmath{x},\boldmath{ε}) is measurable on Ω.(d) The mapping \boldmath{x}↦A(\boldmath{x},\boldmath{0}) belongs to Q.⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭ (2.13)

Similarly, we assume the operator has the following properties:

 (a) There exists LB>0 such that∥B(\boldmath{x},\boldmath{ε}1,ζ1)−B(\boldmath{x},\boldmath{ε}2,ζ2)∥≤LB(∥\boldmath{ε}1−% \boldmath{ε}2∥+|ζ1−ζ2|)∀\boldmath{ε}1,\boldmath{ε}2∈Sd, ζ1,ζ2∈R, a.e. \boldmath{x}∈Ω.(b) For any \boldmath{ε}∈Sd and ζ∈R, \boldmath{x}↦B(\boldmath{x},\boldmath{ε},ζ) is measurable on Ω.(c) The mapping \boldmath{x}↦B(\boldmath{x},\boldmath{0},0) belongs to Q.⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭ (2.14)

As an example of the viscoelastic constitutive law with damage, we consider

 \boldmath{σ}=A\boldmath{ε}(˙% \boldmath{u})+η(\boldmath{ε}(\boldmath{u})−PK(ζ)(\boldmath{ε}(\boldmath{u}))), (2.15)

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

 K(ζ)={\boldmath{τ}∈Sd∣∥\boldmath% {τ}D∥≤ζσY}, (2.16)

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

 (a) There exists Lϕ>0 such that|ϕ(\boldmath{x},\boldmath{ε}1,ζ1)−ϕ(\boldmath{x},\boldmath{ε}2,ζ2)|≤Lϕ(∥\boldmath{ε}1−\boldmath{ε}2∥+|ζ1−ζ2|)∀\boldmath{ε}1,\boldmath{ε}2∈Sd, ζ1,ζ2∈R, a.e. \boldmath{x}∈Ω.(b) For any \boldmath{ε}∈Sd and ζ∈R, \boldmath{x}↦ϕ(\boldmath{x},\boldmath{ε},ζ) is measurable on Ω.(c) The mapping \boldmath{x}↦ϕ(\boldmath{x},\boldmath{0},0) belongs to L2(Ω).⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭ (2.17)

On the friction dissipation pseudopotential ,

 (a) \boldmath{x}↦j(% \boldmath{x},\boldmath{ξ}) is measurable on Γ3 ∀\boldmath{ξ}∈Rdand j(\boldmath{x},\boldmath{0})∈L2(Γ3).(b) \boldmath{ξ}↦j(\boldmath{x},% \boldmath{ξ}) is locally Lipschitz in Rd,a.e. \boldmath{x}∈Γ3.(c) There exist constants c0τ,c1τ≥0 such that∥∂j(\boldmath{x},\boldmath{ξ})∥≤c0τ+c1τ∥\boldmath{ξ}∥ ∀\boldmath{ξ}∈Rd,a.e. \boldmath{x}∈Γ3.(d) There is a constant c2τ≥0 such thatj0(\boldmath{x},\boldmath{ξ}1;% \boldmath{ξ}2−\boldmath{ξ}1)+j0(\boldmath{x},\boldmath{ξ}2;\boldmath{ξ}1−\boldmath{ξ% }2)≤c2τ∥\boldmath{ξ}1−\boldmath{ξ}2∥∀\boldmath{ξ}1,\boldmath{ξ}2∈Rd,a.e. \boldmath{x}∈Γ3.⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎬⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎭ (2.18)

Moreover, we assume

 κ>0 (2.19)

on the microcrack diffusion coefficient,

 \boldmath{f}0∈C([0,T];L2(Ω)d),% \boldmath{f}2∈C([0,T];L2(Γ2)d) (2.20)

on the densities of forces and tractions,

 \boldmath{u}0∈V,ζ0∈K (2.21)

on the initial data. Here represents the set of admissible damage functions defined by

 K={ξ∈Z∣ξ∈[0,1] a.e.\ in Ω}. (2.22)

By the Riesz representation theorem, we can define by

 (\boldmath{f}(t),\boldmath{v})V=∫Ω% \boldmath{f}0(t)⋅\boldmath{v}dx+∫Γ2% \boldmath{f}2(t)⋅\boldmath{v}da∀% \boldmath{v}∈V, t∈[0,T]. (2.23)

Then conditions (2.20) imply

 \boldmath{f}∈C([0,T];V). (2.24)

Let be the bilinear form

 a(ξ,η)=κ∫Ω∇ξ⋅∇ηdx,ξ,η∈Z. (2.25)

Let us introduce a weak formulation of the problem (2.1)–(2.8).

###### Problem 2.2

Find a displacement field , a stress field , and a damage field such that for all ,

 \boldmath{σ}(t)=A\boldmath{ε}(˙\boldmath{u}(t))+B(\boldmath{ε}(\boldmath{u}(t)),ζ(t)), (2.26) (\boldmath{σ}(t),\boldmath{ε}(\boldmath{v}−˙\boldmath{u}(t)))Q+∫Γ3j0(˙\boldmath{u}τ(t);\boldmath{v}τ−˙\boldmath{u}τ(t))da≥(\boldmath{f}(t),% \boldmath{v}−˙\boldmath{u}(t))V∀% \boldmath{v}∈V, (2.27)
 ζ(t)∈K, (˙ζ(t),ξ−ζ(t))Z0+a(ζ(t),ξ−ζ(t)) ≥(ϕ(\boldmath{ε}(\boldmath{u}(t)),ζ(t)),ξ−ζ(t))Z0∀ξ∈K, (2.28)

and

 \boldmath{u}(0)=\boldmath{u}0,ζ(0)=ζ0. (2.29)

Denote by the smallest constant in the trace inequality

 ∥\boldmath{v}τ∥L2(Γ3)d≤cτ∥% \boldmath{v}∥V∀\boldmath{v}∈V. (2.30)

Similar to [12, Theorem 5.1], we can prove the following result.

###### Theorem 2.3

Assume (2.13)–(2.21),

 c2τc2τ

and either or for a constant , for all and a.e. . Then Problem 2.2 has a unique solution , and . Moreover, .

In the numerical solution of Problem 2.2, it will be convenient to introduce the velocity variable

 \boldmath{w}(t)=˙\boldmath{u}(t). (2.32)

Given the velocity and the initial displacement from (2.29), we can recover the displacement by the formula

 \boldmath{u}(t)=\boldmath{u}0+∫t0\boldmath% {w}(s)ds. (2.33)

Then (2.26)–(2.27) can be rewritten as

 \boldmath{σ}(t)=A\boldmath{ε}(\boldmath{w}(t))+B(\boldmath{ε}(\boldmath{u}(t)),ζ(t)), (2.34) (\boldmath{σ}(t),\boldmath{ε}(\boldmath{v}−\boldmath{w}(t)))Q+∫Γ3j0(\boldmath{w}τ(t);\boldmath{v}τ−\boldmath{w}τ(t))da≥(\boldmath{f}(t),\boldmath{v}−% \boldmath{w}(t))V∀\boldmath{v}∈V, (2.35)

for all .

## 3 Numerical analysis of the weak formulation

In this section, we introduce and study a fully discrete numerical scheme to solve Problem 2.2. We assume the conditions stated in Theorem 2.3 are valid so that Problem 2.2 has a unique solution.

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

 ˙ζ(tn)≈δζn:=ζn−ζn−1k,1≤n≤N. (3.1)

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:

 ¯¯¯¯¯¯Γk=∪iki=1Γk,i,1≤k≤3.

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

 Vh={\boldmath{v}h∈C(¯¯¯¯Ω)d∣% \boldmath{v}h|T∈P1(T)d  ∀ T∈Th,\boldmath{v}h=\boldmath{0} on Γ1,vhν=0 on Γ3} (3.2)

for the displacement field, the piecewise constant finite element space

 Qh={\boldmath{τ}h∈Q∣\boldmath{τ}h|T∈Rd×d ∀T∈Th} (3.3)

for the stress field, and the linear finite element space

 Zh={ξh∈C(¯¯¯¯Ω)∣ξh|T∈P1(T) ∀T∈Th} (3.4)

for the damage field. Define the constrained subset of :

 Kh={ξh∈Zh∣ξh|T∈[0,1] ∀T∈Th}. (3.5)

Let and be appropriate approximations of and such that

 ∥\boldmath{u}0−\boldmath{u}h0∥V≤ch,∥ζ0−ζh0∥Z0≤ch. (3.6)

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

 PQh\boldmath{σ}∈Qh,(PQh\boldmath{σ},\boldmath{τ}h)Q=(\boldmath{% σ},\boldmath{τ}h)Q∀\boldmath{σ}∈Q,\boldmath{τ}h∈Qh. (3.7)

Then a fully discrete scheme for Problem 2.2 is the following.

###### Problem 3.1

Find a discrete displacement field , a discrete stress field , and a discrete damage field such that for ,

 \boldmath{σ}hkn=PQhA%\boldmath$ε$(\boldmath{w}hkn)+PQhB(\boldmath{ε}(\boldmath{u}hkn−1),ζhkn−1), (3.8) (\boldmath{σ}hkn,\boldmath{ε}(\boldmath{v}h−\boldmath{w}hkn))Q+∫Γ3j0(\boldmath{w}hknτ;\boldmath{v}hτ−\boldmath{w}hknτ)da≥(\boldmath{f}n,\boldmath{v}h−\boldmath{w}hkn)V∀\boldmath{v}h∈Vh, (3.9) (δζhkn,ξh−ζhkn)Z0+a(ζhkn,ξh−ζhkn)≥(ϕ(\boldmath{ε}(% \boldmath{u}hkn−1),ζhkn−1),ξh−ζhkn)Z0∀ξh∈Kh, (3.10)

and

 \boldmath{u}hk0=\boldmath{u}h0,ζhk0=ζh0. (3.11)

Here and are related by the equalities

 \boldmath{w}hkn=δ\boldmath{u}hkn% and\boldmath{u}hkn=\boldmath{u}h0+kn∑i=1\boldmath{w}hki. (3.12)

Note that for implementation, (3.8) and (3.9) are combined together to give

 (A\boldmath{ε}(\boldmath{w}hkn),\boldmath{ε}(\boldmath{v}h−% \boldmath{w}hkn))Q+(B(\boldmath{ε}(%\boldmath$u$hkn−1),ζhkn−1),\boldmath{ε}(\boldmath{v}h−\boldmath{w}hkn))Q+∫Γ3j0(\boldmath{w}hknτ;\boldmath{v}hτ−\boldmath{w}hknτ)da ≥(\boldmath{f}n,\boldmath{v}h−\boldmath{w}hkn)V∀\boldmath{v}h∈Vh. (3.13)

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:

 \boldmath{u}∈W2,1(0,T;V)∩C(